Rework the AbstractCell interface

This patch reworks the `AbstractCell` interface (refer to #677 for
motivation and discussion for alternatives to this patch). Specifically,
this patch

 - introduces the reference dimension as a type parameter on `AbstractRefShape`,
 - introduces `RefInterval`, `RefTriangle`, `RefQuadrilateral`, and
   `RefHexahedron` as new subtypes of `AbstractRefShape{refdim}` (note
   that interpolations are not updated to this system yet since it can
   be done in a separate patch),
 - deprecates the old `Cell{refdim, nnodes, nfaces}` struct,
 - introduces new structs for all supported cell types, e.g.
   `struct Triangle <: AbstractCell{2, RefTriangle}`,
 - defines the unique vertex/edge/face identifiers for DoF-distribution
   in terms of the reference shape.

The new parametrization makes it possible to dispatch on the reference
shape instead of union of concrete implementations, i.e.
`::AbstractCell{2, RefTriangle}` instead of `::Union{Triangle,
QuadraticTriangle}`, and similarly for the reference dimension.

One drawback is that it is no longer possible to use "anonymous
celltypes" after this patch, i.e. using, e.g., `Cell{3, 20, 6}` and
defining methods for that instead of explicitly naming it
`SerendipityQuadraticHexahedron` or similar. Instead, you now have to
define a struct and some methods. Arguably this is a good thing --
nobody understands the meaning of `Cell{3, 20, 6}`.

For normal usage this patch is not breaking (see e.g. no required
changes for examples), unless one dispatches on the `Cell` struct
directly.

Fixes #677.

src/Export/VTK.jl

@@ -1,17 +1,14 @@
 cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE
-cell_to_vtkcell(::Type{Line2D}) = VTKCellTypes.VTK_LINE
-cell_to_vtkcell(::Type{Line3D}) = VTKCellTypes.VTK_LINE
 cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE
 
 cell_to_vtkcell(::Type{Quadrilateral}) = VTKCellTypes.VTK_QUAD
