Merge branch 'master' of github.com:gridap/Gridap.jl into polytopal

NEWS.md

@@ -5,7 +5,18 @@ 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).
 
-## [Unreleased]
+## [0.18.10] - 2025-01-13
+
+### Added
+
+- Added corresponding function `get_tangent_vector` to `get_normal_vector`. This method calculates the (unique up to sign) tangential unit vector to edges in 2D meshes, by rotating the normal (nx, ny) -> (ny, -nx). Since PR[#1071](https://github.com/gridap/Gridap.jl/pull/1071).
+
+## [0.18.9] - 2025-01-13
+
+
+### Fixed
+
+- BUG in `FineToCoarseFields.jl`. Since PR[#1074](https://github.com/gridap/Gridap.jl/pull/1074)
 
 ### Added
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
-version = "0.18.8"
+version = "0.18.9"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/Adaptivity/FineToCoarseFields.jl

@@ -108,7 +108,7 @@ function Arrays.evaluate!(cache,a::FineToCoarseField,x::AbstractArray{<:Point})
     field_id = id_map[child_id]
     if field_id != 0
       fi = getindex!(fi_cache,fields,field_id)
-      mi = getindex!(mi_cache,cmaps,field_id)
+      mi = getindex!(mi_cache,cmaps,child_id)
       zi = Fields.evaluate!(zi_cache,mi,xi)
       y[i] = Fields.evaluate!(yi_cache,fi,zi)
     end

src/CellData/CellData.jl

@@ -61,6 +61,7 @@ export Integrand
 export ∫
 export CellDof
 export get_normal_vector
+export get_tangent_vector
 export get_cell_measure
 export Interpolable
 export KDTreeSearch

src/CellData/CellFields.jl

@@ -113,11 +113,20 @@ end
 
 function get_normal_vector(trian::Triangulation)
   cell_normal = get_facet_normal(trian)
-  get_normal_vector(trian,cell_normal)
+  get_normal_vector(trian, cell_normal)
 end
 
-function get_normal_vector(trian::Triangulation,cell_normal::AbstractArray)
-  GenericCellField(cell_normal,trian,ReferenceDomain())
+function get_tangent_vector(trian::Triangulation)
+  cell_tangent = get_edge_tangent(trian)
+  get_tangent_vector(trian, cell_tangent)
+end
+
+function get_normal_vector(trian::Triangulation,cell_vectors::AbstractArray)
+  GenericCellField(cell_vectors,trian,ReferenceDomain())
+end
+
+function get_tangent_vector(trian::Triangulation,cell_vectors::AbstractArray)
+  GenericCellField(cell_vectors,trian,ReferenceDomain())
 end
 
 evaluate!(cache,f::Function,x::CellPoint) = CellField(f,get_triangulation(x))(x)
@@ -712,6 +721,12 @@ function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair)
   SkeletonPair(plus,minus)
 end
 
+function get_tangent_vector(trian::Triangulation,cell_tangent::SkeletonPair)
+  plus = get_normal_vector(trian,cell_tangent.plus)
+  minus = get_normal_vector(trian,cell_tangent.minus)
+  SkeletonPair(plus,minus)
+end
+
 for op in (:outer,:*,:dot)
   @eval begin
     ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)

src/Exports.jl

@@ -141,6 +141,7 @@ using Gridap.TensorValues: ⊗; export ⊗
 @publish CellData mean
 @publish CellData update_state!
 @publish CellData get_normal_vector
+@publish CellData get_tangent_vector
 using Gridap.CellData: ∫; export ∫
 @publish CellData get_cell_measure
 @publish CellData get_physical_coordinate

src/Geometry/BoundaryTriangulations.jl

