Use the new reference shapes for Interpolations and QuadratureRule (

#711)

This patch replaces `RefCube{1, 2, 3}` with `RefLine`,
`RefQuadrilateral`, and `RefHexahedron`; and `RefTetrahedron{2}` with
`RefTriangle` in the definitions of interpolations and quadrature.

Since the reference dimension is now encoded in the shape, this patch
also changes the parameter order of interpolations, to move the dim
parameter to the end, and always set it to `Nothing`. It can't be
removed completely right away since that would give cryptic error
messages.

Interpolations are thus now constructed with just the reference shape
and the order, e.g. `Lagrange{RefTriangle, 2}()` instead of `Lagrange{2,
RefTetrahedron, 2}()`.

The same is done for `QuadratureRule`, for constructing cell quadrature
rules. A face quadrature rule is still written as `QuadratureRule{D-1,
shape}(...)`, where `shape` is the reference shape of the *cell*, i.e.
face quadrature for the hexahedron is constructed as `QuadratureRule{2,
RefHexahedron}()`. The new reference shapes make this non-ambiguous
since we can now directly that a face quadrature is constructed.

docs/src/literate/computational_homogenization.jl

@@ -212,8 +212,8 @@ grid = togrid("periodic-rve.msh") #src
 # cellvalues as usual:
 
 dim = 2
-ip = Lagrange{dim, RefTetrahedron, 1}()^dim
-qr = QuadratureRule{dim, RefTetrahedron}(2)
+ip = Lagrange{RefTriangle, 1}()^dim
+qr = QuadratureRule{dim, RefTriangle}(2)
 cellvalues = CellValues(qr, ip);
 
 # We define a dof handler with a displacement field `:u`:

docs/src/literate/heat_equation.jl

@@ -52,12 +52,12 @@ grid = generate_grid(Quadrilateral, (20, 20));
 # test and trial functions (among other things). To define
 # this we need to specify an interpolation space for the shape functions.
 # We use Lagrange functions
-# based on the two-dimensional reference "cube". We also define a quadrature rule based on
+# based on the two-dimensional reference quadrilateral. We also define a quadrature rule based on
 # the same reference cube. We combine the interpolation and the quadrature rule
 # to a `CellValues` object.
 dim = 2
-ip = Lagrange{dim, RefCube, 1}()
-qr = QuadratureRule{dim, RefCube}(2)
+ip = Lagrange{RefQuadrilateral, 1}()
+qr = QuadratureRule{dim, RefQuadrilateral}(2)
 cellvalues = CellValues(qr, ip);
 
 # ### Degrees of freedom

docs/src/literate/helmholtz.jl

@@ -54,9 +54,9 @@ const Δ = Tensors.hessian;
 grid = generate_grid(Quadrilateral, (150, 150))
 
 dim = 2
-ip = Lagrange{dim, RefCube, 1}()
-qr = QuadratureRule{dim, RefCube}(2)
-qr_face = QuadratureRule{dim-1, RefCube}(2)
+ip = Lagrange{RefQuadrilateral, 1}()
+qr = QuadratureRule{dim, RefQuadrilateral}(2)
+qr_face = QuadratureRule{dim-1, RefQuadrilateral}(2)
 cellvalues = CellValues(qr, ip);
 facevalues = FaceValues(qr_face, ip);
 

docs/src/literate/hyperelasticity.jl

@@ -315,7 +315,7 @@ function solve()
     mp = NeoHooke(μ, λ)
 
     ## Finite element base
-    ip = Lagrange{3, RefTetrahedron, 1}()^3
+    ip = Lagrange{RefTetrahedron, 1}()^3
     qr = QuadratureRule{3, RefTetrahedron}(1)
     qr_face = QuadratureRule{2, RefTetrahedron}(1)
     cv = CellValues(qr, ip)

docs/src/literate/incompressible_elasticity.jl

@@ -41,8 +41,8 @@ end;
 # Next we define a function to set up our cell- and facevalues.
 function create_values(interpolation_u, interpolation_p)
     ## quadrature rules
-    qr      = QuadratureRule{2,RefTetrahedron}(3)
-    face_qr = QuadratureRule{1,RefTetrahedron}(3)
+    qr      = QuadratureRule{2,RefTriangle}(3)
+    face_qr = QuadratureRule{1,RefTriangle}(3)
 
     ## cell and facevalues for u
     cellvalues_u = CellValues(qr, interpolation_u)
@@ -207,7 +207,7 @@ function solve(ν, interpolation_u, interpolation_p)
     u = Symmetric(K) \ f;
 
     ## export
-    filename = "cook_" * (isa(interpolation_u, Lagrange{2,RefTetrahedron,1}) ? "linear" : "quadratic") *
+    filename = "cook_" * (isa(interpolation_u, Lagrange{RefTriangle,1}) ? "linear" : "quadratic") *
                          "_linear"
     vtk_grid(filename, dh) do vtkfile
         vtk_point_data(vtkfile, dh, u)
@@ -220,9 +220,9 @@ end
 # vectorize it to 2 dimensions such that we obtain vector shape functions (and 2nd order
 # tensors for the gradients).
 
-linear_p    = Lagrange{2,RefTetrahedron,1}()
-linear_u    = Lagrange{2,RefTetrahedron,1}()^2
-quadratic_u = Lagrange{2,RefTetrahedron,2}()^2
+linear_p    = Lagrange{RefTriangle,1}()
+linear_u    = Lagrange{RefTriangle,1}()^2
+quadratic_u = Lagrange{RefTriangle,2}()^2
 
 # All that is left is to solve the problem. We choose a value of Poissons
 # ratio that is near incompressibility -- $ν = 0.5$ -- and thus expect the

docs/src/literate/landau.jl

@@ -87,7 +87,7 @@ function LandauModel(α, G, gridsize, left::Vec{DIM, T}, right::Vec{DIM, T}, elp
     threadindices = Ferrite.create_coloring(grid)
 
     qr  = QuadratureRule{DIM, RefTetrahedron}(2)
-    ipP = Lagrange{DIM, RefTetrahedron, 1}()^3
+    ipP = Lagrange{RefTetrahedron, 1}()^3
     cvP = CellValues(qr, ipP)
 
     dofhandler = DofHandler(grid)

docs/src/literate/linear_shell.jl

@@ -26,9 +26,9 @@ grid = generate_shell_grid(nels, size)
 # We also create two quadrature rules for the in-plane and out-of-plane directions. Note that we use 
 # under integration for the inplane integration, to avoid shear locking. 
 #+
-ip = Lagrange{2,RefCube,1}()
-qr_inplane = QuadratureRule{2,RefCube}(1)
-qr_ooplane = QuadratureRule{1,RefCube}(2)
+ip = Lagrange{RefQuadrilateral,1}()
+qr_inplane = QuadratureRule{2,RefQuadrilateral}(1)
+qr_ooplane = QuadratureRule{1,RefLine}(2)
 cv = CellValues(qr_inplane, ip)
 
 # Next we distribute displacement dofs,`:u = (x,y,z)` and rotational dofs, `:θ = (θ₁,  θ₂)`.

docs/src/literate/ns_vs_diffeq.jl

@@ -158,11 +158,11 @@ grid = generate_grid(Quadrilateral, (x_cells, y_cells), Vec{2}((0.0, 0.0)), Vec{
 # To ensure stability we utilize the Taylor-Hood element pair Q2-Q1.
 # We have to utilize the same quadrature rule for the pressure as for the velocity, because in the weak form the
 # linear pressure term is tested against a quadratic function.
-ip_v = Lagrange{dim, RefCube, 2}()^dim
-qr = QuadratureRule{dim, RefCube}(4)
+ip_v = Lagrange{RefQuadrilateral, 2}()^dim
+qr = QuadratureRule{dim, RefQuadrilateral}(4)
 cellvalues_v = CellValues(qr, ip_v);
 
-ip_p = Lagrange{dim, RefCube, 1}()
+ip_p = Lagrange{RefQuadrilateral, 1}()
 cellvalues_p = CellValues(qr, ip_p);
 
 dh = DofHandler(grid)

docs/src/literate/plasticity.jl

@@ -269,7 +269,7 @@ function solve()
     P1 = Vec((0.0, 0.0, 0.0))  # start point for geometry
     P2 = Vec((L, w, h))        # end point for geometry
     grid = generate_grid(Tetrahedron, nels, P1, P2)
-    interpolation = Lagrange{3, RefTetrahedron, 1}()^3
+    interpolation = Lagrange{RefTetrahedron, 1}()^3
 
     dh = create_dofhandler(grid, interpolation) # JuaFEM helper function
     dbcs = create_bc(dh, grid) # create Dirichlet boundary-conditions

docs/src/literate/quasi_incompressible_hyperelasticity.jl

@@ -347,8 +347,8 @@ function solve(interpolation_u, interpolation_p)
 end;
 
 # We can now test the solution using the Taylor-Hood approximation
-quadratic_u = Lagrange{3, RefTetrahedron, 2}()^3
-linear_p = Lagrange{3, RefTetrahedron, 1}()
+quadratic_u = Lagrange{RefTetrahedron, 2}()^3
+linear_p = Lagrange{RefTetrahedron, 1}()
 vol_def = solve(quadratic_u, linear_p)
 
 # The deformed volume is indeed close to 1 (as should be for a nearly incompressible material).

docs/src/literate/stokes-flow.jl

@@ -227,10 +227,10 @@ end
 # the boundary when assembling the constraint matrix ``\underline{\underline{C}}``.
 
 function setup_fevalues(ipu, ipp, ipg)
-    qr = QuadratureRule{2,RefTetrahedron}(2)
+    qr = QuadratureRule{2,RefTriangle}(2)
     cvu = CellValues(qr, ipu, ipg)
     cvp = CellValues(qr, ipp, ipg)
-    qr_face = QuadratureRule{1,RefTetrahedron}(2)
+    qr_face = QuadratureRule{1,RefTriangle}(2)
     fvp = FaceValues(qr_face, ipp, ipg)
     return cvu, cvp, fvp
 end
@@ -485,12 +485,12 @@ function main()
     h = 0.05 # approximate element size
     grid = setup_grid(h)
     ## Interpolations
-    ipu = Lagrange{2,RefTetrahedron,2}() ^ 2 # quadratic
-    ipp = Lagrange{2,RefTetrahedron,1}()     # linear
+    ipu = Lagrange{RefTriangle,2}() ^ 2 # quadratic
+    ipp = Lagrange{RefTriangle,1}()     # linear
     ## Dofs
     dh = setup_dofs(grid, ipu, ipp)
     ## FE values
-    ipg = Lagrange{2,RefTetrahedron,1}() # linear geometric interpolation
+    ipg = Lagrange{RefTriangle,1}() # linear geometric interpolation
     cvu, cvp, fvp = setup_fevalues(ipu, ipp, ipg)
     ## Boundary conditions
     ch = setup_constraints(dh, fvp)

docs/src/literate/threaded_assembly.jl

@@ -86,7 +86,7 @@ end;
 
 # Each thread need its own CellValues and FaceValues (although, for this example we don't use
 # the FaceValues)
-function create_values(interpolation_space::Interpolation{dim, refshape}, qr_order::Int) where {dim, refshape}
+function create_values(interpolation_space::Interpolation{refshape}, qr_order::Int) where {dim, refshape<:Ferrite.AbstractRefShape{dim}}
     ## Interpolations and values
     quadrature_rule = QuadratureRule{dim, refshape}(qr_order)
     face_quadrature_rule = QuadratureRule{dim-1, refshape}(qr_order)
@@ -179,7 +179,7 @@ end;
 function run_assemble()
     n = 20
     grid, colors = create_colored_cantilever_grid(Hexahedron, n);
-    ip = Lagrange{3,RefCube,1}()^3
+    ip = Lagrange{RefHexahedron,1}()^3
     dh = create_dofhandler(grid, ip);
 
     K = create_sparsity_pattern(dh);

docs/src/literate/topology_optimization.jl

@@ -88,11 +88,11 @@ end
 
 function create_values()
     ## quadrature rules
-    qr      = QuadratureRule{2,RefCube}(2)
-    face_qr = QuadratureRule{1,RefCube}(2)
+    qr      = QuadratureRule{2,RefQuadrilateral}(2)
+    face_qr = QuadratureRule{1,RefQuadrilateral}(2)
 
     ## cell and facevalues for u
-    ip = Lagrange{2,RefCube,1}()^2
+    ip = Lagrange{RefQuadrilateral,1}()^2
     cellvalues = CellValues(qr, ip)
     facevalues = FaceValues(face_qr, ip)
 
@@ -101,7 +101,7 @@ end
 
 function create_dofhandler(grid)
     dh = DofHandler(grid)
-    add!(dh, :u, Lagrange{2,RefCube,1}()^2) # displacement
+    add!(dh, :u, Lagrange{RefQuadrilateral,1}()^2) # displacement
     close!(dh)
     return dh
 end

docs/src/literate/transient_heat_equation.jl

@@ -65,8 +65,8 @@ grid = generate_grid(Quadrilateral, (100, 100));
 # ### Trial and test functions
 # Again, we define the structs that are responsible for the `shape_value` and `shape_gradient` evaluation.
 dim = 2
-ip = Lagrange{dim, RefCube, 1}()
-qr = QuadratureRule{dim, RefCube}(2)
+ip = Lagrange{RefQuadrilateral, 1}()
+qr = QuadratureRule{dim, RefQuadrilateral}(2)
 cellvalues = CellValues(qr, ip);
 
 # ### Degrees of freedom

src/Dofs/ConstraintHandler.jl

@@ -1167,10 +1167,10 @@ function construct_cornerish(min_x::V, max_x::V) where {T, V <: Vec{3,T}}
     ]
 end
 
-function mirror_local_dofs(_, _, ::Lagrange{1}, ::Int)
+function mirror_local_dofs(_, _, ::Lagrange{RefLine}, ::Int)
     # For 1D there is nothing to do
 end
-function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{2,<:Union{RefCube,RefTetrahedron}}, n::Int)
+function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{<:Union{RefQuadrilateral,RefTriangle}}, n::Int)
     # For 2D we always permute since Ferrite defines dofs counter-clockwise
     ret = collect(1:length(local_face_dofs))
     for (i, f) in enumerate(facedof_indices(ip))
@@ -1188,9 +1188,9 @@ function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange
 end
 
 # TODO: Can probably be combined with the method above.
-function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{3,<:Union{RefCube,RefTetrahedron},O}, n::Int) where O
+function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{<:Union{RefHexahedron,RefTetrahedron},O}, n::Int) where O
     @assert 1 <= O <= 2