-cell_to_vtkcell(::Type{Quadrilateral3D}) = VTKCellTypes.VTK_QUAD
 cell_to_vtkcell(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
 cell_to_vtkcell(::Type{Triangle}) = VTKCellTypes.VTK_TRIANGLE
 cell_to_vtkcell(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE
-cell_to_vtkcell(::Type{Cell{2,8,4}}) = VTKCellTypes.VTK_QUADRATIC_QUAD
+cell_to_vtkcell(::Type{SQuadraticQuadrilateral}) = VTKCellTypes.VTK_QUADRATIC_QUAD
 
 cell_to_vtkcell(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON
-cell_to_vtkcell(::Type{Cell{3,20,6}}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON
+cell_to_vtkcell(::Type{SQuadraticHexahedron}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON
 cell_to_vtkcell(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA
 cell_to_vtkcell(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA
 cell_to_vtkcell(::Type{Wedge}) = VTKCellTypes.VTK_WEDGE

src/Ferrite.jl

@@ -16,11 +16,17 @@ Represents a reference shape which quadrature rules and interpolations are defin
 Currently, the only concrete types that subtype this type are `RefCube` in 1, 2 and 3 dimensions,
 and `RefTetrahedron` in 2 and 3 dimensions.
 """
-abstract type AbstractRefShape end
+abstract type AbstractRefShape{refdim} end
 
-struct RefTetrahedron <: AbstractRefShape end
-struct RefCube <: AbstractRefShape end
-struct RefPrism <: AbstractRefShape end
+struct RefInterval      <: AbstractRefShape{1} end
+struct RefTriangle      <: AbstractRefShape{2} end
+struct RefQuadrilateral <: AbstractRefShape{2} end
+struct RefTetrahedron   <: AbstractRefShape{3} end
+struct RefHexahedron    <: AbstractRefShape{3} end
+struct RefPrism         <: AbstractRefShape{3} end
+
+# TODO: Update interpolation definitions and enable deprecation
+const RefCube = RefHexahedron
 
 """
 Abstract type which has `CellValues` and `FaceValues` as subtypes

src/Grid/grid_generators.jl

@@ -102,13 +102,12 @@ function _generate_2d_nodes!(nodes, nx, ny, LL, LR, UR, UL)
     end
 end
 
-
-function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, X::Vector{Vec{2,T}}) where {M,N,T}
+function generate_grid(C::Type{<:AbstractCell{2}}, nel::NTuple{2,Int}, X::Vector{Vec{2,T}}) where {T}
     @assert length(X) == 4
     generate_grid(C, nel, X[1], X[2], X[3], X[4])
 end
 
-function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, left::Vec{2,T}=Vec{2}((-1.0,-1.0)), right::Vec{2,T}=Vec{2}((1.0,1.0))) where {M,N,T}
+function generate_grid(C::Type{<:AbstractCell{2}}, nel::NTuple{2,Int}, left::Vec{2,T}=Vec{2}((-1.0,-1.0)), right::Vec{2,T}=Vec{2}((1.0,1.0))) where {T}
     LL = left
     UR = right
     LR = Vec{2}((UR[1], LL[2]))
@@ -288,7 +287,7 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}(
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end 
 
-function Ferrite.generate_grid(::Type{Cell{3,20,6}}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
+function Ferrite.generate_grid(::Type{SQuadraticHexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
     nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
     nnode_x = 2nel_x + 1; nnode_y = 2nel_y + 1; nnode_z = 2nel_z + 1 #Note: not the actually number of nodes in x/y/z, just a temporary variables
 
@@ -311,9 +310,9 @@ function Ferrite.generate_grid(::Type{Cell{3,20,6}}, nel::NTuple{3,Int}, left::V
 
 
     # Generate cells
-    cells = Cell{3,20,6}[]
+    cells = SQuadraticHexahedron[]
     for k in 1:2:2nel_z, j in 1:2:2nel_y, i in 1:2:2nel_x     
-        push!(cells, Cell{3,20,6}((
+        push!(cells, SQuadraticHexahedron((
                 node_array[i,j,k], node_array[i+2,j,k], node_array[i+2,j+2,k], node_array[i,j+2,k], # vertices bot 
                 node_array[i,j,k+2], node_array[i+2,j,k+2], node_array[i+2,j+2,k+2], node_array[i,j+2,k+2], # vertices top
                 node_array[i+1,j,k], node_array[i+2,j+1,k], node_array[i+1,j+2,k], node_array[i,j+1,k], # edges horizontal bottom

src/deprecations.jl

@@ -10,3 +10,45 @@ import Base: push!
 @deprecate nfields(dh::AbstractDofHandler) length(getfieldnames(dh)) false
 
 @deprecate add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet) add!(ch, dbc)
+
+
+struct Cell{refdim, nnodes, nfaces}
+    function Cell{refdim, nnodes, nfaces}(nodes) where {refdim, nnodes, nfaces}
+        params = (refdim, nnodes, nfaces)
+        replacement = nothing
+        if params == (1, 2, 2) || params == (2, 2, 1) || params == (3, 2, 0)
+            replacement = Line
+        elseif params == (1, 3, 2)
+            replacement = QuadraticLine
+        elseif params == (2, 3, 3)
+            replacement = Triangle
+        elseif params == (2, 6, 3)
+            replacement = QuadraticTriangle
+        elseif params == (2, 4, 4) || params == (3, 4, 1)
+            replacement = Quadrilateral
+        elseif params == (2, 9, 4)
+            replacement = QuadraticQuadrilateral
+        elseif params == (2, 8, 4)
+            replacement = SQuadraticQuadrilateral
+        elseif params == (3, 4, 4)
+            replacement = Tetrahedron
+        elseif params == (3, 10, 4)
+            replacement = QuadraticTetrahedron
+        elseif params == (3, 8, 6)
+            replacement = Hexahedron
+        elseif params == (3, 27, 6)
+            replacement = QuadraticHexahedron
+        elseif params == (3, 20, 6)
+            replacement = SQuadraticHexahedron
+        elseif params == (3, 6, 5)
+            replacement = Wedge
+        end
+        error("Fool of a Took!")
+        if replacement === nothing
+            error("?")
+        else
+            Base.depwarn("Use `$(replacement)(nodes)` instead of `Cell{$refdim, $nnodes, $nfaces}(nodes)`.", :Cell)
+            return replacement(nodes)
+        end
+    end
+end

src/exports.jl

@@ -52,10 +52,12 @@ export
     Quadrilateral,
     Quadrilateral3D,
     QuadraticQuadrilateral,
+    SQuadraticQuadrilateral,
     Tetrahedron,
     QuadraticTetrahedron,
     Hexahedron,
-    #QuadraticHexahedron,
+    QuadraticHexahedron,
+    SQuadraticHexahedron,
     Wedge,
     CellIndex,
     FaceIndex,

src/iterators.jl

@@ -97,7 +97,7 @@ cellid(cc::CellCache) = cc.cellid[]
 celldofs!(v::Vector, cc::CellCache) = copyto!(v, cc.dofs) # celldofs!(v, cc.dh, cc.cellid[])
 
 # TODO: These should really be replaced with something better...
-nfaces(cc::CellCache) = nfaces(eltype(cc.grid.cells))
+nfaces(cc::CellCache) = nfaces(cc.grid.cells[cc.cellid[]])
 onboundary(cc::CellCache, face::Int) = cc.grid.boundary_matrix[face, cc.cellid[]]
 
 ##################

test/test_apply_analytical.jl

@@ -11,7 +11,7 @@
         return B{Dim,RefShape,order}()
     end
     getcellorder(CT) = Ferrite.getorder(Ferrite.default_interpolation(CT))
-    getcelltypedim(::Type{<:Cell{dim}}) where dim = dim
+    getcelltypedim(::Type{<:Ferrite.AbstractCell{dim}}) where dim = dim
 
     # Functions to create dof handlers for testing
     function testdh(CT, ip_order_u, ip_order_p)
@@ -40,14 +40,14 @@
     function testdh_subdomains(dim, ip_order_u, ip_order_p)
         if dim == 1
             nodes = Node.([Vec{1}((x,)) for x in 0.0:1.0:3.0])
-            cell1 = Cell{1,2,2}((1,2))
-            cell2 = Cell{1,3,2}((2,3,4))
+            cell1 = Line((1,2))
+            cell2 = QuadraticLine((2,3,4))
             grid = Grid([cell1, cell2], nodes; cellsets=Dict("A"=>Set(1:1), "B"=>Set(2:2)))
         elseif dim == 2
             nodes = Node.([Vec{2}((x,y)) for y in 0.0:2 for x in 0.0:1])
-            cell1 = Cell{2,3,3}((1,2,3))
-            cell2 = Cell{2,3,3}((2,4,3))
-            cell3 = Cell{2,4,4}((3,4,6,5))
+            cell1 = Triangle((1,2,3))
+            cell2 = Triangle((2,4,3))
+            cell3 = Quadrilateral((3,4,6,5))
             grid = Grid([cell1,cell2,cell3], nodes; cellsets=Dict("A"=>Set(1:2), "B"=>Set(3:3)))
         else
             error("Only dim=1 & 2 supported")
@@ -88,7 +88,14 @@
     end
 
     @testset "DofHandler" begin
-        for CT in Ferrite.implemented_celltypes
+        for CT in (
+            Line, QuadraticLine,
+            Triangle, QuadraticTriangle,
+            Quadrilateral, QuadraticQuadrilateral,
+            Tetrahedron, QuadraticTetrahedron,
+            Hexahedron, # SQuadraticHexahedron,
+            Wedge,
+        )
             for ip_order_u in 1:2
                 for ip_order_p in 1:2
                     dh = testdh(CT, ip_order_u, ip_order_p)

test/test_grid_dofhandler_vtk.jl

@@ -16,7 +16,7 @@ end
                             (Triangle,               2),
                             (QuadraticTriangle,      2),
                             (Hexahedron,             3),
-                            (Cell{3,20,6},           3),
+                            (SQuadraticHexahedron,   3),
                             (Tetrahedron,            3))
 
         # create test grid, do some operations on it and then test
@@ -284,7 +284,7 @@ end
     @test getneighborhood(topology,hexgrid,FaceIndex((3,3))) == [FaceIndex((4,5))]
     @test getneighborhood(topology,hexgrid,FaceIndex((4,2))) == [FaceIndex((2,4))]
     @test getneighborhood(topology,hexgrid,FaceIndex((4,5))) == [FaceIndex((3,3))]
-    serendipitygrid = generate_grid(Cell{3,20,6},(2,2,1))
+    serendipitygrid = generate_grid(SQuadraticHexahedron,(2,2,1))
     stopology = ExclusiveTopology(serendipitygrid)
     # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
     @test all(stopology.face_neighbor .== topology.face_neighbor)