Decoupling of Geometrical, Topological and Algebraical Interfaces

[termi-official]

This PR explores possibilities to simplify and generalize the code base. Legacy code is removed and ambiguities are resolved on the way. Note that the naming is not final and I can sweep unwanted renamings back to the current naming scheme. The renaming is also there because I want to be able to test different combinations of the old and new interface and its equivalence in tests.
Addresses #191 #167 #395 #394 and is a foundation for #272 and 4D elements. Will be the basis to simplify #486 . Will also be a foundation to merge MixedDofHandler into DofHandler.

TODOs

• Introducing data types for the pieces
• Improve topology interface
• Type stability
• Arbitrary order mesh generators
• 4D Poisson problem with tesseract elements
• Test battery with manufactured solutions
• Add missing todos to todo list.

[kimauth]

This looks like a big effort! I'm not sure I'm following all the points that you want to make here, this is what I can see:

1. Being more consistent in using the grid interface instead of assuming Grid inside the code base. That's a good point and could probably be merged right away if it was it's own PR (in my opinion).
2. A new grid data structure. Could you maybe expand a little bit on what problems you're tackling with that? E.g. from what I understand it doesn't work any better for mixed grids than the current Grid, does it?
3. A new DofHandler. Can you explain what exactly it improves compared to DofHandler? It says, it is slightly faster than MixedDofHandler for a single concrete element type in the mesh, but in that case we would currently use DofHandler.

Most likely I didn't capture everything in this PR, feel free to add to my list!
However, it looks like a combination of different ideas where some could be instant improvements to the current code base, while others are extremely braking and would need very careful consideration. Perhaps it would be good to split it in several PRs if possible?

One of the main reasons that we haven't merged the DofHandlers yet is that it turned out hard to do it without breaking basic functionality. I'm all ears for suggestions how to solve this problem! #445 was an attempt as unifying things but it turned out to be more breaking than it initially looked, so it is stalling at the moment.

[termi-official]

Thanks for taking the time to look at the proposal Kim! Appreciate it. :)

Yes, this PR is not intended to be merged as is, but I will consecutively cherry pick changes into separate PRs. I am mostly exploring ideas to generalize and decouple things that should be separate.

My work here aims at several things.

1. I want to decouple the topology of the element (i.e. being a quadrilateral or a line or smth) from its geoemtric realization (e.g. being a quadrilateral in 2d or 3d). This separation was incomplete up to now. I do not think I can fully address this, but will try to address this "just" for all common geometries. I have to explore this further, but it should allow us e.g. to simplify the constraint handler.
2. The new grid data structure is there for two reasons. First, I have an easiert time to detect where the grid is used directly and not in the interface, and introduce the interface at relevant points. Second, "cell" was always confusing for me, because Cell is usually defined as a 3d volumetric entity. Note that the change from "element" to "cell" is a breaking change, because the interface is different. Also, the current cell interface is ambiguous - the new one should always allow unique element definitions for the geometry due to the coupling to the interpolation. Further, having a guarantee that the interpolation for the geometry is the correct one (which, as far as I understand, is currently not guaranteed) is helpful down the pipeline, e.g. for visualization.
3. The NewDofHandler is intended to merge both the DofHandler and the MixedDofHandler. At the current state this one is not even close to being done yet. A first improvement is that the code will be shorter (after refactoring), because we can separate out the topological routine in the beginning (marked with "Phase 1"), which just determines the connectivity and does some partial materialization. Note that this interface still is not super general. It should almous allow assembly in 4D tho. I try to be minially breaking over the current DofHandler.

More things that I am trying to address is the introduction of 4D and space-time elements, routines to change the interpolation for the mesh geometry (which might be handy for e.g. higher order FEM in mechanics), better integration for sub- and superparametric elements and finally to simplify the dof management for distributed code. Also, having a good topology interface will be crucial for highly efficient space-time adaptive codes, which is still the main goal in the long run. Hope that explais it better.

Edit 1: So another big issue is that we do not strictly separate the algebraic interface, e.g. if I remember correctly, then we assume in several parts in the code that the solution approximation is nodal and in tight relation with the geometry of each element. This is a blocker for most H(div) and H(curl) approximations. See e.g. https://defelement.com/elements/examples/triangle-Nedelec-1.html for an example and e.g. https://defelement.com/ciarlet.html for a short introduction to the math behind FEM with cleaner definitions.

[KristofferC]

if I remember correctly, then we assume in several parts in the code that the solution approximation is nodal and in tight relation with the geometry of each element

I think that should only be in the VTK output where something along those lines is assumed.

[termi-official]

Thanks for the pointer Kristoffer. Indeed, VTK export is one of the portions where I saw this assumption.

So, if I see it correctly then the "nodesset" really is intended to be a constraint on the geometric nodes of the grid and not on the algebraic nodes of the solution approximation?

Edit: I will propose something to remove this constraint.

[kimauth]

Yes, constraints on the nodes are constraints only on the geometric nodes. The ConstraintHandler even gives you a warning if you put constraints on node sets of sub-/superparameteric approximations.

[termi-official]