@@ -237,6 +237,24 @@ function get_facet_normal(trian::BoundaryTriangulation, face_glue::FaceToCellGlu
   Fields.MemoArray(face_s_n)
 end
 
+function get_edge_tangent(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue)
+  face_s_n = get_facet_normal(trian, boundary_trian_glue)
+  rotate_normal_2d(n) = VectorValue(n[2], -n[1])
+  face_s_t = lazy_map(Operation(rotate_normal_2d), face_s_n)
+  return Fields.MemoArray(face_s_t)
+end
+
+function get_edge_tangent(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A}
+  Dp != 2 && error("get_edge_tangent only implemented for 2D")
+  glue = trian.glue
+  get_edge_tangent(trian, glue)
+end
+
+function get_facet_normal(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A}
+  glue = trian.glue
+  get_facet_normal(trian, glue)
+end
+
 function push_normal(invJt,n)
   v = invJt⋅n
   m = sqrt(inner(v,v))

src/Geometry/Geometry.jl

@@ -73,6 +73,7 @@ import Gridap.ReferenceFEs: num_cell_dims
 import Gridap.ReferenceFEs: num_point_dims
 import Gridap.ReferenceFEs: simplexify
 import Gridap.ReferenceFEs: get_facet_normal
+import Gridap.ReferenceFEs: get_edge_tangent
 import Gridap.ReferenceFEs: Quadrature
 
 export GridTopology
@@ -110,6 +111,7 @@ export get_cell_ref_coordinates
 export get_cell_reffe
 export get_cell_shapefuns
 export get_facet_normal
+export get_edge_tangent
 export test_triangulation
 export get_cell_map
 

src/Geometry/SkeletonTriangulations.jl

@@ -94,6 +94,12 @@ function get_facet_normal(trian::SkeletonTriangulation)
   SkeletonPair(plus,minus)
 end
 
+function get_edge_tangent(trian::SkeletonTriangulation)
+  plus = get_edge_tangent(trian.plus)
+  minus = get_edge_tangent(trian.minus)
+  SkeletonPair(plus,minus)
+end
+
 # Related with CompositeTriangulation
 function _compose_glues(rglue::FaceToFaceGlue,dglue::SkeletonPair)
   plus = _compose_glues(rglue,dglue.plus)

test/AdaptivityTests/FineToCoarseFieldsTests.jl

@@ -100,4 +100,91 @@ xHh = change_domain(xH,get_triangulation(modelh),ReferenceDomain())
 evaluate(Gridap.CellData.get_data(xHh)[1],[Point(0.0,0.0),Point(0.5,0.5)])
 evaluate(Gridap.CellData.get_data(xHh)[1],Point(0.5,0.5))
 
+function test_PR_1074()
+
+    function setup_coarse_discrete_model()
+        ptr  = [ 1, 5 ]
+        data = [ 1,2,3,4  ]
+        cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+        node_coordinates = Vector{Point{2,Float64}}(undef,4)
+        node_coordinates[1]=Point{2,Float64}(0.0,0.0)
+        node_coordinates[2]=Point{2,Float64}(1.0,0.0)
+        node_coordinates[3]=Point{2,Float64}(0.0,1.0)
+        node_coordinates[4]=Point{2,Float64}(1.0,1.0)
+
+        polytope=QUAD
+        scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
+        cell_types=collect(Fill(1,length(cell_vertex_lids)))
+        cell_reffes=[scalar_reffe]
+        grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
+                                                cell_vertex_lids,
+                                                cell_reffes,
+                                                cell_types,
+                                                Gridap.Geometry.NonOriented())
+        Gridap.Geometry.UnstructuredDiscreteModel(grid)
+    end
+
+    function setup_adapted_discrete_model(parent)
+        ptr  = [ 1, 5, 9 ]
+        data = [ 1,2,3,4, 2,5,4,6 ] 
+        cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+        node_coordinates = Vector{Point{2,Float64}}(undef,6)
+        node_coordinates[1]=Point{2,Float64}(0.0,0.5)
+        node_coordinates[2]=Point{2,Float64}(0.5,0.5)
+        node_coordinates[3]=Point{2,Float64}(0.0,1.0)
+        node_coordinates[4]=Point{2,Float64}(0.5,1.0)
+        node_coordinates[5]=Point{2,Float64}(1.0,0.5)
+        node_coordinates[6]=Point{2,Float64}(1.0,1.0)
+
+        polytope=QUAD
+        scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
+        cell_types=collect(Fill(1,length(cell_vertex_lids)))
+        cell_reffes=[scalar_reffe]
+        grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
+                                                cell_vertex_lids,
+                                                cell_reffes,
+                                                cell_types,
+                                                Gridap.Geometry.NonOriented())
+        model=Gridap.Geometry.UnstructuredDiscreteModel(grid)
+
+        n2o_faces_map=Vector{Vector{Int64}}(undef,3)
+        n2o_faces_map[3]=[1,1]
+        n2o_cell_to_child_id=[2,3]
+
+        ref_rules = [Gridap.Adaptivity.RefinementRule(Gridap.ReferenceFEs.LagrangianRefFE(Float64,QUAD,1),2)]
+
+        glue=Gridap.Adaptivity.AdaptivityGlue(Gridap.Adaptivity.RefinementGlue(),
+                                        n2o_faces_map,
+                                        n2o_cell_to_child_id,
+                                        ref_rules)
+        Gridap.Adaptivity.AdaptedDiscreteModel(model,parent,glue)
+    end
+
+    coarse_model=setup_coarse_discrete_model()
+    fine_model=setup_adapted_discrete_model(coarse_model)
+
+    order=0
+
+    VH = FESpace(coarse_model,
+                ReferenceFE(raviart_thomas,Float64,order),
+                conformity=:Hdiv)
+
+    UH = TrialFESpace(VH)
+
+    Vh = FESpace(fine_model,
+                ReferenceFE(raviart_thomas,Float64,order),
+                conformity=:Hdiv)
+
+    Uh = TrialFESpace(Vh)
+
+    u(x) = VectorValue(x[1],x[2])
+    uh  = interpolate(u,Uh)
+    uH  = interpolate(u,UH)
+    uhH = interpolate(uh,UH) 
+
+    @test get_free_dof_values(uhH)[2] ≈ 1.0
+end
+
+test_PR_1074()
+
 end

