From c0be5344a487adad123bb2c24d14f49af8e1798c Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Tue, 18 Apr 2023 14:59:04 +0200 Subject: [PATCH] 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 | 7 +- src/Ferrite.jl | 14 +- src/Grid/grid.jl | 274 ++++++++++++++++++++----------- src/Grid/grid_generators.jl | 11 +- src/deprecations.jl | 42 +++++ src/exports.jl | 4 +- src/iterators.jl | 2 +- test/test_apply_analytical.jl | 21 ++- test/test_grid_dofhandler_vtk.jl | 4 +- 9 files changed, 260 insertions(+), 119 deletions(-) diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 9a1cf6de6e..d7f36c8a64 100644 --- a/src/Export/VTK.jl +++ b/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 diff --git a/src/Ferrite.jl b/src/Ferrite.jl index bb378b8a2b..4570b3b329 100644 --- a/src/Ferrite.jl +++ b/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 diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index cb1665545d..066f06e2b2 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -22,27 +22,16 @@ Get the data type of the components of the nodes coordinate. """ get_coordinate_eltype(::Node{dim,T}) where {dim,T} = T -abstract type AbstractCell{dim,N,M} end -""" - Cell{dim,N,M} <: AbstractCell{dim,N,M} +########################## +# AbstractCell interface # +########################## -A `Cell` is a subdomain defined by a collection of `Node`s. -The parameter `dim` refers here to the geometrical/ambient dimension, i.e. the dimension of the `nodes` in the grid and **not** the topological dimension of the cell (i.e. the dimension of the reference element obtained by default_interpolation). -A `Cell` has `N` nodes and `M` faces. -Note that a `Cell` is not defined geometrically by node coordinates, but rather topologically by node indices into the node vector of some grid. +abstract type AbstractCell{refdim, refshape <: AbstractRefShape{refdim}} end -# Fields -- `nodes::Ntuple{N,Int}`: N-tuple that stores the node ids. The ordering defines a cell's and its subentities' orientations. -""" -struct Cell{dim,N,M} <: AbstractCell{dim,N,M} - nodes::NTuple{N,Int} -end -nfaces(c::C) where {C<:AbstractCell} = nfaces(typeof(c)) -nfaces(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = M -nedges(c::C) where {C<:AbstractCell} = length(edges(c)) -nvertices(c::C) where {C<:AbstractCell} = length(vertices(c)) -nnodes(c::C) where {C<:AbstractCell} = nnodes(typeof(c)) -nnodes(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = N +nvertices(c::AbstractCell) = length(vertices(c)) +nedges( c::AbstractCell) = length(edges(c)) +nfaces( c::AbstractCell) = length(faces(c)) +nnodes( c::AbstractCell) = length(get_node_ids(c)) """ Ferrite.vertices(::AbstractCell) @@ -51,7 +40,7 @@ Returns a tuple with the node indices (of the nodes in a grid) for each vertex i This function induces the [`VertexIndex`](@ref), where the second index corresponds to the local index into this tuple. """ -vertices(::Ferrite.AbstractCell) +vertices(::AbstractCell) """ Ferrite.edges(::AbstractCell) @@ -62,7 +51,7 @@ the vertices that define an *oriented edge*. This function induces the Note that the vertices are sufficient to define an edge uniquely. """ -edges(::Ferrite.AbstractCell) +edges(::AbstractCell) """ Ferrite.faces(::AbstractCell) @@ -73,37 +62,176 @@ the vertices that define an *oriented face*. This function induces the Note that the vertices are sufficient to define a face uniquely. """ -faces(::Ferrite.AbstractCell) +faces(::AbstractCell) """ Ferrite.default_interpolation(::AbstractCell)::Interpolation Returns the interpolation which defines the geometry of a given cell. """ -default_interpolation(::Ferrite.AbstractCell) +default_interpolation(::AbstractCell) -# Typealias for commonly used cells -const implemented_celltypes = ( - (const Line = Cell{1,2,2}), - (const Line2D = Cell{2,2,1}), - (const Line3D = Cell{3,2,0}), - (const QuadraticLine = Cell{1,3,2}), - - (const Triangle = Cell{2,3,3}), - (const QuadraticTriangle = Cell{2,6,3}), - - (const Quadrilateral = Cell{2,4,4}), - (const Quadrilateral3D = Cell{3,4,1}), - (const QuadraticQuadrilateral = Cell{2,9,4}), - - (const Tetrahedron = Cell{3,4,4}), - (const QuadraticTetrahedron = Cell{3,10,4}), - - (const Hexahedron = Cell{3,8,6}), - (Cell{2,20,6}), +""" + Ferrite.get_node_ids(c::AbstractCell) - (const Wedge = Cell{3,6,5}) -) +Return the node id's for cell `c` in the order determined by the cell's reference cell. + +Default implementation: `c.nodes`. +""" +get_node_ids(c::AbstractCell) = c.nodes + +# Default implementations of vertices/edges/faces that work as long as get_node_ids is +# correctly implemented for the cell. + +# RefInterval (refdim = 1): vertices for vertexdofs, faces for BC +function vertices(c::AbstractCell{1, RefInterval}) + ns = get_node_ids(c) + return (ns[1], ns[2]) # v1, v2 +end +function faces(c::AbstractCell{1, RefInterval}) + ns = get_node_ids(c) + return (ns[1], ns[2]) # f1, f2 +end + +# RefTriangle (refdim = 2): vertices for vertexdofs, faces for facedofs (edgedofs) and BC +function vertices(c::AbstractCell{2, RefTriangle}) + ns = get_node_ids(c) + return (ns[1], ns[2], ns[3]) # v1, v2, v3 +end +function faces(c::AbstractCell{2, RefTriangle}) + ns = get_node_ids(c) + return ( + (ns[1], ns[2]), (ns[2], ns[3]), (ns[3], ns[1]), # f1, f2, f3 + ) +end + +# RefQuadrilateral (refdim = 2): vertices for vertexdofs, faces for facedofs (edgedofs) and BC +function vertices(c::AbstractCell{2, RefQuadrilateral}) + ns = get_node_ids(c) + return (ns[1], ns[2], ns[3], ns[4]) # v1, v2, v3, v4 +end +function faces(c::AbstractCell{2, RefQuadrilateral}) + ns = get_node_ids(c) + return ( + (ns[1], ns[2]), (ns[2], ns[3]), (ns[3], ns[4]), (ns[4], ns[1]), # f1, f2, f3, f4 + ) +end + +# RefTetrahedron (refdim = 3): vertices for vertexdofs, edges for edgedofs, faces for facedofs and BC +function vertices(c::AbstractCell{3, RefTetrahedron}) + ns = get_node_ids(c) + return (ns[1], ns[2], ns[3], ns[4]) # v1, v2, v3, v4 +end +function edges(c::AbstractCell{3, RefTetrahedron}) + ns = get_node_ids(c) + return ( + (ns[1], ns[2]), (ns[2], ns[3]), (ns[3], ns[1]), # e1, e2, e3 + (ns[1], ns[4]), (ns[2], ns[4]), (ns[3], ns[4]), # e4, e5, e6 + ) +end +function faces(c::AbstractCell{3, RefTetrahedron}) + ns = get_node_ids(c) + return ( + (ns[1], ns[3], ns[2]), (ns[1], ns[2], ns[4]), # f1, f2 + (ns[2], ns[3], ns[4]), (ns[1], ns[4], ns[3]), # f3, f4 + ) +end + +# RefHexahedron (refdim = 3): vertices for vertexdofs, edges for edgedofs, faces for facedofs and BC +function vertices(c::AbstractCell{3, RefHexahedron}) + ns = get_node_ids(c) + return ( + ns[1], ns[2], ns[3], ns[4], ns[5], ns[6], ns[7], ns[8], # v1, ..., v8 + ) +end +function edges(c::AbstractCell{3, RefHexahedron}) + ns = get_node_ids(c) + return ( + (ns[1], ns[2]), (ns[2], ns[3]), (ns[3], ns[4]), (ns[4], ns[1]), # e1, e2, e3, e4 + (ns[5], ns[6]), (ns[6], ns[7]), (ns[7], ns[8]), (ns[8], ns[5]), # e5, e6, e7, e8 + (ns[1], ns[5]), (ns[2], ns[6]), (ns[3], ns[7]), (ns[4], ns[8]), # e9, e10, e11, e12 + ) +end +function faces(c::AbstractCell{3, RefHexahedron}) + ns = get_node_ids(c) + return ( + (ns[1], ns[4], ns[3], ns[2]), (ns[1], ns[2], ns[6], ns[5]), # f1, f2 + (ns[2], ns[3], ns[7], ns[6]), (ns[3], ns[4], ns[8], ns[7]), # f3, f4 + (ns[1], ns[5], ns[8], ns[4]), (ns[5], ns[6], ns[7], ns[8]), # f5, f6 + ) +end + +# RefPrism (refdim = 3): vertices for vertexdofs, edges for edgedofs, faces for facedofs and BC +function vertices(c::AbstractCell{3, RefPrism}) + ns = get_node_ids(c) + return (ns[1], ns[2], ns[3], ns[4], ns[5], ns[6]) # v1, ..., v6 +end +function edges(c::AbstractCell{3, RefPrism}) + ns = get_node_ids(c) + return ( + (ns[2], ns[1]), (ns[1], ns[3]), (ns[1], ns[4]), (ns[3], ns[2]), # e1, e2, e3, e4 + (ns[2], ns[5]), (ns[3], ns[6]), (ns[4], ns[5]), (ns[4], ns[6]), # e5, e6, e7, e8 + (ns[6], ns[5]), # e9 + ) +end +function faces(c::AbstractCell{3, RefPrism}) + ns = get_node_ids(c) + return ( + (ns[1], ns[3], ns[2]), (ns[1], ns[2], ns[5], ns[4]), # f1, f2 + (ns[3], ns[1], ns[4], ns[6]), (ns[2], ns[3], ns[6], ns[5]), # f3, f4 + (ns[4], ns[5], ns[6]), # f5 + ) +end + + +###################################################### +# Concrete implementations of AbstractCell interface # +###################################################### + +# Lagrange interpolation based cells +struct Line <: AbstractCell{1, RefInterval} nodes::NTuple{ 2, Int} end +struct QuadraticLine <: AbstractCell{1, RefInterval} nodes::NTuple{ 3, Int} end +struct Triangle <: AbstractCell{2, RefTriangle} nodes::NTuple{ 3, Int} end +struct QuadraticTriangle <: AbstractCell{2, RefTriangle} nodes::NTuple{ 6, Int} end +struct Quadrilateral <: AbstractCell{2, RefQuadrilateral} nodes::NTuple{ 4, Int} end +struct QuadraticQuadrilateral <: AbstractCell{2, RefQuadrilateral} nodes::NTuple{ 9, Int} end +struct Tetrahedron <: AbstractCell{3, RefTetrahedron} nodes::NTuple{ 4, Int} end +struct QuadraticTetrahedron <: AbstractCell{3, RefTetrahedron} nodes::NTuple{10, Int} end +struct Hexahedron <: AbstractCell{3, RefHexahedron} nodes::NTuple{ 8, Int} end +struct QuadraticHexahedron <: AbstractCell{3, RefHexahedron} nodes::NTuple{27, Int} end +struct Wedge <: AbstractCell{3, RefPrism} nodes::NTuple{ 6, Int} end + +default_interpolation(::Type{Line}) = Lagrange{1, RefCube, 1}() +default_interpolation(::Type{QuadraticLine}) = Lagrange{1, RefCube, 2}() +default_interpolation(::Type{Triangle}) = Lagrange{2, RefTetrahedron, 1}() +default_interpolation(::Type{QuadraticTriangle}) = Lagrange{2, RefTetrahedron, 2}() +default_interpolation(::Type{Quadrilateral}) = Lagrange{2, RefCube, 1}() +default_interpolation(::Type{QuadraticQuadrilateral}) = Lagrange{2, RefCube, 2}() +default_interpolation(::Type{Tetrahedron}) = Lagrange{3, RefTetrahedron, 1}() +default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3, RefTetrahedron, 2}() +default_interpolation(::Type{Hexahedron}) = Lagrange{3, RefCube, 1}() +default_interpolation(::Type{QuadraticHexahedron}) = Lagrange{3, RefCube, 2}() +default_interpolation(::Type{Wedge}) = Lagrange{3, RefPrism, 1}() + +# TODO: Deprecate +const Line2D = Line +const Line3D = Line +const Quadrilateral3D = Quadrilateral +edges(c::Quadrilateral3D) = faces(c) + +# Serendipity interpolation based cells +# TODO: Naming? +struct SQuadraticQuadrilateral <: AbstractCell{2, RefQuadrilateral} nodes::NTuple{ 8, Int} end +struct SQuadraticHexahedron <: AbstractCell{3, RefHexahedron} nodes::NTuple{20, Int} end + +default_interpolation(::Type{SQuadraticQuadrilateral}) = Serendipity{2, RefCube, 2}() +default_interpolation(::Type{SQuadraticHexahedron}) = Serendipity{3, RefCube, 2}() + + +############ +# Topology # +############ +# TODO: Move topology stuff to src/Grid/Topology.jl or something struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}} neighbor_info::Vector{T} @@ -135,21 +263,21 @@ function Base.show(io::IO, ::MIME"text/plain", n::EntityNeighborhood) end """ - face_npoints(::AbstractCell{dim,N,M) + face_npoints(::AbstractCell) Specifies for each subtype of AbstractCell how many nodes form a face """ -face_npoints(::Cell{2,N,M}) where {N,M} = 2 -face_npoints(::Cell{3,4,1}) = 4 #not sure how to handle embedded cells e.g. Quadrilateral3D +face_npoints(::AbstractCell{2}) = 2 +# face_npoints(::Cell{3,4,1}) = 4 #not sure how to handle embedded cells e.g. Quadrilateral3D """ - edge_npoints(::AbstractCell{dim,N,M) + edge_npoints(::AbstractCell) Specifies for each subtype of AbstractCell how many nodes form an edge """ -edge_npoints(::Cell{3,4,1}) = 2 #not sure how to handle embedded cells e.g. Quadrilateral3D -face_npoints(::Cell{3,N,6}) where N = 4 -face_npoints(::Cell{3,N,4}) where N = 3 -edge_npoints(::Cell{3,N,M}) where {N,M} = 2 +# edge_npoints(::Cell{3,4,1}) = 2 #not sure how to handle embedded cells e.g. Quadrilateral3D +face_npoints(::AbstractCell{3, RefHexahedron}) = 4 +face_npoints(::AbstractCell{3, RefTetrahedron}) = 3 +edge_npoints(::AbstractCell{3}) = 2 -getdim(::Cell{dim}) where dim = dim +getdim(::AbstractCell{refdim}) where {refdim} = refdim abstract type AbstractTopology end @@ -810,46 +938,6 @@ function Base.show(io::IO, ::MIME"text/plain", grid::Grid) print(io, " cells and $(getnnodes(grid)) nodes") end -# Functions to uniquely identify vertices, edges and faces, used when distributing -# dofs over a mesh. For this we can ignore the nodes on edged, faces and inside cells, -# we only need to use the nodes that are vertices. -# 1D: vertices -faces(c::Union{Line,QuadraticLine}) = (c.nodes[1], c.nodes[2]) -vertices(c::Union{Line,Line2D,Line3D,QuadraticLine}) = (c.nodes[1], c.nodes[2]) -# 2D: vertices, faces -faces(c::Line2D) = ((c.nodes[1],c.nodes[2]),) -vertices(c::Union{Triangle,QuadraticTriangle}) = (c.nodes[1], c.nodes[2], c.nodes[3]) -faces(c::Union{Triangle,QuadraticTriangle}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1])) -vertices(c::Union{Quadrilateral,Quadrilateral3D,QuadraticQuadrilateral}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4]) -faces(c::Union{Quadrilateral,QuadraticQuadrilateral}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1])) -# 3D: vertices, edges, faces -edges(c::Line3D) = ((c.nodes[1],c.nodes[2]),) -vertices(c::Union{Tetrahedron,QuadraticTetrahedron}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4]) -edges(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]), (c.nodes[1],c.nodes[4]), (c.nodes[2],c.nodes[4]), (c.nodes[3],c.nodes[4])) -faces(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[4]), (c.nodes[2],c.nodes[3],c.nodes[4]), (c.nodes[1],c.nodes[4],c.nodes[3])) -vertices(c::Union{Hexahedron,Cell{3,20,6}}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6], c.nodes[7], c.nodes[8]) -edges(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]), (c.nodes[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8])) -faces(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[6],c.nodes[5]), (c.nodes[2],c.nodes[3],c.nodes[7],c.nodes[6]), (c.nodes[3],c.nodes[4],c.nodes[8],c.nodes[7]), (c.nodes[1],c.nodes[5],c.nodes[8],c.nodes[4]), (c.nodes[5],c.nodes[6],c.nodes[7],c.nodes[8])) -edges(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1])) -faces(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2],c.nodes[3],c.nodes[4]),) - -vertices(c::Wedge) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6]) -edges(c::Wedge) = ((c.nodes[2],c.nodes[1]), (c.nodes[1],c.nodes[3]), (c.nodes[1],c.nodes[4]), (c.nodes[3],c.nodes[2]), (c.nodes[2],c.nodes[5]), (c.nodes[3],c.nodes[6]), (c.nodes[4],c.nodes[5]), (c.nodes[4],c.nodes[6]), (c.nodes[6],c.nodes[5])) -faces(c::Wedge) = ((c.nodes[1],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[5],c.nodes[4]), (c.nodes[3],c.nodes[1],c.nodes[4],c.nodes[6]), (c.nodes[2],c.nodes[3],c.nodes[6],c.nodes[5]), (c.nodes[4],c.nodes[5],c.nodes[6])) - -# random stuff -default_interpolation(::Union{Type{Line},Type{Line2D},Type{Line3D}}) = Lagrange{1,RefCube,1}() -default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}() -default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}() -default_interpolation(::Type{QuadraticTriangle}) = Lagrange{2,RefTetrahedron,2}() -default_interpolation(::Union{Type{Quadrilateral},Type{Quadrilateral3D}}) = Lagrange{2,RefCube,1}() -default_interpolation(::Type{QuadraticQuadrilateral}) = Lagrange{2,RefCube,2}() -default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}() -default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3,RefTetrahedron,2}() -default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}() -default_interpolation(::Type{Cell{3,20,6}}) = Serendipity{3,RefCube,2}() -default_interpolation(::Type{Wedge}) = Lagrange{3,RefPrism,1}() - """ boundaryfunction(::Type{<:BoundaryIndex}) diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl index d5ca1c8f36..5ea475d214 100644 --- a/src/Grid/grid_generators.jl +++ b/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 diff --git a/src/deprecations.jl b/src/deprecations.jl index 202fe6faa7..eb93162076 100644 --- a/src/deprecations.jl +++ b/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 diff --git a/src/exports.jl b/src/exports.jl index 8933912339..2e28d9755a 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -52,10 +52,12 @@ export Quadrilateral, Quadrilateral3D, QuadraticQuadrilateral, + SQuadraticQuadrilateral, Tetrahedron, QuadraticTetrahedron, Hexahedron, - #QuadraticHexahedron, + QuadraticHexahedron, + SQuadraticHexahedron, Wedge, CellIndex, FaceIndex, diff --git a/src/iterators.jl b/src/iterators.jl index fefe12c40e..ed8470ecbc 100644 --- a/src/iterators.jl +++ b/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[]] ################## diff --git a/test/test_apply_analytical.jl b/test/test_apply_analytical.jl index 5e58242585..d1f8658078 100644 --- a/test/test_apply_analytical.jl +++ b/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) diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl index 67c7e12d1c..a100e86e2f 100644 --- a/test/test_grid_dofhandler_vtk.jl +++ b/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)