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