Closed fredrikekre closed 1 year ago
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?~
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.
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.
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
).
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.
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
.
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. https://github.com/Ferrite-FEM/Ferrite.jl/blob/a917a24ca0437d2d4b30777c0f7c4d637d13ade0/src/interpolations.jl#L242 and https://github.com/Ferrite-FEM/Ferrite.jl/blob/a917a24ca0437d2d4b30777c0f7c4d637d13ade0/src/interpolations.jl#L223 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.
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: https://github.com/Ferrite-FEM/Ferrite.jl/blob/d5f18d9e838e91122a95bd9442ad91a433494131/src/interpolations.jl#L1125-L1140
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 ofdim == 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:
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).
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.
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 https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/FEValues/face_values.jl#L245 (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.
Regarding https://github.com/Ferrite-FEM/Ferrite.jl/issues/676#issuecomment-1503114298, 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.
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.
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)
?
taylor_hood_ip = Lagrange{2,RefTriangle,2}()^2 × Lagrange{2,RefTriangle,1}()
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}()
.
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)
what should value return in this construct?
I imagine this would return a "multivalues" (#674)
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
?.
With #679 we can eliminate dim
, because it currently creates some weird inconsistencies.
Or maybe even
fields = (:u,) ∈ Lagrange{rdim=2,order=2}()^2
where
Lagrange{2,2}()
returns some function which mapsrefshape -> 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.
Sure. It also does not bother me to be explicit here, but I thought that it might be a nice convenience feature.
Resolved by #694, mostly.
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:
CellVectorValues
: https://github.com/Ferrite-FEM/Ferrite.jl/blob/247d7c3efdd6325a78c7274180d4a439921dcc78/src/FEValues/cell_values.jl#L103dim
argument: https://github.com/Ferrite-FEM/Ferrite.jl/blob/247d7c3efdd6325a78c7274180d4a439921dcc78/src/Dofs/DofHandler.jl#L217This becomes a bit awkward for vector-valued interpolations such as e.g. Nédélec or Raviart–Thomas (xref https://github.com/Ferrite-FEM/Ferrite.jl/issues/569). Since these shouldn't be vectorized it isn't clear if one should use
CellScalarValues
orCellVectorValues
, and it isn't clear what to pass as thedim
parameter to the DofHandler.For these reasons I suggest adding abstract types
ScalarInterpolation
andVectorInterpolation
(where all current Ferrite interpolations would be<: ScalarInterpolation
):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:
(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)
, orip::ScalarInterpolation ^ ncomponents::Int)
as alternativeVectorizedInterpolation
constructors, similar to what deal.ii does in the FESystem constructor.With the type hierarchy above we can implement dispatches for
CellValues
directly:and we can (probably) remove the FieldTrait: https://github.com/Ferrite-FEM/Ferrite.jl/blob/aa261308181a5d856742385eafaa88e87481ad16/src/FEValues/common_values.jl#L5-L7 in favor of dispatching on
CellValues{<:ScalarInterpolation}
andCellValues{<:VectorInterpolation}
.For the Taylor-Hood element this would then look like: