GitHub

Interpolations

All Interpolations should subtype Interpolation{shape, order}, where shape <: AbstractRefShape is the reference shape for which the interpolation is defined and order is the characteristic interpolation order. The how-to at bottom of this page describes how to implement a new interpolation.

Methods to be implemented for a new interpolation

Always

Ferrite.vertexdof_indices — Method
vertexdof_indices(ip::Interpolation)

A tuple containing tuples of local dof indices for the respective vertex in local enumeration on a cell defined by vertices(::Cell). The vertex enumeration must match the vertex enumeration of the corresponding geometrical cell.

Note

The dofs appearing in the tuple must be continuous and increasing! The first dof must be the 1, as vertex dofs are enumerated first.

source
Ferrite.facedof_indices — Method
facedof_indices(ip::Interpolation)

A tuple containing tuples of all local dof indices for the respective face in local enumeration on a cell defined by faces(::Cell). The face enumeration must match the face enumeration of the corresponding geometrical cell.

source
Ferrite.facedof_interior_indices — Method
facedof_interior_indices(ip::Interpolation)

A tuple containing tuples of the local dof indices on the interior of the respective face in local enumeration on a cell defined by faces(::Cell). The face enumeration must match the face enumeration of the corresponding geometrical cell. Note that the vertex and edge dofs are included here.

Note

The dofs appearing in the tuple must be continuous and increasing! The first dof must be the computed via "last edge interior dof index + 1", if face dofs exist.

source
Ferrite.edgedof_indices — Method
edgedof_indices(ip::Interpolation)

A tuple containing tuples of local dof indices for the respective edge in local enumeration on a cell defined by edges(::Cell). The edge enumeration must match the edge enumeration of the corresponding geometrical cell.

The dofs are guaranteed to be aligned with the local ordering of the entities on the oriented edge. Here the first entries are the vertex dofs, followed by the edge interior dofs.

source
Ferrite.edgedof_interior_indices — Method
edgedof_interior_indices(ip::Interpolation)

A tuple containing tuples of the local dof indices on the interior of the respective edge in local enumeration on a cell defined by edges(::Cell). The edge enumeration must match the edge enumeration of the corresponding geometrical cell. Note that the vertex dofs are included here.

Note

The dofs appearing in the tuple must be continuous and increasing! The first dof must be computed via "last vertex dof index + 1", if edge dofs exist.

source
Ferrite.volumedof_interior_indices — Method
volumedof_interior_indices(ip::Interpolation)

Tuple containing the dof indices associated with the interior of a volume.

Note

The dofs appearing in the tuple must be continuous and increasing, volumedofs are enumerated last.

source
Ferrite.adjust_dofs_during_distribution — Method
adjust_dofs_during_distribution(::Interpolation)

This function must return true if the dofs should be adjusted (i.e. permuted) during dof distribution. This is in contrast to i) adjusting the dofs during reinit! in the assembly loop, or ii) not adjusting at all (which is not needed for low order interpolations, generally).

source

For special interpolations

Discontinuous interpolations

For discontinuous interpolations, implementing the following methods might be required to apply Dirichlet boundary conditions.

Ferrite.is_discontinuous — Method
is_discontinuous(::Interpolation)
is_discontinuous(::Type{<:Interpolation})

Checks whether the interpolation is discontinuous (i.e. DiscontinuousLagrange)

source
Ferrite.dirichlet_vertexdof_indices — Method
dirichlet_vertexdof_indices(ip::Interpolation)

A tuple containing tuples of local dof indices for the respective vertex in local enumeration on a cell defined by vertices(::Cell). The vertex enumeration must match the vertex enumeration of the corresponding geometrical cell. Used internally in ConstraintHandler and defaults to vertexdof_indices(ip::Interpolation) for continuous interpolation.

Note

The dofs appearing in the tuple must be continuous and increasing! The first dof must be the 1, as vertex dofs are enumerated first.

source
Ferrite.dirichlet_edgedof_indices — Method
dirichlet_edgedof_indices(ip::Interpolation)

A tuple containing tuples of local dof indices for the respective edge in local enumeration on a cell defined by edges(::Cell). The edge enumeration must match the edge enumeration of the corresponding geometrical cell. Used internally in ConstraintHandler and defaults to edgedof_indices(ip::Interpolation) for continuous interpolation.

