Adding support to DiracDelta for generic points in the domain

NEWS.md

@@ -10,6 +10,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `lastindex` for `MultiValue`s for consistent usage of `[end]` as per `length`. Since PR [#834](https://github.com/gridap/Gridap.jl/pull/834)
 - BDM (Brezzi-Douglas-Marini) ReferenceFEs in PR [#823](https://github.com/gridap/Gridap.jl/pull/823)
 - Implemented `ConstantFESpace`. This space allows, e.g., to impose the mean value of the pressure via an additional unknown/equation (a Lagrange multiplier). Since PR [#836](https://github.com/gridap/Gridap.jl/pull/836)
+- Methods to construct `DiracDelta` at generic `Point`(s) in the domain. Since PR [#837](https://github.com/gridap/Gridap.jl/pull/837)
+- some key missing `lastindex` and `end` tests belonging to PR [#834]. Since PR [#837](https://github.com/gridap/Gridap.jl/pull/837)
 
 ## [0.17.14] - 2022-07-29
 

src/CellData/DiracDeltas.jl

@@ -1,9 +1,14 @@
+abstract type DiracDeltaSupportType <: GridapType end
+abstract type IsGridEntity <: DiracDeltaSupportType end
+abstract type NotGridEntity <: DiracDeltaSupportType end
 
-struct DiracDelta{D} <: GridapType
-  Γ::Triangulation{D}
+struct GenericDiracDelta{D,Dt,S<:DiracDeltaSupportType} <: GridapType
+  Γ::Triangulation{Dt}
   dΓ::Measure
 end
 
+const DiracDelta{D} = GenericDiracDelta{D,D,IsGridEntity}
+
 function DiracDelta{D}(
   model::DiscreteModel,
   face_to_bgface::AbstractVector{<:Integer},
@@ -24,7 +29,7 @@ function DiracDelta{D}(
   trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)
   Γ = BoundaryTriangulation(trian,glue)
   dΓ = Measure(Γ,degree)
-  DiracDelta(Γ,dΓ)
+  DiracDelta{D}(Γ,dΓ)
 end
 
 function DiracDelta{D}(
@@ -65,14 +70,65 @@ function DiracDelta{0}(model::DiscreteModel;tags)
   DiracDelta{0}(model,degree;tags=tags)
 end
 
-function (d::DiracDelta)(f)
+function (d::GenericDiracDelta)(f)
   evaluate(d,f)
 end
 
-function evaluate!(cache,d::DiracDelta,f)
+function evaluate!(cache,d::GenericDiracDelta,f)
   ∫(f)*d.dΓ
 end
 
-function get_triangulation(d::DiracDelta)
+function get_triangulation(d::GenericDiracDelta)
   d.Γ
 end
+
+# For handling DiracDelta at a generic Point in the domain #
+
+function _cell_to_pindices(pvec::Vector{<:Point},trian::Triangulation)
+  cache = _point_to_cell_cache(KDTreeSearch(),trian)
+  cell_to_pindex = Dict{Int, Vector{Int32}}()
+  for i in 1:length(pvec)
+    cell = _point_to_cell!(cache, pvec[i])
+    push!(get!(() -> valtype(cell_to_pindex)[], cell_to_pindex, cell), i)
+  end
+  cell_to_pindex
+end
+
+function DiracDelta(model::DiscreteModel{D}, p::Point{D,T}) where {D,T}
+  trian = Triangulation(model)
+  cache = _point_to_cell_cache(KDTreeSearch(),trian)
+  cell = _point_to_cell!(cache, p)
+  trianv = TriangulationView(trian,[cell])
+  point = [p]
+  weight = [one(T)]
+  pquad = GenericQuadrature(point,weight)
+  pmeas = Measure(CellQuadrature([pquad],[pquad.coordinates],[pquad.weights],trianv,PhysicalDomain(),PhysicalDomain()))
+  GenericDiracDelta{0,D,NotGridEntity}(trianv,pmeas)
+end
+
+function DiracDelta(model::DiscreteModel{D}, pvec::Vector{Point{D,T}}) where {D,T}
+  trian = Triangulation(model)
+  cell_to_pindices = _cell_to_pindices(pvec,trian)
+  cell_ids = collect(keys(cell_to_pindices))
+  cell_points = collect(values(cell_to_pindices))
+  points = map(i->pvec[cell_points[i]], 1:length(cell_ids))
+  weights_x_cell = collect.(Fill.(one(T),length.(cell_points)))
+  pquad = map(i -> GenericQuadrature(points[i],weights_x_cell[i]), 1:length(cell_ids))
+  trianv = TriangulationView(trian,cell_ids)
+  pmeas = Measure(CellQuadrature(pquad,points,weights_x_cell,trianv,PhysicalDomain(),PhysicalDomain()))
+  GenericDiracDelta{0,D,NotGridEntity}(trianv,pmeas)
+end
+
+function evaluate!(cache,d::GenericDiracDelta{0,Dt,NotGridEntity},f::Function) where Dt
+  dc = DomainContribution()
+  quad_points = d.dΓ.quad.cell_point
+  weights = d.dΓ.quad.cell_weight
+  dc_Γ = zeros(eltype(eltype(weights)), length(weights))
+  for i in eachindex(weights)
+    for j in eachindex(weights[i])
+      @inbounds dc_Γ[i] += f(quad_points[i][j])*weights[i][j]
+    end
+  end
+  add_contribution!(dc,d.Γ,dc_Γ)
+  dc
+end

src/Geometry/Geometry.jl

@@ -102,6 +102,7 @@ export compress_ids
 export UnstructuredGridTopology
 
 export Triangulation
+export TriangulationView
 export get_reffes
 export get_cell_coordinates
 export get_cell_ref_coordinates

test/CellDataTests/DiracDeltasTests.jl

@@ -8,6 +8,11 @@ using Gridap.ReferenceFEs
 using Gridap.Geometry
 using Gridap.CellData
 
+using Gridap.Algebra
+using Gridap.Fields
+using Gridap.FESpaces
+using Gridap.ReferenceFEs
+
 domain = (0,1,0,1)
 cells = (30,30)
 model = CartesianDiscreteModel(domain,cells)
@@ -37,6 +42,117 @@ degree = 2
 δ = DiracDelta{0}(model,tags=[2,4])
 @test sum(δ(v)) ≈ 5
 
+# Tests for DiracDelta at a generic Point in the domain #
+
+p = Point(0.2,0.3)
+δ = DiracDelta(model,p)
+@test sum(δ(v)) ≈ v(p)
+
+pvec = [p,π*p]
+δ = DiracDelta(model,pvec)
+@test sum(δ(v)) ≈ sum(v(pvec))
+
+p = Point(0.2,0.3)
+δ = DiracDelta(model,p)
+reffe = ReferenceFE(lagrangian,Float64,3)
+V = FESpace(model,reffe,conformity=:L2)
+V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
+U0 = TrialFESpace(V0)
+uh = FEFunction(U0,rand(num_free_dofs(U0)))
+vh = FEFunction(V,rand(num_free_dofs(V)))
+fe_basis = get_fe_basis(U0)
+dg_basis = get_fe_basis(V)
+
+@test sum(δ(uh)) ≈ uh(p)
+@test sum(δ(vh)) ≈ vh(p)
+@test norm(sum(δ(fe_basis)) .- fe_basis(p)) ≈ 0
+@test norm(sum(δ(dg_basis)) .- dg_basis(p)) ≈ 0
+
+pvec = [p,π*p]
+δ = DiracDelta(model,pvec)
+
+@test sum(δ(uh)) ≈ sum(map(uh,pvec))
+@test sum(δ(vh)) ≈ sum(map(vh,pvec))
+@test norm(sum(δ(fe_basis)) .- sum(map(fe_basis,pvec))) ≈ 0
+@test norm(sum(δ(dg_basis)) .- sum(map(dg_basis,pvec))) ≈ 0
+
+# comparing exact solution with that of GenericDiracDelta
+
+# 1D Test
+model = CartesianDiscreteModel((-1.,1.),(2,))
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2)
+
+# exact solution
+function u(x)
+  if x[1] < 0.0
+    return 0.0
+  else
+    return -x[1]
+  end
+end
+
+ucf = CellField(u,Ω)
+
+order = 1
+reffe = ReferenceFE(lagrangian,Float64,order)
+V = FESpace(model,reffe,conformity=:L2)
+V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
+U0 = TrialFESpace(V0,u)
+
+p = Point(0.0)
+
+δ_p = DiracDelta(model,p)
+δ_tag = DiracDelta{0}(model,tags=3)
+
+a(u,v) = ∫(∇(u)⋅∇(v))*dΩ
+l(v) = ∫(0.0*v)*dΩ + δ_p(1.0*v)
+
+op = AffineFEOperator(a,l,U0,V0)
+uh_p = solve(op)
+
+l(v) = ∫(0.0*v)*dΩ + δ_tag(1.0*v)
+op = AffineFEOperator(a,l,U0,V0)
+uh_tag = solve(op)
+
+e_p = u - uh_p
+e = uh_p - uh_tag
+err_p = ∫(e_p*e_p)*dΩ
+err = ∫(e*e)*dΩ
+@test sum(err_p) < eps()
+@test sum(err) < eps()
+
+# 2D Test - comparing with tag version and point version
+
+model = CartesianDiscreteModel((-1.,1.,-1.,1.),(3,3))
+Ω = Triangulation(model)
+dΩ = Measure(Ω,2)
+
+pvec = [Point(-1.0/3,-1.0/3), Point(-1.0/3,1.0/3), Point(1.0/3,1.0/3), Point(1.0/3,-1.0/3)]
+δ_p = DiracDelta(model,pvec)
+
+a(u,v) = ∫(∇(u)⋅∇(v))*dΩ
+l(v) = δ_p(v)
+V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
+U0 = TrialFESpace(V0,u)
+op = AffineFEOperator(a,l,U0,V0)
+uh_p = solve(op)
+
+δ_tag = DiracDelta{0}(model,tags=9)
+l(v) = δ_tag(v)
+op = AffineFEOperator(a,l,U0,V0)
+uh_tag = solve(op)
+
+err = uh_p - uh_tag
+@test sum(∫(err*err)*dΩ) < eps()
+
+# testing the faster and direct functional evalulation method
+
+f(x) = sin(norm(x))
+fcf = CellField(f,Ω)
+@test sum(δ_p(f)) ≈ sum(δ_p(fcf))
+
+
 #using Gridap
 #
 #order = 2

test/TensorValuesTests/IndexingTests.jl

@@ -14,6 +14,7 @@ v = VectorValue{4}(a)
 @test size(v) == (4,)
 @test length(v) == 4
 @test lastindex(v) == length(v)
+@test v[end] == a[end]
 
 for (k,i) in enumerate(eachindex(v))
   @test v[i] == a[k]
@@ -24,6 +25,7 @@ t = TensorValue{2}(a)
 @test size(t) == (2,2)
 @test length(t) == 4
 @test lastindex(t) == length(t)
+@test t[end] == a[end]
 
 for (k,i) in enumerate(eachindex(t))
   @test t[i] == a[k]
@@ -43,6 +45,7 @@ t = TensorValue(convert(SMatrix{2,2,Int},s))
 @test size(s) == (2,2)
 @test length(s) == 4
 @test lastindex(s) == length(s)
+@test s[end] == 22 
 
 for (k,i) in enumerate(eachindex(t))
     @test s[i] == t[k]