More first-order CR type elements

docs/src/assets/references.bib

@@ -123,3 +123,23 @@ @article{Mu:2014:IP
 keywords = {Discontinuous Galerkin, Finite element, Interior penalty, Second-order elliptic equations, Hybrid mesh},
 abstract = {This paper provides a theoretical foundation for interior penalty discontinuous Galerkin methods for second-order elliptic equations on very general polygonal or polyhedral meshes. The mesh can be composed of any polygons or polyhedra that satisfy certain shape regularity conditions characterized in a recent paper by two of the authors, Wang and Ye (2012) [11]. The usual H1-conforming finite element methods on such meshes are either very complicated or impossible to implement in practical computation. The interior penalty discontinuous Galerkin method provides a simple and effective alternative approach which is efficient and robust. Results with such general meshes have important application in computational sciences.}
 }
+@article{RanTur:1992:snq,
+  title={Simple nonconforming quadrilateral Stokes element},
+  author={Rannacher, Rolf and Turek, Stefan},
+  journal={Numerical Methods for Partial Differential Equations},
+  volume={8},
+  number={2},
+  pages={97--111},
+  year={1992},
+  publisher={Wiley Online Library}
+}
+@article{CroRav:1973:cnf,
+  title={Conforming and nonconforming finite element methods for solving the stationary Stokes equations I},
+  author={Crouzeix, Michel and Raviart, P-A},
+  journal={Revue fran{\c{c}}aise d'automatique informatique recherche op{\'e}rationnelle. Math{\'e}matique},
+  volume={7},
+  number={R3},
+  pages={33--75},
+  year={1973},
+  publisher={EDP Sciences}
+}

docs/src/reference/interpolations.md

@@ -12,3 +12,14 @@ getrefdim(::Interpolation)
 getrefshape
 getorder
 ```
+
+Implemented interpolations:
+
+```@docs
+Lagrange
+Serendipity
+DiscontinuousLagrange
+BubbleEnrichedLagrange
+CrouzeixRaviart
+RannacherTurek
+```

src/Quadrature/gaussquad_prism_table.jl

@@ -59,7 +59,7 @@ function _get_gauss_prismdata_polyquad(n::Int)
              -0.79828210803458293816870337538129563058   -0.79828210803458293816870337538129563058    0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
         ]
     elseif n == 6
-        wx = [
+        xw = [
             -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.52662270480497475869798007748379244322    0.22269313409222502250038629788629976598
             -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.52662270480497475869798007748379244322    0.22269313409222502250038629788629976598
             -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.99083630081924474869286718300205167763    0.13839610937412451172575590081432159364
@@ -90,7 +90,7 @@ function _get_gauss_prismdata_polyquad(n::Int)
             -0.59485220654953844609181984939099585377   -0.93064127511790609994722013739838094335    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
         ]
     elseif n == 7
-        wx = [
+        xw = [
             -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.98022806959089160171504882914128008603    0.11378272445075715563392222208600913026
             -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.98022806959089160171504882914128008603    0.11378272445075715563392222208600913026
             -0.98332480906795705560590831565479581293    0.96664961813591411121181663130959162585                                           0   0.024472968692380158977782543188506640362

src/exports.jl

@@ -13,6 +13,7 @@ export
     RefPyramid,
     BubbleEnrichedLagrange,
     CrouzeixRaviart,
+    RannacherTurek,
     Lagrange,
     DiscontinuousLagrange,
     Serendipity,

src/interpolations.jl

@@ -21,6 +21,9 @@
 * `Lagrange{RefTriangle,5}`
 * `BubbleEnrichedLagrange{RefTriangle,1}`
 * `CrouzeixRaviart{RefTriangle, 1}`
+* `CrouzeixRaviart{RefTetrahedron, 1}`
+* `RannacherTurek{RefQuadrilateral, 1}`
+* `RannacherTurek{RefHexahedron, 1}`
 * `Lagrange{RefHexahedron,1}`
 * `Lagrange{RefHexahedron,2}`
 * `Lagrange{RefTetrahedron,1}`
@@ -477,6 +480,11 @@
 ############
 # Lagrange #
 ############
+"""
+    Lagrange{refshape, order} <: ScalarInterpolation
+
+Standard continuous Lagrange polynomials with equidistant node placement.
+"""
 struct Lagrange{shape, order, unused} <: ScalarInterpolation{shape, order}
     function Lagrange{shape, order}() where {shape <: AbstractRefShape, order}
         new{shape, order, Nothing}()
@@ -1298,6 +1306,11 @@
 ###############
 # Serendipity #
 ###############
+"""
+    Serendipity{refshape, order} <: ScalarInterpolation
+
+Serendipity element on hypercubes. Currently only second order variants are implemented.
+"""
 struct Serendipity{shape, order, unused} <: ScalarInterpolation{shape,order}
     function Serendipity{shape, order}() where {shape <: AbstractRefShape, order}
         new{shape, order, Nothing}()
@@ -1439,27 +1452,33 @@
 # Crouzeix–Raviart Elements #
 #############################
 """
+    CrouzeixRaviart{refshape, order} <: ScalarInterpolation
+
 Classical non-conforming Crouzeix–Raviart element.
 
-For details we refer to the original paper:
-M. Crouzeix and P. Raviart. "Conforming and nonconforming finite element
-methods for solving the stationary Stokes equations I." ESAIM: Mathematical Modelling
-and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 7.R3 (1973): 33-75.
+For details we refer to the original paper [CroRav:1973:cnf](@cite).
 """
 struct CrouzeixRaviart{shape, order, unused} <: ScalarInterpolation{shape, order}
     CrouzeixRaviart{RefTriangle, 1}() = new{RefTriangle, 1, Nothing}()
+    CrouzeixRaviart{RefTetrahedron, 1}() = new{RefTetrahedron, 1, Nothing}()
 end
 
+# CR elements are characterized by not having vertex dofs
+vertexdof_indices(ip::CrouzeixRaviart) = ntuple(i->(), nvertices(ip))
+
+#################################################
+# Non-conforming Crouzeix-Raviart dim 2 order 1 #
+#################################################
+getnbasefunctions(::CrouzeixRaviart{RefTriangle,1}) = 3
+
 adjust_dofs_during_distribution(::CrouzeixRaviart) = true
 adjust_dofs_during_distribution(::CrouzeixRaviart{<:Any, 1}) = false
 
-getnbasefunctions(::CrouzeixRaviart) = 3
-
-edgedof_indices(::CrouzeixRaviart) = ((1,), (2,), (3,))
-edgedof_interior_indices(::CrouzeixRaviart) = ((1,), (2,), (3,))
-facedof_indices(ip::CrouzeixRaviart) = (ntuple(i->i, getnbasefunctions(ip)),)
+edgedof_indices(::CrouzeixRaviart{RefTriangle,1}) = ((1,), (2,), (3,))
+edgedof_interior_indices(::CrouzeixRaviart{RefTriangle,1}) = ((1,), (2,), (3,))
+facedof_indices(ip::CrouzeixRaviart{RefTriangle,1}) = (ntuple(i->i, getnbasefunctions(ip)),)
 
-function reference_coordinates(::CrouzeixRaviart)
+function reference_coordinates(::CrouzeixRaviart{RefTriangle,1})
     return [Vec{2, Float64}((0.5, 0.5)),
             Vec{2, Float64}((0.0, 0.5)),
             Vec{2, Float64}((0.5, 0.0))]
@@ -1474,6 +1493,106 @@
     throw(ArgumentError("no shape function $i for interpolation $ip"))
 end
 
+#################################################
+# Non-conforming Crouzeix-Raviart dim 3 order 1 #
+#################################################
+getnbasefunctions(::CrouzeixRaviart{RefTetrahedron,1}) = 4
+
+facedof_indices(::CrouzeixRaviart{RefTetrahedron,1}) = ((1,), (2,), (3,), (4,))
+facedof_interior_indices(::CrouzeixRaviart{RefTetrahedron,1}) = ((1,), (2,), (3,), (4,))
+
+function reference_coordinates(::CrouzeixRaviart{RefTetrahedron,1})
+    return [
+            Vec{3, Float64}((1/3, 1/3, 0.0)),
+            Vec{3, Float64}((1/3, 0.0, 1/3)),
+            Vec{3, Float64}((1/3, 1/3, 1/3)),
+            Vec{3, Float64}((0.0, 1/3, 1/3)),
+            ]
+end
+
+function reference_shape_value(ip::CrouzeixRaviart{RefTetrahedron,1}, ξ::Vec{3}, i::Int)
+    (x,y,z) = ξ
+    i == 1 && return 1 -3z
+    i == 2 && return 1 -3y
+    i == 3 && return 3x +3y +3z -2
+    i == 4 && return 1 -3x
+    return throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+"""
+    RannacherTurek{refshape, order} <: ScalarInterpolation
+
+Classical non-conforming Rannacher-Turek element.
+
+This element is basically the idea from Crouzeix and Raviart applied to
+hypercubes. For details see the original paper [RanTur:1992:snq](@cite).
+"""
+struct RannacherTurek{shape,order} <: ScalarInterpolation{shape,order} end
+
+# CR-type elements are characterized by not having vertex dofs
+vertexdof_indices(ip::RannacherTurek) = ntuple(i->(), nvertices(ip))
+
+adjust_dofs_during_distribution(::RannacherTurek) = true
+adjust_dofs_during_distribution(::RannacherTurek{<:Any, 1}) = false
+
+#################################
+# Rannacher-Turek dim 2 order 1 #
+#################################
+getnbasefunctions(::RannacherTurek{RefQuadrilateral,1}) = 4
+
+edgedof_indices(::RannacherTurek{RefQuadrilateral,1}) = ((1,), (2,), (3,), (4,))
+edgedof_interior_indices(::RannacherTurek{RefQuadrilateral,1}) = ((1,), (2,), (3,), (4,))
+facedof_indices(ip::RannacherTurek{RefQuadrilateral,1}) = (ntuple(i->i, getnbasefunctions(ip)),)
+
+function reference_coordinates(::RannacherTurek{RefQuadrilateral,1})
+    return [Vec{2, Float64}(( 0.0, -1.0)),
+            Vec{2, Float64}(( 1.0,  0.0)),
+            Vec{2, Float64}(( 0.0,  1.0)),
+            Vec{2, Float64}((-1.0,  0.0))]
+end
+
+function reference_shape_value(ip::RannacherTurek{RefQuadrilateral,1}, ξ::Vec{2,T},  i::Int) where T
+    (x,y) = ξ
+
+    i == 1 && return -(x+1)^2/4 +(y+1)^2/4 +(x+1)/2 -(y+1)   +T(3)/4
+    i == 2 && return  (x+1)^2/4 -(y+1)^2/4          +(y+1)/2 -T(1)/4
+    i == 3 && return -(x+1)^2/4 +(y+1)^2/4 +(x+1)/2          -T(1)/4
+    i == 4 && return  (x+1)^2/4 -(y+1)^2/4 -(x+1)   +(y+1)/2 +T(3)/4
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+#################################
+# Rannacher-Turek dim 3 order 1 #
+#################################
+getnbasefunctions(::RannacherTurek{RefHexahedron,1}) = 6
+
+edgedof_indices(ip::RannacherTurek{RefHexahedron,1}) = ntuple(i->(), nedges(ip))
+edgedof_interior_indices(ip::RannacherTurek{RefHexahedron,1}) = ntuple(i->(), nedges(ip))
+facedof_indices(::RannacherTurek{RefHexahedron,1}) = ((1,), (2,), (3,), (4,), (5,), (6,))
+facedof_interior_indices(::RannacherTurek{RefHexahedron,1}) = ((1,), (2,), (3,), (4,), (5,), (6,))
+
+function reference_coordinates(::RannacherTurek{RefHexahedron,1})
+    return [Vec{3, Float64}(( 0.0,  0.0, -1.0)),
+            Vec{3, Float64}(( 0.0, -1.0,  0.0)),
+            Vec{3, Float64}(( 1.0,  0.0,  0.0)),
+            Vec{3, Float64}(( 0.0,  1.0,  0.0)),
+            Vec{3, Float64}((-1.0,  0.0,  0.0)),
+            Vec{3, Float64}(( 0.0,  0.0,  1.0)),]
+end
+
+function reference_shape_value(ip::RannacherTurek{RefHexahedron,1}, ξ::Vec{3,T}, i::Int) where T
+    (x,y,z) = ξ
+
+    i == 1 && return -2((x+1))^2/12 +1(x+1)/3 -2((y+1))^2/12 +1(y+1)/3 +4((z+1))^2/12 -7(z+1)/6 + T(2)/3
+    i == 2 && return -2((x+1))^2/12 +1(x+1)/3 +4((y+1))^2/12 -7(y+1)/6 -2((z+1))^2/12 +1(z+1)/3 + T(2)/3
+    i == 3 && return  4((x+1))^2/12 -1(x+1)/6 -2((y+1))^2/12 +1(y+1)/3 -2((z+1))^2/12 +1(z+1)/3 - T(1)/3
+    i == 4 && return -2((x+1))^2/12 +1(x+1)/3 +4((y+1))^2/12 -1(y+1)/6 -2((z+1))^2/12 +1(z+1)/3 - T(1)/3
+    i == 5 && return  4((x+1))^2/12 -7(x+1)/6 -2((y+1))^2/12 +1(y+1)/3 -2((z+1))^2/12 +1(z+1)/3 + T(2)/3
+    i == 6 && return -2((x+1))^2/12 +1(x+1)/3 -2((y+1))^2/12 +1(y+1)/3 +4((z+1))^2/12 -1(z+1)/6 - T(1)/3
+
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
 ##################################################
 # VectorizedInterpolation{<:ScalarInterpolation} #
 ##################################################

test/integration/test_simple_scalar_convergence.jl

@@ -14,19 +14,27 @@ get_geometry(::Ferrite.Interpolation{RefHexahedron}) = Hexahedron
 get_geometry(::Ferrite.Interpolation{RefTetrahedron}) = Tetrahedron
 get_geometry(::Ferrite.Interpolation{RefPyramid}) = Pyramid
 
-get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = 2*order
-get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = 2*order
-get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = 2*order+1
-get_quadrature_order(::BubbleEnrichedLagrange{shape, order}) where {shape, order} = 2*order
+get_quadrature_order(::Lagrange{shape, order}) where {shape, order} = max(2*order-1,2)
+get_quadrature_order(::Lagrange{RefPrism, order}) where order = 2*order # Don't know why
+get_quadrature_order(::Serendipity{shape, order}) where {shape, order} = max(2*order-1,2)
+get_quadrature_order(::CrouzeixRaviart{shape, order}) where {shape, order} = max(2*order-1,2)
+get_quadrature_order(::RannacherTurek{shape, order}) where {shape, order} = max(2*order-1,2)
+get_quadrature_order(::BubbleEnrichedLagrange{shape, order}) where {shape, order} = max(2*order-1,2)
 
 get_num_elements(::Ferrite.Interpolation{shape, 1}) where {shape} = 21
 get_num_elements(::Ferrite.Interpolation{shape, 2}) where {shape} = 7
 get_num_elements(::Ferrite.Interpolation{RefHexahedron, 1}) = 11
+get_num_elements(::Ferrite.RannacherTurek{RefQuadrilateral, 1}) = 15
+get_num_elements(::Ferrite.RannacherTurek{RefHexahedron, 1}) = 13
 get_num_elements(::Ferrite.Interpolation{RefHexahedron, 2}) = 4
 get_num_elements(::Ferrite.Interpolation{shape, 3}) where {shape} = 8
 get_num_elements(::Ferrite.Interpolation{shape, 4}) where {shape} = 5
 get_num_elements(::Ferrite.Interpolation{shape, 5}) where {shape} = 3
 
+get_test_tolerance(ip) = 1e-2
+get_test_tolerance(ip::RannacherTurek) = 4e-2
+get_test_tolerance(ip::CrouzeixRaviart) = 4e-2
+
 analytical_solution(x) = prod(cos, x*π/2)
 analytical_rhs(x) = -Tensors.laplace(analytical_solution,x)
 
@@ -140,7 +148,7 @@ end # module ConvergenceTestHelper
 
 # These test only for convergence within margins
 @testset "convergence analysis" begin
-    @testset "$interpolation" for interpolation in (
+    @testset failfast=true "$interpolation" for interpolation in (
         Lagrange{RefTriangle, 3}(),
         Lagrange{RefTriangle, 4}(),
         Lagrange{RefTriangle, 5}(),
@@ -154,7 +162,10 @@ end # module ConvergenceTestHelper
         #
         BubbleEnrichedLagrange{RefTriangle, 1}(),
         #
-        CrouzeixRaviart{RefTriangle, 1}(),
+        CrouzeixRaviart{RefTriangle,1}(),
+        CrouzeixRaviart{RefTetrahedron,1}(),
+        RannacherTurek{RefQuadrilateral,1}(),
+        RannacherTurek{RefHexahedron,1}(),
     )
         # Generate a grid ...
         geometry = ConvergenceTestHelper.get_geometry(interpolation)
@@ -173,7 +184,7 @@ end
 
 # These test also for correct convergence rates
 @testset "convergence rate" begin
-    @testset "$interpolation" for interpolation in (
+    @testset failfast=true  "$interpolation" for interpolation in (
         Lagrange{RefLine, 1}(),
         Lagrange{RefLine, 2}(),
         Lagrange{RefQuadrilateral, 1}(),
@@ -184,6 +195,10 @@ end
         Lagrange{RefHexahedron, 2}(),
         Lagrange{RefTetrahedron, 2}(),
         Lagrange{RefPrism, 2}(),
+        CrouzeixRaviart{RefTriangle,1}(),
+        CrouzeixRaviart{RefTetrahedron,1}(),
+        RannacherTurek{RefQuadrilateral,1}(),
+        RannacherTurek{RefHexahedron,1}(),
     )
         # Generate a grid ...
         geometry = ConvergenceTestHelper.get_geometry(interpolation)

test/test_interpolations.jl

@@ -40,7 +40,10 @@ using Ferrite: reference_shape_value, reference_shape_gradient
                       #
                       BubbleEnrichedLagrange{RefTriangle, 1}(),
                       #
-                      CrouzeixRaviart{RefTriangle, 1}(),
+                      CrouzeixRaviart{RefTriangle,1}(),
+                      CrouzeixRaviart{RefTetrahedron,1}(),
+                      RannacherTurek{RefQuadrilateral,1}(),
+                      RannacherTurek{RefHexahedron,1}(),
     )
         # Test of utility functions
         ref_dim = Ferrite.getrefdim(interpolation)

test/test_quadrules.jl

@@ -2,6 +2,7 @@ using Ferrite: reference_shape_value
 
 @testset "Quadrature testing" begin
     ref_tet_vol(dim) = 1 / factorial(dim)
+    ref_prism_vol() = 1 / 2
     ref_square_vol(dim) = 2^dim
 
     function integrate(qr::QuadratureRule, f::Function)
@@ -55,6 +56,13 @@ using Ferrite: reference_shape_value
     @test_throws ArgumentError QuadratureRule{RefTetrahedron}(:einstein, 2)
     @test_throws ArgumentError QuadratureRule{RefTetrahedron}(0)
 
+    g = (x) -> sqrt(x[1] + x[2])*x[3]^2
+    for order in 1:10
+        qr = QuadratureRule{RefPrism}(:polyquad, order)
+        @test integrate(qr, g) - 2/15 < 0.01
+        @test sum(qr.weights) ≈ ref_prism_vol()
+    end
+
     @testset "Quadrature rules for $ref_cell" for ref_cell in (
         Line,
         Quadrilateral,