Interpolation parameterization

[fredrikekre]

Currently Ferrite only implements scalar interpolations, i.e. interpolations where the shape functions are scalars. When using these for vector-valued function approximations we automatically "vectorize" these in two places:

1. Implicitly when passing scalar interpolations to CellVectorValues:

   Ferrite.jl/src/FEValues/cell_values.jl

   Line 103 in 247d7c3

   function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)

2. Explicitly when adding fields to the DofHandler as the dim argument:

   Ferrite.jl/src/Dofs/DofHandler.jl

   Line 217 in 247d7c3

   function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation)

This becomes a bit awkward for vector-valued interpolations such as e.g. Nédélec or Raviart–Thomas (xref #569). Since these shouldn't be vectorized it isn't clear if one should use CellScalarValues or CellVectorValues, and it isn't clear what to pass as the dim parameter to the DofHandler.

For these reasons I suggest adding abstract types ScalarInterpolation and VectorInterpolation (where all current Ferrite interpolations would be <: ScalarInterpolation):

abstract type Interpolation{dim,shape,order} end
abstract type ScalarInterpolation{dim,shape,order} <: Interpolation{dim,shape,order} end
abstract type VectorInterpolation{dim,shape,order} <: Interpolation{dim,shape,order} end

struct Lagrange{dim,shape,order} <: ScalarInterpolation{dim,shape,order} end
struct Nedelec{dim,refshape,order} <: VectorInterpolation{dim,refshape,order} end

To use a scalar interpolation like Lagrange we then also add a "vectorization layer" so that this can be encoded already in the interpolation, for example something like:

struct VectorizedInterpolation{dim,shape,order,base} <: VectorInterpolation{dim,shape,order}
    base::ScalarInterpolation{dim,shape,order}
    # The `ncomponents` parameter is currently redundant since we currently only can
    # vectorize such that `ncomponents === dim`.
    ncomponents::Int
end

(For some prior art, DefElement actually considers scalar Lagrange and vector Lagrange to be independent interpolations, and in deal.ii one also have to explicitly vectorize when creating an FESystem: https://www.dealii.org/current/doxygen/deal.II/classFESystem.html.)

Potentially we can implement ncomponents::Int * ip::ScalarInterpolation), or ip::ScalarInterpolation ^ ncomponents::Int) as alternative VectorizedInterpolation constructors, similar to what deal.ii does in the FESystem constructor.

With the type hierarchy above we can implement dispatches for CellValues directly:

CellValues(ip::ScalarInterpolation) = CellScalarValues(ip)
CellValues(ip::VectorInterpolation) = CellVectorValues(ip)

and we can (probably) remove the FieldTrait:

Ferrite.jl/src/FEValues/common_values.jl

Lines 5 to 7 in aa26130

	abstract type FieldTrait end
	struct VectorValued <: FieldTrait end
	struct ScalarValued <: FieldTrait end

in favor of dispatching on CellValues{<:ScalarInterpolation} and CellValues{<:VectorInterpolation}.

For the Taylor-Hood element this would then look like:

# Interpolations
interpolation_u = VectorizedInterpolation(Lagrange{dim,RefTetra,2}()) # or Lagrange{dim,RefTetra,2}() ^ dim
interpolation_p = Lagrange{dim,RefTetra,1}()

# CellValues
qr = QuadratureRule(...)
cellvalues_u = CellValues(qr, cellvalues_u)
cellvalues_p = CellValues(qr, cellvalues_p)

# DofHandler
dh = DofHandler(grid)
add!(dh, :u, interpolation_u) # no need to specify dim here
add!(dh, :p, interpolation_p)
close!(dh)

[KnutAM]

I like this direction!

One question/thought: Wouldn't ncomponents need to be a type parameter (if ncomponents != dim is implemented this would lead to MixedTensor I think)
What's the advantage of wrapping the interpolation instead of calling it e.g. VectorLagrange directly?

[fredrikekre]