The dofs are guaranteed to be aligned with the local ordering of the entities on the oriented edge. Here the first entries are the vertex dofs, followed by the edge interior dofs.

source

Non-identity mapping

For interpolations that have a non-identity mapping (see Mapping of finite elements), the mapping type must be specified.

Ferrite.mapping_type — Function
mapping_type(ip::Interpolation)

Get the type of mapping from the reference cell to the real cell for an interpolation ip. Subtypes of ScalarInterpolation and VectorizedInterpolation return IdentityMapping(), but other non-scalar interpolations may request different mapping types.

source
Ferrite.get_direction — Function
get_direction(interpolation::Interpolation, shape_nr::Int, cell::AbstractCell)

Return the direction, ±1, of the cell entity (e.g. facet or edge) associated with the interpolation's shape function nr. shape_nr. This is only required for interpolations with non-identity mappings, where the direction is required during the mapping of the shape values.

TODO: Move the following description to get_edge_direction and get_face_direction following #1162

The direction of entities are defined as following the node numbers of the entity's vertices, vnodes. For an edge, vnodes[2] > vnodes[1] implies positive direction.

For a face, we first find index, i, of the smallest value in vnodes. Considering circular indexing, then a positive face has vnodes[i-1] > vnodes[i+1].

source

The defaults should always work

The following functions are defined such that they should work for any interpolations that defines the required functions specified above.

Ferrite.getrefshape — Method
Ferrite.getrefshape(::Interpolation)::AbstractRefShape

Return the reference element shape of the interpolation.

source
Ferrite.reference_shape_gradient — Method
reference_shape_gradient(ip::Interpolation, ξ::Vec, i::Int)

Evaluate the gradient of the ith shape function of the interpolation ip in reference coordinate ξ.

source
Ferrite.boundarydof_indices — Function
boundarydof_indices(::Type{<:BoundaryIndex})

Helper function to generically dispatch on the correct dof sets of a boundary entity.

source
Ferrite.reference_shape_values! — Function
reference_shape_values!(values::AbstractArray{T}, ip::Interpolation, ξ::Vec)

Evaluate all shape functions of ip at once at the reference point ξ and store them in values.

source
Ferrite.reference_shape_gradients! — Function
reference_shape_gradients!(gradients::AbstractArray, ip::Interpolation, ξ::Vec)

Evaluate all shape function gradients of ip at once at the reference point ξ and store them in gradients.

source
Ferrite.reference_shape_gradients_and_values! — Function
reference_shape_gradients_and_values!(gradients::AbstractArray, values::AbstractArray, ip::Interpolation, ξ::Vec)

Evaluate all shape function gradients and values of ip at once at the reference point ξ and store them in values.

source
Ferrite.reference_shape_hessians_gradients_and_values! — Function
reference_shape_hessians_gradients_and_values!(hessians::AbstractVector, gradients::AbstractVector, values::AbstractVector, ip::Interpolation, ξ::Vec)

Evaluate all shape function hessians, gradients and values of ip at once at the reference point ξ and store them in hessians, gradients, and values.

source
Ferrite.shape_value_type — Method
shape_value_type(ip::Interpolation, ::Type{T}) where T<:Number

Return the type of shape_value(ip::Interpolation, ξ::Vec, ib::Int).

source

How to implement a new interpolation

Warning

The API for implementing a new interpolation is not fully stable yet.

As an example, we will implement a Quadratic Triangle Interpolation (QTI), corresponding to Lagrange{RefTriangle, 2}(). Since this interpolation gives scalar-valued shape functions, it is a subtype of ScalarInterpolation,

struct QTI <: ScalarInterpolation{RefTriangle, 2}
end

Before we jump in and define the shape functions associated with this interpolation, we first need to consider the RefTriangle entities,

Ferrite.RefTriangle — Type
RefTriangle <: AbstractRefShape{2}

Reference triangle, reference dimension 2.

----------------+--------------------
Vertex numbers: | Vertex coordinates:
    2           |
    | \         | v1: 𝛏 = (1.0, 0.0)
    |   \       | v2: 𝛏 = (0.0, 1.0)
ξ₂^ |     \     | v3: 𝛏 = (0.0, 0.0)
  | 3-------1   |
  +--> ξ₁       |