test/GeometryTests/BoundaryTriangulationsTests.jl

@@ -110,6 +110,11 @@ face_to_nvec_s = lazy_map(evaluate,face_to_nvec,face_to_s)
 test_array(face_to_nvec_s,collect(face_to_nvec_s))
 @test isa(face_to_nvec_s,Geometry.FaceCompressedVector)
 
+face_to_tvec = get_edge_tangent(btrian)
+face_to_tvec_s = lazy_map(evaluate,face_to_tvec,face_to_s)
+test_array(face_to_tvec_s,collect(face_to_tvec_s))
+@test isa(face_to_tvec_s,Geometry.FaceCompressedVector)
+
 #print_op_tree(face_shapefuns_q)
 #print_op_tree(face_grad_shapefuns_q)
 #print_op_tree(face_to_nvec_s)
@@ -147,6 +152,16 @@ r = Vector{Point{2,Float64}}[
   [(0.0,1.0),(0.0,1.0)],[(1.0,-0.0),(1.0,-0.0)]]
 test_array(nvec_s,r)
 
+tvec = get_edge_tangent(btrian)
+
+tvec_s = lazy_map(evaluate,tvec,s)
+
+rotate_90(x) = [Point(x[1][2], -x[1][1]), Point(x[2][2], -x[2][1])]
+
+rr = rotate_90.(r)
+
+test_array(rr, tvec_s)
+
 cellids = collect(1:num_cells(model))
 
 face_to_cellid = lazy_map(Reindex(cellids),glue.tface_to_mface)

test/GeometryTests/SkeletonTriangulationsTests.jl

@@ -1,6 +1,7 @@
 module SkeletonTriangulationsTests
 
 using Test
+using Gridap
 using Gridap.TensorValues
 using Gridap.Arrays
 using Gridap.Fields
@@ -52,6 +53,7 @@ vn = get_facet_normal(vtrian)
 Ω = Triangulation(model)
 Γ = BoundaryTriangulation(model)
 Λ = SkeletonTriangulation(Γ)
+normal = get_normal_vector(Λ)
 @test Λ.rtrian === Γ
 @test isa(Λ.dtrian,SkeletonTriangulation)
 glue = get_glue(Λ,Val(3))
@@ -110,13 +112,27 @@ itrian = InterfaceTriangulation(Ω_in,Ω_out)
 #writevtk(trian,"trian",celldata=["inout"=>cell_to_inout])
 #writevtk(itrian,"itrian",nsubcells=10,cellfields=["ni"=>ni,"nl"=>nl,"nr"=>nr])
 
+ti = get_tangent_vector(itrian)
+tl = get_tangent_vector(ltrian)
+tr = get_tangent_vector(rtrian)
+
+@test ti isa SkeletonPair
+@test tl isa Gridap.CellData.GenericCellField
+@test tr isa Gridap.CellData.GenericCellField 
+
 reffe = LagrangianRefFE(Float64,QUAD,(2,2))
 conf = CDConformity((CONT,DISC))
 face_own_dofs = get_face_own_dofs(reffe,conf)
 strian = SkeletonTriangulation(model,reffe,face_own_dofs)
 test_triangulation(strian)
 ns = get_facet_normal(strian)
+ts = get_edge_tangent(strian)
 @test length(ns.⁺) == num_cells(strian)
+@test length(ts.⁺) == num_cells(strian)
+# Test orthogonality
+should_be_zero = lazy_map(Broadcasting(Operation(dot)), ns.plus, ts.plus)
+ndot_t_evaluated = evaluate(should_be_zero, VectorValue.(0:0.05:1))
+@test maximum(ndot_t_evaluated) < 1e-15
 
 #using Gridap.Visualization
 #writevtk(strian,"strian",cellfields=["normal"=>ns])