Merge pull request #691 from gridap/fix_issue_689

Fix issue 689

Author: fverdugo

NEWS.md

@@ -4,6 +4,11 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+### Fixed
+
+- Laplacian `Δ` operator on unstructured linear grids for quantities defined in the reference space (i.e. shape funcitons in standard FEM). Since PR [#691](https://github.com/gridap/Gridap.jl/pull/691).
+- Laplacian `Δ` operator on triangulations using `GridView` (e.g., when interpolating functions in a sub-domain or on the boundary). Since PR [#691](https://github.com/gridap/Gridap.jl/pull/691).
+
 ## [0.17.0] - 2021-10-22
 
 ### Added

src/Fields/AffineMaps.jl

@@ -53,7 +53,8 @@ function push_∇∇(∇∇a::Field,ϕ::AffineMap)
 end
 
 function push_∇∇(∇∇a::Number,Jt_inv::MultiValue{Tuple{D,D}} where D)
-  Jt_inv⋅Jt_inv⋅∇∇a
+  #Jt_inv⋅Jt_inv⋅∇∇a
+  Jt_inv⋅∇∇a⋅transpose(Jt_inv)
 end
 
 function push_∇∇(∇∇a::Number,Jt_inv::MultiValue{Tuple{D1,D2}} where {D1,D2})

src/Geometry/GridPortions.jl

@@ -180,6 +180,10 @@ function get_cell_type(grid::GridView)
   lazy_map(Reindex(get_cell_type(grid.parent)),grid.cell_to_parent_cell)
 end
 
+function get_cell_map(grid::GridView)
+  lazy_map(Reindex(get_cell_map(grid.parent)),grid.cell_to_parent_cell)
+end
+
 function get_facet_normal(grid::GridView)
   lazy_map(Reindex(get_facet_normal(grid.parent)),grid.cell_to_parent_cell)
 end

test/GridapTests/issue_689.jl

@@ -0,0 +1,40 @@
+module Issue689
+
+using Gridap
+using Gridap.ReferenceFEs
+using Gridap.Geometry
+using Gridap.CellData
+using Test
+
+function main(transf,T)
+
+  𝒯 = CartesianDiscreteModel((0,1,0,1),(10,10)) |> transf
+  labels = get_face_labeling(𝒯)
+  add_tag_from_tags!(labels,"surface",6)
+  Ω = Interior(𝒯)
+  Γ = Boundary(𝒯,tags="surface")
+
+  k = 2
+  dΩ = Measure(Ω,2*k)
+  dΓ = Measure(Γ,2*k)
+
+  reffe = ReferenceFE(lagrangian,T,k)
+  V = FESpace(Ω,reffe)
+  Q = FESpace(Γ,reffe)
+
+  u((x,y)) = x^k + y^k
+  uh = interpolate(u,V)
+  qh = interpolate(u,Q)
+  Δ_Γ(f) = x->∇∇(f)(x)[1,1]
+
+  tol = 1.0e-9
+  @test sqrt(sum( ∫( abs2( Δ(u) - Δ(uh) )  )dΩ )) < tol
+  @test sqrt(sum( ∫( abs2( Δ(u) - Δ(uh) )  )dΓ )) < tol
+  @test sqrt(sum( ∫( abs2( Δ_Γ(u) - Δ(qh) )  )dΓ )) < tol
+
+end
+
+main(simplexify,Float64)
+main(identity,Float64)
+
+end

test/GridapTests/runtests.jl

@@ -30,4 +30,6 @@ using Test
 
 @time @testset "Issue614" begin include("issue_614.jl") end
 
+@time @testset "Issue689" begin include("issue_689.jl") end
+
 end # module