Merge pull request #133 from gridap/dirac_delta

Support for Dirac delta defined by tags

src/CellData.jl

@@ -346,12 +346,12 @@ function Gridap.Arrays.evaluate!(cache,s::DistributedCellDof,f::DistributedCellF
 end
 
 function Gridap.Arrays.evaluate!(cache, ::DistributedCellField, ::DistributedCellDof)
-  @unreachable """\n
-  CellField (f) objects cannot be evaluated at CellDof (s) objects.
-  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)?
-  """
+@unreachable """\n
+CellField (f) objects cannot be evaluated at CellDof (s) objects.
+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
 
 # Interpolation at arbitrary points (returns -Inf if the point is not found)
@@ -423,3 +423,63 @@ function Gridap.Arrays.evaluate!(caches,I::DistributedInterpolable,v::AbstractVe
     end
   # reduce((v,w)->broadcast(max,v,w),y)
 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

test/CellDataTests.jl

@@ -58,7 +58,7 @@ function main(distribute,parts)
   u3 = CellField(2.0,Ω)
   u = _my_op∘(u1,u2,u3)
 
-  order = 1
+ order = 1
   reffe = ReferenceFE(lagrangian,Float64,order)
   V = TestFESpace(model,reffe)
   uh = interpolate_everywhere(x->x[1]+x[2],V)
@@ -77,6 +77,19 @@ function main(distribute,parts)
     @test v ==[0.2,1.0,1.8]
   end
 
+  # 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