Adding support to DiracDelta for generic points in the domain

NEWS.md

@@ -10,7 +10,7 @@ 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)
-
+- Added support to construct `DiracDelta` at generic `Point`(s) in the domain. Since PR [#837](https://github.com/gridap/Gridap.jl/pull/837)
 ## [0.17.14] - 2022-07-29
 
 ### Added

src/CellData/DiracDeltas.jl

@@ -1,5 +1,8 @@
+abstract type DiracDeltaSupportType <: GridapType end
+struct IsGridEntity <: DiracDeltaSupportType end
+struct NotGridEntity <: DiracDeltaSupportType end
 
-struct DiracDelta{D} <: GridapType
+struct DiracDelta{D,S<:DiracDeltaSupportType} <: GridapType
   Γ::Triangulation{D}
   dΓ::Measure
 end
@@ -24,7 +27,7 @@ function DiracDelta{D}(
   trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)
   Γ = BoundaryTriangulation(trian,glue)
   dΓ = Measure(Γ,degree)
-  DiracDelta(Γ,dΓ)
+  DiracDelta{D,IsGridEntity}(Γ,dΓ)
 end
 
 function DiracDelta{D}(
@@ -76,3 +79,33 @@ end
 function get_triangulation(d::DiracDelta)
   d.Γ
 end
+
+# For handling DiracDelta at a generic Point in the domain #
+
+function DiracDelta(x::Point{D,T}, model::DiscreteModel{D}) where {D,T}
+  trian = Triangulation(model)
+  cache1 = _point_to_cell_cache(KDTreeSearch(),trian)
+  cell = _point_to_cell!(cache1, x)
+  point_grid = UnstructuredGrid([x])
+  point_model = UnstructuredDiscreteModel(point_grid)
+  point_trian = Triangulation(point_model)
+  dx = Measure(point_trian,1)
+  DiracDelta{0,NotGridEntity}(point_trian,dx)
+end
+
+function DiracDelta(v::Vector{Point{D,T}},model::DiscreteModel{D}) where {D,T}
+  trian = Triangulation(model)
+  cache1 = _point_to_cell_cache(KDTreeSearch(),trian)
+  cell = map(x -> _point_to_cell!(cache1, x), v)
+  point_grid = UnstructuredGrid(v)
+  point_model = UnstructuredDiscreteModel(point_grid)
+  point_trian = Triangulation(point_model)
+  dx = Measure(point_trian,1)
+  DiracDelta{0,NotGridEntity}(point_trian,dx)
+end
+
+function evaluate!(cache,d::DiracDelta{0,NotGridEntity},f::CellField)
+  d_f = lazy_map(f,d.Γ.grid.node_coordinates)
+  dc = DomainContribution()
+  add_contribution!(dc, d.Γ, d_f)
+end

test/CellDataTests/DiracDeltasTests.jl

@@ -37,6 +37,22 @@ degree = 2
 δ = DiracDelta{0}(model,tags=[2,4])
 @test sum(δ(v)) ≈ 5
 
+# Tests for DiracDelta at a generic Point in the domain #
+
+p = VectorValue(0.2,0.3)
+δ = DiracDelta(p,model)
+@test sum(δ(v)) == v(p)
+
+pvec = [p,π*p]
+δ = DiracDelta(pvec,model)
+@test sum(δ(v)) == sum(v(pvec))
+
+# below passes but cannot use FESpace and FEFunction here
+# reffe = ReferenceFE(lagrangian,Float64,3)
+# V = FESpace(model,reffe,conformity=:L2)
+# uh = FEFunction(V,rand(num_free_dofs(V)))
+# @test sum(d(uh)) == uh(p)
+
 #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]