-    N = ip isa Lagrange{3,RefCube} ? 4 : 3
+    N = ip isa Lagrange{RefHexahedron} ? 4 : 3
     ret = collect(1:length(local_face_dofs))
 
     # Mirror by changing from counter-clockwise to clockwise
@@ -1234,12 +1234,12 @@ end
 circshift!(args...) = Base.circshift!(args...)
 
 
-function rotate_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{2}, ncomponents)
+function rotate_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{<:Union{RefQuadrilateral,RefTriangle}}, ncomponents)
     return collect(1:length(local_face_dofs)) # TODO: Return range?
 end
-function rotate_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{3,<:Union{RefCube,RefTetrahedron}, O}, ncomponents) where O
+function rotate_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{<:Union{RefHexahedron,RefTetrahedron}, O}, ncomponents) where O
     @assert 1 <= O <= 2
-    N = ip isa Lagrange{3,RefCube} ? 4 : 3
+    N = ip isa Lagrange{RefHexahedron} ? 4 : 3
     ret = similar(local_face_dofs, length(local_face_dofs), N)
     ret[:, :] .= 1:length(local_face_dofs)
     for f in 1:length(local_face_dofs_offset)-1

src/Dofs/apply_analytical.jl

@@ -50,7 +50,7 @@ end
 
 function _apply_analytical!(
     a::AbstractVector, dh::AbstractDofHandler, celldofinds, field_dim,
-    ip_fun::Interpolation{dim,RefShape}, ip_geo::Interpolation, f::Function, cellset) where {dim, RefShape}
+    ip_fun::Interpolation{RefShape}, ip_geo::Interpolation, f::Function, cellset) where {dim, RefShape<:AbstractRefShape{dim}}
 
     coords = getcoordinates(dh.grid, first(cellset))
     ref_points = reference_coordinates(ip_fun)