Sorry to ask, but what is the use-case for this? Wouldn't it make more sense to just directly constrain the dofs here instead of the nodes?

[kimauth]

E.g. on the cohesive elements that I presented at FerriteCon, there is no face on two of the sides, but it still makes sense to apply Dirichlet boundary conditions there (not Neumann though, thus no face). I usually constrain them via constraining these nodes (which I can get from the meshing routine I use, opposed to the dofs).

I would also say it is fairly common to use the same approximation on the geometry and the function and to constrain node sets then.

[kimauth]

Edit: I will propose something to remove this constraint.

What do you suggest to remove?
The restriction that a constraint on nodes does not work well if the function approximation is of a higher order than the geometry approximation?
Or the possibility to add constraint on node sets?

[fredrikekre]

It is also quite common to have to prescribe some corner in order to avoid rigid body motion. Prescribing the node is the only current way of reaching for those dofs.

[termi-official]

It is also quite common to have to prescribe some corner in order to avoid rigid body motion. Prescribing the node is the only current way of reaching for those dofs.

If I see this correctly, then wouldn't this be a constraint on the vertex?

What do you suggest to remove?
The restriction that a constraint on nodes does not work well if the function approximation is of a higher order than the geometry approximation?

I want to remove the constraint in the VTK export, that the geometrical nodes have to match the algebraic nodes in some way. This should be rather easy, as, in my understanding, VTK only needs the solution at the geometric nodes, which can be interpolated from the algebraic nodes via some <: Interpolation.

I would also say it is fairly common to use the same approximation on the geometry and the function and to constrain node sets then.

If we only talk about displacement fields then yes, it is common to use the same ansatz for geometry and solution.

Edit: Thanks for elaborating!

[kimauth]

I want to remove the constraint in the VTK export, that the geometrical nodes have to match the algebraic nodes in some way. This should be rather easy, as, in my understanding, VTK only needs the solution at the geometric nodes, which can be interpolated from the algebraic nodes via some <: Interpolation.

In that sense I think we have a solution already in terms of being able to use L2Projector and PointEvalHandler for obtaining values in arbitrary positions - if I understand you correctly that could achieve what you would want?
There is of course no convenience layer in the vtk export at the moment that does the work for the user in this case.

[termi-official]

I mean, yes, we could do it this way, but it will not be super efficient. You already have all information on the cell which you want to export, so there is no need to project the solution globally and search the node. I almost have a general solution in FerriteVis.jl which I would just adapt.

[codecov-commenter]

Codecov Report

Base: 92.34% // Head: 78.85% // Decreases project coverage by -13.48% ⚠️

Coverage data is based on head (353b1c3) compared to base (f6c7e23).
Patch coverage: 12.79% of modified lines in pull request are covered.

❗ Current head 353b1c3 differs from pull request most recent head fee7fb0. Consider uploading reports for the commit fee7fb0 to get more accurate results

Additional details and impacted files

@@             Coverage Diff             @@
##           master     #519       +/-   ##
===========================================
- Coverage   92.34%   78.85%   -13.49%     
===========================================
  Files          22       27        +5     
  Lines        3762     4455      +693     
===========================================
+ Hits         3474     3513       +39     
- Misses        288      942      +654     

Impacted Files	Coverage Δ
src/Dofs/NewDofHandler.jl	0.00% <0.00%> (ø)
src/Mesh/elements.jl	0.00% <0.00%> (ø)
src/Mesh/mesh.jl	0.00% <0.00%> (ø)
src/Mesh/mesh_generators.jl	0.00% <0.00%> (ø)
src/Mesh/utils.jl	0.00% <0.00%> (ø)
src/utils.jl	100.00% <ø> (ø)
src/Export/VTK.jl	68.83% <28.57%> (-22.55%)	⬇️
src/iterators.jl	81.57% <50.00%> (-15.35%)	⬇️
src/Grid/grid.jl	86.70% <86.36%> (-0.24%)	⬇️
src/Dofs/ConstraintHandler.jl	94.84% <100.00%> (+0.11%)	⬆️
... and 7 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[termi-official]

So the first set of changes which I want to pick out of this branch is the breaking change of the element interface. Here I have different suggestions

1. Rename of Cell to Element. I do not have a strong opinion here. However, Cell is a bit confusing in my opinion, because it often is used for something 3d volumetric, but we use it for 3d,2d,1d. I would provide type bindings with warnings if we introduce this change.
2. Changing the type parametrization of Cell.
   • Here I would like to change the dim parameter referring to the spatial dimension into which the element is embedded to the dimension of the reference geometry.
   • The second change which I propose is to remove the number of faces as a type parameter because it makes things ambiguous. Instead, I propose to replace the number of faces with either the reference geometry, or directly including the utilized nodal interpolation as a type parameters. Technically speaking only the latter gives us truly unique geometries, because the former does technically allow to move nodes around. No hard opinion about which one to choose in the end.
