Merge pull request #485 from gridap/diracs_deltas

Added DiracDelta

src/CellData/CellData.jl

@@ -74,6 +74,8 @@ export get_physical_coordinate
 export CellState
 export update_state!
 
+export DiracDelta
+
 include("CellDataInterface.jl")
 
 include("CellFields.jl")
@@ -84,6 +86,8 @@ include("CellStates.jl")
 
 include("DomainContributions.jl")
 
+include("DiracDeltas.jl")
+
 include("CellDofs.jl")
 
 include("AttachDirichlet.jl")

src/CellData/DiracDeltas.jl

@@ -0,0 +1,77 @@
+
+struct DiracDelta{D} <: GridapType
+  Γ::Triangulation{D}
+  dΓ::LebesgueMeasure
+end
+
+function DiracDelta{D}(
+  model::DiscreteModel,
+  face_to_bgface::AbstractVector{<:Integer},
+  degree::Integer) where D
+
+  @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))).
+  """
+  topo = get_grid_topology(model)
+  bgface_grid = Grid(ReferenceFE{D},model)
+  face_trian = RestrictedTriangulation(bgface_grid,face_to_bgface)
+  cell_trian = Triangulation(model)
+  bgface_to_lcell = Fill(1,num_faces(model,D))
+  glue = Geometry.FaceToCellGlue(topo,cell_trian,face_trian,face_to_bgface,bgface_to_lcell)
+  Γ = BoundaryTriangulation(face_trian,cell_trian,glue)
+  dΓ = LebesgueMeasure(Γ,degree)
+  DiracDelta(Γ,dΓ)
+end
+
+function DiracDelta{D}(
+  model::DiscreteModel,
+  bgface_to_mask::AbstractVector{Bool},
+  degree::Integer) where D
+
+  face_to_bgface = findall(bgface_to_mask)
+  DiracDelta{D}(model,face_to_bgface,degree)
+end
+
+function DiracDelta{D}(
+  model::DiscreteModel,
+  labeling::FaceLabeling,
+  degree::Integer;tags) where D
+
+  @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))).
+  """
+  face_to_mask = get_face_mask(labeling,tags,D)
+  DiracDelta{D}(model,face_to_mask,degree)
+end 
+
+function DiracDelta{D}(model::DiscreteModel,degree::Integer;tags) where D
+  labeling = get_face_labeling(model)
+  DiracDelta{D}(model,labeling,degree,tags=tags)
+end
+
+function DiracDelta{0}(model::DiscreteModel,labeling::FaceLabeling;tags)
+  degree = 0
+  DiracDelta{0}(model,labeling,degree;tags=tags)
+end
+
+function DiracDelta{0}(model::DiscreteModel;tags)
+  degree = 0
+  DiracDelta{0}(model,degree;tags=tags)
+end
+
+function (d::DiracDelta)(f::CellField)
+  evaluate(d,f)
+end
+
+function evaluate!(cache,d::DiracDelta,f::CellField)
+  ∫(f)*d.dΓ
+end
+
+function get_triangulation(d::DiracDelta)
+  d.Γ
+end
+

src/Exports.jl

@@ -130,6 +130,7 @@ using Gridap.TensorValues: ⊗; export ⊗
 using Gridap.CellData: ∫; export ∫
 @publish CellData get_cell_measure
 @publish CellData get_physical_coordinate
+@publish CellData DiracDelta
 
 @publish FESpaces FESpace
 @publish FESpaces TrialFESpace

src/Geometry/BoundaryTriangulations.jl

@@ -21,9 +21,10 @@ function FaceToCellGlue(
   bgface_to_lcell::AbstractVector)
 
   D = num_cell_dims(cell_trian)
-  bgface_to_cells = get_faces(topo,D-1,D)
-  cell_to_bgfaces = get_faces(topo,D,D-1)
-  cell_to_lface_to_pindex = Table(get_cell_permutations(topo,D-1))
+  cD = num_cell_dims(face_trian)
+  bgface_to_cells = get_faces(topo,cD,D)
+  cell_to_bgfaces = get_faces(topo,D,cD)
+  cell_to_lface_to_pindex = Table(get_cell_permutations(topo,cD))
 
   bgface_to_cell = lazy_map(getindex,bgface_to_cells, bgface_to_lcell)
   bgface_to_lface = find_local_index(bgface_to_cell, cell_to_bgfaces)
