Merge branch 'gridap:master' into q_tensor

NEWS.md

@@ -5,13 +5,20 @@ 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]
+
+### Added
+
+- Added MacroFElements. These are defined as having the basis/dof-basis of a FESpace created on top of a RefinementRule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
+- Added Barycentric refinement rule in 2D and 3D. Added Simplexify refinement rule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
+
 ## [0.18.6] - 2024-08-29
 
 ### Fixed
 
 - Improved performance of PR[#967](https://github.com/gridap/Gridap.jl/pull/967). Along the way, opened the door to Triangulations of different type in SkeletonTriangulation. Since PR[#1026](https://github.com/gridap/Gridap.jl/pull/1026).
 
-## [0.18.5] - 2024-08-28 
+## [0.18.5] - 2024-08-28
 
 ### Changed
 

Project.toml

@@ -42,7 +42,7 @@ FastGaussQuadrature = "0.4.2, 1"
 FileIO = "1.2.2, 1.3, 1.4"
 FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 0.13, 1"
 ForwardDiff = "0.10.10"
-JLD2 = "0.1.11, 0.3, 0.4"
+JLD2 = "0.1.11, 0.3, 0.4, 0.5"
 JSON = "0.21.0"
 LineSearches = "7.0.1"
 NLsolve = "4.3.0"

README.md

@@ -63,7 +63,7 @@ pkg> add Gridap
 
 - [GridapDistributed](https://github.com/gridap/GridapDistributed.jl) Distributed-memory extension of Gridap.
 - [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) Embedded finite elements in Julia.
-- [GridapGmsh](https://github.com/gridap/GridapGmsh.jl) Generate a FE mesh with [GMSH](www.gmsh.info) and use it in Gridap.
+- [GridapGmsh](https://github.com/gridap/GridapGmsh.jl) Generate a FE mesh with [GMSH](https://www.gmsh.info) and use it in Gridap.
 - [GridapMakie](https://github.com/gridap/GridapMakie.jl) Makie plotting recipes for Gridap.
 - [GridapPardiso](https://github.com/gridap/GridapPardiso.jl) Use the [Intel Pardiso MKL direct sparse solver](https://software.intel.com/en-us/mkl-developer-reference-fortran-intel-mkl-pardiso-parallel-direct-sparse-solver-interface) in Gridap.
 - [GridapPETSc](https://github.com/gridap/GridapPETSc.jl) Use [PETSc](https://petsc.org/) linear and nonlinear solvers in Gridap.

src/Adaptivity/CompositeQuadratures.jl

@@ -71,8 +71,6 @@ the quadrature `quad` into the subcells of the RefinementRule.
 function CompositeQuadrature(
   quad::Quadrature,rr::RefinementRule
 )
-  @check ReferenceFEs.get_polytope(quad) === ReferenceFEs.get_polytope(rr)
-
   weights = get_weights(quad)
   cpoints = get_coordinates(quad)
 

src/Adaptivity/MacroFEs.jl

@@ -348,6 +348,7 @@ function MacroReferenceFE(
 
   ndofs = num_free_dofs(space)
   poly = get_polytope(rrule)
+  # This is a hack to be able to compute the orders
   prebasis = FineToCoarseArray(rrule,collect(map(get_prebasis,reffes)))
   metadata = (rrule,conn,face_own_dofs,face_own_perms)
 

test/AdaptivityTests/FineToCoarseFieldsTests.jl

@@ -91,13 +91,13 @@ eh2 = u_c - u_fc
 eh3 = u_c - u_fc2 
 @test sum(∫(eh3⋅eh3)*dΩ_c) < 1.e-12
 
-modelH=CartesianDiscreteModel((0,1,0,1),(1,1))
-modelh=refine(modelH,2)
-reffe=LagrangianRefFE(Float64,QUAD,1)
-XH  = TestFESpace(modelH,reffe)
-xH  = get_fe_basis(XH)
+modelH = CartesianDiscreteModel((0,1,0,1),(1,1))
+modelh = refine(modelH,2)
+reffe = LagrangianRefFE(Float64,QUAD,1)
+XH = TestFESpace(modelH,reffe)
+xH = get_fe_basis(XH)
 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.0,0.0),Point(0.5,0.5)])
+evaluate(Gridap.CellData.get_data(xHh)[1],Point(0.5,0.5))
 
 end

test/AdaptivityTests/MacroFEStokesTests.jl

@@ -11,15 +11,14 @@ function main(Dc,reftype)
   @assert reftype ∈ [:barycentric,:powellsabin]
 
   u_sol(x) = (Dc == 2) ? VectorValue(x[1],-x[2]) : VectorValue(x[1],-x[2],0.0)
-  p_sol(x) = (x[1] - 1.0/2.0)
+  p_sol(x) = x[1] - 1.0/2.0
 
   domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1)
   nc = (Dc == 2) ? (2,2) : (1,1,1)
   model = simplexify(CartesianDiscreteModel(domain,nc))
 
-  min_order = (reftype == :barycentric) ? Dc : Dc-1
-  order = max(2,min_order)
   poly  = (Dc == 2) ? TRI : TET
+  order = (reftype == :barycentric) ? Dc : Dc-1
   rrule = (reftype == :barycentric) ? Adaptivity.BarycentricRefinementRule(poly) : Adaptivity.PowellSabinRefinementRule(poly)
 
   subreffes_u = Fill(LagrangianRefFE(VectorValue{Dc,Float64},poly,order),Adaptivity.num_subcells(rrule))
@@ -28,7 +27,7 @@ function main(Dc,reftype)
   subreffes_p = Fill(LagrangianRefFE(Float64,poly,order-1),Adaptivity.num_subcells(rrule))
   reffe_p = Adaptivity.MacroReferenceFE(rrule,subreffes_p;conformity=L2Conformity())
 
-  qdegree = 2*order
+  qdegree = 2*(order-1)
   quad  = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,qdegree)
 
   V = FESpace(model,reffe_u,dirichlet_tags=["boundary"])
@@ -54,15 +53,13 @@ function main(Dc,reftype)
   uh, ph = xh
   eh_u = uh - u_sol
   eh_p = ph - p_sol
-  println(sum(∫(eh_u⋅eh_u)dΩ))
-  println(sum(∫(eh_p⋅eh_p)dΩ))
   @test sum(∫(eh_u⋅eh_u)dΩ) < 1.e-10
   if reftype != :powellsabin
     @test sum(∫(eh_p*eh_p)dΩ) < 1.e-10
   end
 end
 
-# NOTE: Powell-Sabin split not working yet. The issue is we woudl need a global cellmap 
+# NOTE: Powell-Sabin split not working yet. The issue is we would need a global cellmap 
 # directly from the sub-cells to the physical domain (due to how the split is built). 
 # This is something I may do in the future. 
 

test/AdaptivityTests/MacroFETests.jl

@@ -13,7 +13,11 @@ function test_macro_reffe(model,fmodel,rrule,order)
   reffe = LagrangianRefFE(Float64,poly,order)
   sub_reffes = Fill(reffe,num_subcells(rrule))
   macro_reffe = Adaptivity.MacroReferenceFE(rrule,sub_reffes)
+
   macro_quad  = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,2*order)
+  macro_quad_bis = Adaptivity.CompositeQuadrature(Quadrature(poly,4*order),rrule)
+  ReferenceFEs.test_quadrature(macro_quad)
+  ReferenceFEs.test_quadrature(macro_quad_bis)
 
   Ω = Triangulation(model)
   Ωf = Triangulation(fmodel)

test/AdaptivityTests/RefinementRulesTests.jl

@@ -24,4 +24,13 @@ for poly in polys
   end
 end
 
+rr_bc2 = Adaptivity.BarycentricRefinementRule(TRI)
+Adaptivity.test_refinement_rule(rr_bc2)
+
+rr_bc3 = Adaptivity.BarycentricRefinementRule(TET)
+Adaptivity.test_refinement_rule(rr_bc3)
+
+rr_ps2 = Adaptivity.PowellSabinRefinementRule(TRI)
+rr_ps3 = Adaptivity.PowellSabinRefinementRule(TET)
+
 end