Yea, I suppose it could be, but I think it could also be handled in the CellValues constructor perhaps. The suggested VectorInterpolations would have the same problem. I think it could be nice to keep the interpolations independent of this though, it is only when you couple it with a spacedim from a grid that you have that problem, I think.

What's the advantage of wrapping the interpolation instead of calling it e.g. VectorLagrange directly?

Thats one option, but doesn't generalize as nicely IMO.

[fredrikekre]

For example, I don't think you could include the MixedTensor stuff in the current CellValue struct layout anyway, so that would have to be separate.

[fredrikekre]

We could also consider removing the order parameter for the abstract types, and only have it for the specific implementations where it make sense (e.g. Lagrange).

[fredrikekre]

And also related to this issue, perhaps it would be nice to introduce RefLine, RefQuadrilateral, RefTriangle instead of having RefCube, dim=1, RefCube, dim=2, and RefTetrahedron, dim=2, respectively.

interpolation_u  = Lagrange{RefTriangle,2}()
interpolation_p  = Lagrange{RefTriangle,1}()

looks better than

interpolation_u  = Lagrange{2,RefTetrahedron,2}()
interpolation_p  = Lagrange{2,RefTetrahedron,1}()

IMO.

This might make it more difficult to write a dimension independent code, but... is that of much use in practice? We already have different cells (i.e. Triangle isn't Tetrahedron{2}) so you already have to branch on the dimension in a couple of places.

[KnutAM]

Thats one option, but doesn't generalize as nicely IMO.

Agreed (that's why I crossed it out:))

Yea, I suppose it could be, but I think it could also be handled in the CellValues constructor perhaps.

True, it would make it type-unstable, I think, but that shouldn't matter since it is during setup.

I don't think you could include the MixedTensor stuff in the current CellValue struct layout anyway

With the intended implementation, having MixedTensors in that would infer zero additional costs even for cases with regular Tensors. Would just need the added double dimensions. But since it isn't certain that will be implemented, perhaps it isn't relevant anyways it this stage.

introduce RefLine, RefQuadrilateral, RefTriangle

For my applications, that would not be an issue. It would also be nice for the quadrature. Furthermore, it clarifies that it is the dim of the shape and not the embedded dimension we are talking about (I think I confused that before).
If needed later, I suppose one could implement the traits RefBox and RefSimplex.

[termi-official]

Thanks for summarizing the discussion Fredrik! I already started implementing our ideas, but hit a hard wall when trying to fix up the dof distribution for high order elements. This is already an issue in the current api, when we have more than 1 distinct dof locations on the interior of an edge/face with more than one entry. I.e.

Ferrite.jl/src/interpolations.jl

Line 242 in a917a24

facedof_interior_indices(::Interpolation)

and

Ferrite.jl/src/interpolations.jl

Line 223 in a917a24

edgedof_interior_indices(::Interpolation)

are ambiguous. Another problem here is that we also need to separate nodal (e.g. Lagrange) from non-nodal (e.g. Nedelec) interpolations. This raises the question what might better better here, a trait system or subtyping or a combination of both.

Before this is lost, we likely can make the change non-breaking (for the user-facing api) by adding the vectorized versions to the dof handler with the corresponding add! call and also creating the vectorized interpolation in the CellVectorValues constructor.

[fredrikekre]

True, it would make it type-unstable, I think, but that shouldn't matter since it is during setup.

I was imagining something like EmbeddedCellValues{spacedim}(ip::Interpolation{dim}) which should not be type-unstable.

With the intended implementation, having MixedTensors in that would infer zero additional costs even for cases with regular Tensors. Would just need the added double dimensions. But since it isn't certain that will be implemented, perhaps it isn't relevant anyways it this stage.

I see, would you cast it to the standard tensor types in shape_gradient etc for the common case of dim == spacedim? Perhaps it can at least be prototyped as a separate CellValue struct regardless.

If needed later, I suppose one could implement the traits RefBox and RefSimplex.

There is already RefPrism . (Incidentally, RefPrism dim = 2 might be more appropriate for the triangle instead of RefTetrahedron dim = 2 since the former is just an extruded triangle.

Another problem here is that we also need to separate nodal (e.g. Lagrange) from non-nodal (e.g. Nedelec) interpolations.

Why do you need to do that? I assumed for Nedelec we could just use facedof_interior_indices like this:

Ferrite.jl/src/interpolations.jl

Lines 1125 to 1140 in d5f18d9

	###################################################################
	# Nédélec dim 2 RefTetrahedron order 1 #
	# https://defelement.com/elements/examples/triangle-N1curl-1.html #
	###################################################################
	getnbasefunctions(::Nedelec{2,RefTetrahedron,1}) = 3
	facedof_interior_indices(::Nedelec{2,RefTetrahedron,1}) = ((1,), (2,), (3,))
	function value(ip::Nedelec{2,RefTetrahedron,1}, i::Int, ξ::Vec{2,T}) where T
	x, y = ξ
	i == 1 && return Vec{2,T}(( - y, x ))
	i == 2 && return Vec{2,T}(( y, 1 - x ))
	i == 3 && return Vec{2,T}(( 1 - y, x ))
	throw(ArgumentError("no shape function $i for interpolation $ip"))
	end

[termi-official]

True, it would make it type-unstable, I think, but that shouldn't matter since it is during setup.

I was imagining something like EmbeddedCellValues{spacedim}(ip::Interpolation{dim}) which should not be type-unstable.

With the intended implementation, having MixedTensors in that would infer zero additional costs even for cases with regular Tensors. Would just need the added double dimensions. But since it isn't certain that will be implemented, perhaps it isn't relevant anyways it this stage.

I see, would you cast it to the standard tensor types in shape_gradient etc for the common case of dim == spacedim? Perhaps it can at least be prototyped as a separate CellValue struct regardless.

Note that users also sometimes vectorize variables which might note be inherently vectors (see e.g. the angle in the shell example) leading to dim != spacedim so EmbeddedCellValues might not be suitable as a name.

Another problem here is that we also need to separate nodal (e.g. Lagrange) from non-nodal (e.g. Nedelec) interpolations.

Why do you need to do that? I assumed for Nedelec we could just use facedof_interior_indices like this:

Ferrite.jl/src/interpolations.jl

Lines 1125 to 1140 in d5f18d9

	###################################################################
	# Nédélec dim 2 RefTetrahedron order 1 #
	# https://defelement.com/elements/examples/triangle-N1curl-1.html #
	###################################################################
	getnbasefunctions(::Nedelec{2,RefTetrahedron,1}) = 3
	facedof_interior_indices(::Nedelec{2,RefTetrahedron,1}) = ((1,), (2,), (3,))
	function value(ip::Nedelec{2,RefTetrahedron,1}, i::Int, ξ::Vec{2,T}) where T
	x, y = ξ
	i == 1 && return Vec{2,T}(( - y, x ))
	i == 2 && return Vec{2,T}(( y, 1 - x ))
	i == 3 && return Vec{2,T}(( 1 - y, x ))
	throw(ArgumentError("no shape function $i for interpolation $ip"))
	end

Linear elements are no problems. For high order stuff I am not sure how it should be handled and also integral elements (like Nedelec or Raviart-Thomas) do not have reference coordinates, because full geometric entities are the reference thing (because they are defined as integral moments over subentities).

[KnutAM]

I see, would you cast it to the standard tensor types in shape_gradient etc for the common case of dim == spacedim? Perhaps it can at least be prototyped as a separate CellValue struct regardless.

Yes, (and if not it would convert on the first tensor operation that it is used in.)

But yes, it should be prototyped separately and I think there is quite some work to get the same performance as regular tensors for the mixed cases.

[fredrikekre]

For high order stuff I am not sure how it should be handled and also integral elements (like Nedelec or Raviart-Thomas) do not have reference coordinates, because full geometric entities are the reference thing (because they are defined as integral moments over subentities).

Yea but reference_coordinates is used in a single place only

Ferrite.jl/src/FEValues/face_values.jl

Line 245 in aa26130

interpolation_coords = reference_coordinates(func_interpol)

(for a very long time it was implemented in test/ for testing purposes only). For applying boundary conditions I think another approach would be needed anyway.

[fredrikekre]

Regarding #676 (comment), an alternative would be to add the reference dimension as a parameter instead:

struct RefCube{rdim} end

This, or the proposal to add RefQuadrilateral would remove the need to always keep two type parameters for the reference entity. For example, in Lagrange{2,RefCube,2}() it isn't so clear that the first 2 is related to the dimension, and the last 2 for the order. Lagrange{RefCube{2},2}() or Lagrange{RefQuadrilateral,2}() are both better options then, IMO.

[fredrikekre]

For e.g. CellValues we need rdim as it's own parameter for the tensor storage anyway so maybe it has to be kept. We could still add

julia> struct RefTriangle <: Ferrite.AbstractRefShape end

julia> function Lagrange{RefTriangle,order}() where order
           return Lagrange{2,RefTetrahedron,order}()
       end

julia> Lagrange{RefTriangle,1}()
Lagrange{2, RefTetrahedron, 1}()

to simplify interpolation construction a bit, perhaps.

[KnutAM]

For e.g. CellValues we need rdim as it's own parameter for the tensor storage anyway so maybe it has to be kept.

Could that be solved by defining e.g. get_reference_dim(::Type{RefTriangle}) = Val(2) ?

[fredrikekre]

taylor_hood_ip = Lagrange{2,RefTriangle,2}()^2 × Lagrange{2,RefTriangle,1}()

[termi-official]

Although I really like such syntactic sugar to construct elements, what should value return in this construct? And how do we interface this with the dof distribution?

With such syntax we could also think about tensor-product elements like for example space_time_interpolation = Lagrange{3,RefCube,1}() ⊗ DiscontinuousLagrange{1,RefLine,1}() or ip_quad = Lagrange{1,RefLine,3}() ⊗ Lagrange{1,RefLine,3}().

[fredrikekre]

Just pointing out that we can make this work:

fields = (:u, :p) ∈ Lagrange{2,RefTriangle,2}()^2 × Lagrange{2,RefTriangle,1}()

qr = QuadratureRule(...)
cellvalues = CellValues(qr, fields)

dh = DofHandler(grid)
add!(dh, fields)
close!(dh)

[fredrikekre]

what should value return in this construct?

I imagine this would return a "multivalues" (#674)

[termi-official]

What about something like

fields = (:u,) ∈ Lagrange{2,RefTriangle,2}()^2 ⊕ Lagrange{2,RefCube,2}()^2

for mixed grids? Or maybe even

fields = (:u,) ∈ Lagrange{rdim=2,order=2}()^2

where Lagrange{2,2}() returns some function which maps refshape -> Lagrange{2,ref_shape,2}()
to simplify things further.

Edit: Analogue treatment for the quadrature rules.

Edit 2: MixedMultiValues?.

[termi-official]

With #679 we can eliminate dim, because it currently creates some weird inconsistencies.

[fredrikekre]

Or maybe even

fields = (:u,) ∈ Lagrange{rdim=2,order=2}()^2

where Lagrange{2,2}() returns some function which maps refshape -> Lagrange{2,ref_shape,2}()

One of the main motivations for the subregion stuff is that it also allows you to specify different fields on different domains, so the extra explicitness for the case when you have the same field everywhere, but different shapes, doesn't bother me that much.

One alternative could be to have RefUnknown that simply does runtime case-switching in a DynamicCellValues implementation.

[termi-official]

Sure. It also does not bother me to be explicit here, but I thought that it might be a nice convenience feature.

[fredrikekre]

Resolved by #694, mostly.