3. Rename compute_vertex_values to compute_nodal_values, because this function computes nodal values, not values on the vertices. I can also provide a proper compute_vertex_values tho, but since the behavior of this function changes this is also breaking.
4. Splitting up the grid.jl file into an cell.jl/element.jl file and a grid.jl/mesh.jl file for cleaner separation of the interfaces.
5. Rename RefTetrahedron to RefSimplex with deprecation binding. This way we have consistent naming with RefCube (which is not called RefHexahedron, to elaborate the naming issue)

Edit: The related files with a sketch of the ideas can be found in the Mesh folder.

[KnutAM]

Cool!

1. I was also confused by cell at first, but I've converged into liking it as the geometric entity. In my head, an element now refers to the cell + defined interpolations etc. Even https://defelement.com/ciarlet.html refer to cells as well... So I don't really see a sufficiently strong benefit to make such a breaking change.
2. Removing the num faces parameter and replacing it with a reference geometry sounds good to me! (But I guess it needs some example showing the need). Changing the definition of dim sounds like it could break many codes as some might rely on that dim deep inside their code. What are the main advantages/disadvantages of doing this?
3. Probably makes sense to rename and have a depreciation warning on the old name, and introduce the new function some months later. (Could provide an alternative name in the meantime if required with a warning that it will be replaced in the future)
4. +1
5. I agree that it is a bit strange with RefTetrahedron, but I would claim it is equivalent to RefCube as: Tetrahedron= 3D Simplex, Cube=3D Hypercube? In that case, we should rename to RefSimplex and RefHypercube (or RefBox)

[termi-official]

Thanks for taking time Knut! Let me answer in order.

1. Yes, my solution is not very clean either yet, because I still carry around some geometry and not only the topology in the Element. Fair point, but we can also use the Element as one way to e.g. materialize faces and edges. E.g. a face for a Hexahedron (Element{3,RefCube,8}) could be a Quadrilateral (Element{2,RefCube,4}). Not a huge point, but can possibly avoid future confusion.
2. Benefit here is that we can start removing ambiguities in the codebase like e.g. this (note that there is more) and I want guarantees that I can always easily access the geometric interpolation - note that this is not guaranteed to be the correct one. Edit 1: Please also note that we have not properly decoupled the different kinds of dimensions yet. E.g. we assume that field dimensions > 1 add extra dofs to vectorial finite elements (see e.g. __close!).
3. What about interpolate_to_nodes and interpolate_to_vertices then? It is also a bit more precise than compute, because compute could also be interpreted as some kind of projection.
4. :)
5. Indeed, good point. I will rethink the naming.

[KnutAM]

2. Yes, the ability to get the geometric interpolation does make a lot of sense - I was just a bit worried about how breaking it would be and thus think it would be nice to have a few clear reasons (e.g. required for MPI:)).
I thought the field dimension == number of dofs per geometric dof location, when would it have a different meaning?
3. Looking at the functions, I would say evaluate_at_* would make the most sense

[termi-official]

Please note that none of these first changes will help with improving the distributed assembly. I also play a bit around with partial materialization in this branch, but nothing is really satisfactory yet.

I thought the field dimension == number of dofs per geometric dof location, when would it have a different meaning?

For all H(curl) and H(div) conforming elements. They are vector-valued (e.g. 3d vector in 3d space) but have only 1 dof per vector and not 3. See e.g. https://defelement.com/elements/examples/triangle-Nedelec-1.html . I remember that there have been some issues when trying to implement the Nedelec elements (https://github.com/Ferrite-FEM/Ferrite.jl/tree/mk/nedelec).

[termi-official]

Made a full 180° and found my example from back then, in which the current AbstractCell interface is quite problematic.

It was Hermite elements, which often cannot be distinguished from some Lagrangian counterparts with the current interface (ran into this problem because one of my students was tasked to implement some Hermite formulation). Two representative examples here:

1. Cell{1,2,2},Cell{2,2,1},Cell{3,2,0} all represent elements with linear Lagrange ansatz as well as with cubic Hermite ansatz
2. Cell{3,4,1} represents a quadrilateral with bilinear ansatz as well as a triangle with cubic Hermite ansatz.

We can argue now that interpolations different from Lagrangian must require custom cell types, which I think is fine. However, all embedded elements are currently not out of the box functional, because we are missing an implementation for the Jacobian on these. Which is now dependent on the interpretation of these elements - are they shells/beams/... or surface/edge/... elements? To elaborate, with shells I refer to elements which are common in continuum mechanics, which come with some artificial treatment in their thickness direction. With surface elements I refer to elements which behave as expected, for example reaction systems defined on the surface of a sphere in 3D.

Edit: To summarize, my argument why we should default for the latter behavior is that shells require special treatment depending on their formulation, while regular surface/edge/... elements do not care what the surrounding dimension is and they behave consistently.

Edit 2: I still do not have a really satisfactory design, because my proposal now couples geometry and topology... for which I am not sure if this is good or bad.

[termi-official]

Almost all changes proposed here have been redistributed to separate PRs and merged in improved versions after more discussion on Slack (compare #651 ,#679, #694, #708 -- I think these are all). The only open remaining proposed change is the introduction of a dimension independent dof distribution algorithm (see __close! in NewDofHandler.jl) which is very raw and hard to read at this point. We may or may not cherry pick this and make also an improved PR out of it if needed.