From fa3584e3ccce09c413374b808200a1ba13060c0e Mon Sep 17 00:00:00 2001
From: principe <principe@cimne.upc.edu>
Date: Thu, 2 Nov 2023 07:52:17 +0100
Subject: [PATCH] Support for Dirac delta defined by tags

---
 src/CellData.jl       | 62 ++++++++++++++++++++++++++++++++++++++++++-
 test/CellDataTests.jl | 13 +++++++++
 2 files changed, 74 insertions(+), 1 deletion(-)

diff --git a/src/CellData.jl b/src/CellData.jl
index e0f93c99..78bd9b95 100644
--- a/src/CellData.jl
+++ b/src/CellData.jl
@@ -357,4 +357,64 @@ However, CellDofs objects can be evaluated at CellField objects.
 Did you mean evaluate(f,s) instead of evaluate(s,f), i.e.
 f(s) instead of s(f)?
 """
-end
\ No newline at end of file
+end
+
+# Support for distributed Dirac deltas
+struct DistributedDiracDelta{D} <: GridapType
+  Γ::DistributedTriangulation
+  dΓ::DistributedMeasure
+end
+# This code is from Gridap, repeated in BoundaryTriangulations.jl and DiracDelta.jl.
+# We could refactor...
+import Gridap.Geometry: FaceToCellGlue
+function BoundaryTriangulation{D}(
+  model::DiscreteModel,
+  face_to_bgface::AbstractVector{<:Integer}) where D
+
+  bgface_to_lcell = Fill(1,num_faces(model,D))
+
+  topo = get_grid_topology(model)
+  bgface_grid = Grid(ReferenceFE{D},model)
+  face_grid = view(bgface_grid,face_to_bgface)
+  cell_grid = get_grid(model)
+  glue = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell)
+  trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)
+  BoundaryTriangulation(trian,glue)
+end
+
+function BoundaryTriangulation{D}(model::DiscreteModel;tags) where D
+  labeling = get_face_labeling(model)
+  bgface_to_mask = get_face_mask(labeling,tags,D)
+  face_to_bgface = findall(bgface_to_mask)
+  BoundaryTriangulation{D}(model,face_to_bgface)
+end
+
+function DiracDelta{D}(model::DistributedDiscreteModel{Dc},degree::Integer;kwargs...) where {D,Dc}
+
+  @assert 0 <= D && D < num_cell_dims(model) """\n
+  Incorrect value of D=$D for building a DiracDelta{D} on a model with $(num_cell_dims(model)) cell dims.
+
+  D should be in [0,$(num_cell_dims(model))).
+  """
+
+  gids   = get_face_gids(model,Dc)
+  trians = map(local_views(model),partition(gids)) do model, gids
+      trian = BoundaryTriangulation{D}(model;kwargs...)
+      filter_cells_when_needed(no_ghost,gids,trian)
+  end
+  Γ=DistributedTriangulation(trians,model)
+  dΓ=Measure(Γ,degree)
+  DistributedDiracDelta{D}(Γ,dΓ)
+end
+
+# Following functions can be eliminated introducing an abstract delta in Gridap.jl
+function DiracDelta{0}(model::DistributedDiscreteModel;tags)
+  degree = 0
+  DiracDelta{0}(model,degree;tags=tags)
+end
+function (d::DistributedDiracDelta)(f)
+ evaluate(d,f)
+end
+function Gridap.Arrays.evaluate!(cache,d::DistributedDiracDelta,f)
+ ∫(f)*d.dΓ
+end
diff --git a/test/CellDataTests.jl b/test/CellDataTests.jl
index 3bc15911..cf22202b 100644
--- a/test/CellDataTests.jl
+++ b/test/CellDataTests.jl
@@ -58,6 +58,19 @@ function main(distribute,parts)
   u3 = CellField(2.0,Ω)
   u = _my_op∘(u1,u2,u3)
 
+  # Point δ
+  δ=DiracDelta{0}(model;tags=2)
+  @test sum(δ(f)) ≈ 4.0
+  @test sum(δ(3.0)) ≈ 3.0
+  @test sum(δ(x->2*x)) ≈ VectorValue(8.0,0.0)
+  
+  # Line δ
+  degree = 2
+  δ = DiracDelta{1}(model,degree,tags=5)
+  @test sum(δ(f)) ≈ 8.0
+  @test sum(δ(3.0)) ≈ 12.0
+  @test sum(δ(x->2*x)) ≈ VectorValue(16.0,0.0)
+
 end
 
 end # module