----------------+--------------------
Edge numbers:   | Edge identifiers:
    +           |
    | \         | e1: (v1, v2)
    2   1       | e2: (v2, v3)
    |     \     | e3: (v3, v1)
    +---3---+   |
----------------+--------------------
Face numbers:   | Face identifiers:
    +           |
    | \         |
    |   \       | f1: (v1, v2, v3)
    |  1  \     |
    +-------+   |
----------------+--------------------
source

For this particular interpolation, we have one degree of freedom associated with each vertex, and one degree of freedom associated with each edge. Following the Ferrite numbering rules, we start by enumerating the vertices first, followed by the edges. The numbering is based on the RefTriangle shown above, and the actual shape functions are taken from defelement.org.

function Ferrite.reference_shape_value(ip::QTI, ξ::Vec{2}, shape_number::Int)
    ξ₁ = ξ[1]
    ξ₂ = ξ[2]
    γ = 1 - ξ₁ - ξ₂ # Helper
    shape_number == 1 && return ξ₁ * (2ξ₁ - 1) # v1: 1 at ξ = (1, 0)
    shape_number == 2 && return ξ₂ * (2ξ₂ - 1) # v2: 1 at ξ = (0, 1)
    shape_number == 3 && return γ * (2γ - 1)   # v3: 1 at ξ = (0, 0)
    shape_number == 4 && return 4ξ₁ * ξ₂ # e1: 1 at ξ = (½, ½)
    shape_number == 5 && return 4ξ₂ * γ  # e2: 1 at ξ = (0, ½)
    shape_number == 6 && return 4ξ₁ * γ  # e3: 1 at ξ = (½, 0)
    throw(ArgumentError("no shape function $shape_number for interpolation $ip"))
end

Having defined the actual interpolation function, we must now provide the information in comments above to Ferrite. We start by providing the reference coordinate for each shape number,

function Ferrite.reference_coordinates(::QTI)
    return [
        Vec{2, Float64}((1.0, 0.0)), # v1
        Vec{2, Float64}((0.0, 1.0)), # v2
        Vec{2, Float64}((0.0, 0.0)), # v3
        Vec{2, Float64}((0.5, 0.5)), # center of e1
        Vec{2, Float64}((0.0, 0.5)), # center of e2
        Vec{2, Float64}((0.5, 0.0)), # center of e3
    ]
end

We move on to defining which dof indices belong to each vertex. As we have 3 vertices on the RefTriangle, this should be a tuple with length 3, and each element contains all dof indices for the particular vertex. In this case, we only have a single dof per vertex,

Ferrite.vertexdof_indices(::QTI) = ((1,), (2,), (3,))

Note that the dof index numbering must be in order, first all in the first vertex, then all in the next vertex, and so on. (If we had two indices per vertex, we would have ((1,2), (3,4), (5,6))). Next, we number the edge dofs, following the same convention, matching how we defined it above,

Ferrite.edgedof_interior_indices(::QTI) = ((4,), (5,), (6,))
Ferrite.edgedof_indices(::QTI) = ((1, 2, 4,), (2, 3, 5,), (3, 1, 6,))

But here we need two functions, one for the interior indices (those that have not yet been included in lower-dimensional entities (vertices in this case)), and one for all indices for dofs that belong to the edge.

For the triangle, we only have a single face. However, all the dofs that belong to the face, also belongs to either the vertices or edges, hence we have no "interior" face dofs. So we get,

Ferrite.facedof_interior_indices(::QTI) = ((),)
Ferrite.facedof_indices(::QTI) = ((1, 2, 3, 4, 5, 6),)

Finally, since this is a 2d element, we have no volumedofs, and thus

Ferrite.volumedof_interior_indices(::QTI) = ()

It is necessary to tell Ferrite the total number of base functions, e.g.,

Ferrite.getnbasefunctions(::QTI) = 6

For distributing the degrees of freedom, higher order interpolations require that we account for the ordering on their entity. For example, if we have two interior dofs associated with an edge, we must match them the edges of the for the cells that share the edge, to make sure we consider the same ordering. Since we only have a single interior dof per edge, we don't need to adjust these, hence,

Ferrite.adjust_dofs_during_distribution(::QTI) = false
Tip

The function test_interpolation_properties in test/test_interpolations.jl can be used when implementation to check that some basic properties are fullfilled.