src/FEValues/cell_values.jl

@@ -29,8 +29,8 @@ values of nodal functions, gradients and divergences of nodal functions etc. in
 """
 CellValues
 
-function default_geometric_interpolation(::Interpolation{dim,shape}) where {dim, shape}
-    return Lagrange{dim,shape,1}()
+function default_geometric_interpolation(::Interpolation{shape}) where {shape}
+    return Lagrange{shape,1}()
 end
 
 struct CellValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, GIP} <: AbstractCellValues
@@ -52,10 +52,10 @@ function CellValues(qr::QuadratureRule, ip::Interpolation,
 end
 # TODO: This doesn't actually work for T != Float64
 function CellValues(::Type{T}, qr::QR, ip::IP, gip::GIP = default_geometric_interpolation(ip)) where {
-    dim, shape, T,
+    dim, shape <: AbstractRefShape{dim}, T,
     QR  <: QuadratureRule{dim, shape},
-    IP  <: Union{ScalarInterpolation{dim, shape}, VectorInterpolation{dim, dim, shape}},
-    GIP <: ScalarInterpolation{dim, shape}
+    IP  <: Union{ScalarInterpolation{shape}, VectorInterpolation{dim, shape}},
+    GIP <: ScalarInterpolation{shape}
 }
     n_qpoints = length(getweights(qr))