From 5a85738cf51dd62d2cc2240f363821e247ca4179 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Fri, 27 Nov 2020 18:03:57 +0100 Subject: [PATCH] Added DiracDelta Fixes issue #408 --- src/CellData/CellData.jl | 4 ++ src/CellData/DiracDeltas.jl | 77 ++++++++++++++++++++++++++ src/Exports.jl | 1 + src/Geometry/BoundaryTriangulations.jl | 11 ++-- src/TensorValues/Operations.jl | 6 +- test/CellDataTests/DiracDeltasTests.jl | 61 ++++++++++++++++++++ test/CellDataTests/runtests.jl | 2 + 7 files changed, 154 insertions(+), 8 deletions(-) create mode 100644 src/CellData/DiracDeltas.jl create mode 100644 test/CellDataTests/DiracDeltasTests.jl diff --git a/src/CellData/CellData.jl b/src/CellData/CellData.jl index da2920f18..541634b96 100644 --- a/src/CellData/CellData.jl +++ b/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") diff --git a/src/CellData/DiracDeltas.jl b/src/CellData/DiracDeltas.jl new file mode 100644 index 000000000..5cda2974a --- /dev/null +++ b/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 + diff --git a/src/Exports.jl b/src/Exports.jl index 726ddedae..ba366e291 100644 --- a/src/Exports.jl +++ b/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 diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl index ae46efd52..d8035db73 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/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) diff --git a/src/TensorValues/Operations.jl b/src/TensorValues/Operations.jl index ce12ef77a..8418f38b3 100644 --- a/src/TensorValues/Operations.jl +++ b/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] diff --git a/test/CellDataTests/DiracDeltasTests.jl b/test/CellDataTests/DiracDeltasTests.jl new file mode 100644 index 000000000..516e81a88 --- /dev/null +++ b/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 diff --git a/test/CellDataTests/runtests.jl b/test/CellDataTests/runtests.jl index dfe5f490a..8e3245605 100644 --- a/test/CellDataTests/runtests.jl +++ b/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