@@ -32,7 +33,7 @@ function FaceToCellGlue(
   face_to_lface = collect(Int8,lazy_map(Reindex(bgface_to_lface), face_to_bgface))
   face_to_lcell = collect(Int8,lazy_map(Reindex(bgface_to_lcell), face_to_bgface))
 
-  f = (p)->fill(Int8(UNSET),num_faces(p,D-1))
+  f = (p)->fill(Int8(UNSET),num_faces(p,cD))
   ctype_to_lface_to_ftype = map( f, get_reffes(cell_trian) )
   face_to_ftype = get_cell_type(face_trian)
   cell_to_ctype = get_cell_type(cell_trian)
@@ -96,7 +97,7 @@ struct BoundaryTriangulation{Dc,Dp,Gf,Gc,G} <: Triangulation{Dc,Dp}
 
     @assert TriangulationStyle(cell_trian) == BackgroundTriangulation()
     @assert num_point_dims(face_trian) == num_point_dims(cell_trian)
-    @assert num_cell_dims(face_trian) == num_cell_dims(cell_trian) - 1
+    #@assert num_cell_dims(face_trian) == num_cell_dims(cell_trian) - 1
 
     Dc = num_cell_dims(face_trian)
     Dp = num_point_dims(face_trian)

src/TensorValues/Operations.jl

@@ -374,9 +374,9 @@ end
     Meta.parse("TensorValue{$D,$Z}($str)")
 end
 
-function outer(a::VectorValue{0,Ta},b::VectorValue{1,Tb}) where {Ta,Tb}
+function outer(a::VectorValue{0,Ta},b::VectorValue{D,Tb}) where {Ta,Tb,D}
   T = promote_type(Ta,Tb)
-  TensorValue{0,1,T}()
+  TensorValue{0,D,T}()
 end
 
 function outer(a::VectorValue{0,Ta},b::Tb) where {Ta,Tb<:Number}
@@ -483,7 +483,7 @@ end
 """
 meas(a::MultiValue{Tuple{D}}) where D = sqrt(inner(a,a))
 meas(a::MultiValue{Tuple{D,D}}) where D = abs(det(a))
-meas(a::TensorValue{0,1,T}) where T = one(T)
+meas(a::TensorValue{0,D,T}) where {T,D} = one(T)
 
 function meas(v::MultiValue{Tuple{1,2}})
   n1 = v[1,2]

test/CellDataTests/DiracDeltasTests.jl

@@ -0,0 +1,61 @@
+module DiracDeltasTests
+
+using Test
+using Gridap.Helpers
+using Gridap.Arrays
+using Gridap.ReferenceFEs
+using Gridap.Geometry
+using Gridap.CellData
+
+domain = (0,1,0,1)
+cells = (30,30)
+model = CartesianDiscreteModel(domain,cells)
+
+Ω = Triangulation(model)
+
+vfun(x) = x[1] + 3*x[2]
+v = CellField(vfun,Ω)
+
+degree = 2
+δ = DiracDelta{1}(model,degree,tags=5)
+@test sum(δ(v)) ≈ 0.5
+
+#using Gridap.Visualization
+#writevtk(get_triangulation(δ),"trian_d")
+#writevtk(Ω,"trian")
+
+δ = DiracDelta{0}(model,tags=2)
+@test sum(δ(v)) ≈ 1
+
+δ = DiracDelta{0}(model,tags=4)
+@test sum(δ(v)) ≈ 4
+
+δ = DiracDelta{0}(model,tags=[2,4])
+@test sum(δ(v)) ≈ 5
+
+#using Gridap
+#
+#order = 2
+#reffe = ReferenceFE(:Lagrangian,VectorValue{2,Float64},order)
+#V = FESpace(model,reffe,dirichlet_tags=[1,2,5])
+#
+#dΩ = LebesgueMeasure(Ω,2*order)
+#
+#δ = DiracDelta{0}(model,tags=4)
+#
+#const E = 1
+#const ν = 0.33
+#const λ = (E*ν)/((1+ν)*(1-2*ν))
+#const μ = E/(2*(1+ν))
+#σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
+#
+#a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
+#l(v) = δ( VectorValue(-1,0)⋅v )
+#
+#op = AffineFEOperator(a,l,V,V)
+#uh = solve(op)
+#
+#writevtk(Ω,"results",cellfields=["uh"=>uh])
+#writevtk(δ.Γ,"d",cellfields=["uh"=>uh])
+
+end # module

test/CellDataTests/runtests.jl

@@ -8,6 +8,8 @@ using Test
 
 @testset "DomainContributions" begin include("DomainContributionsTests.jl") end
 
+@testset "DiracDeltas" begin include("DiracDeltasTests.jl") end
+
 @testset "CellStates" begin include("CellStatesTests.jl") end
 
 @testset "CellDofsTests" begin include("CellDofsTests.jl") end