FEM IGA interface

[TTitscher]

IGA and FEM have a common interface regarding nodal values and (Derivative)ShapeFunctions. When it comes to mesh algorithms (class MeshGenerator, GetNodeAt(), ConvertToInterpolationType(), PrismsCreate()), they differ.

We currently have an ElementInterface that is an interface for nodal vales and shape functions. Thus, it is useless for mesh algorithms that depend on the internals of FEM or IGA or whatever may come in the future. Since the cell works on an ElementCollection, we can move the polymorphism there. So I propose:

class ElementFem; // no common base class and ...
class ElementIga; // ... very specific methods for mesh algorithms

class ElementCollection
{
    virtual NodalValues ExtractNodalValues() const = 0;
    virtual ShapeFunctions GetShapeFunctions(...) const = 0;
    // ...
};

class ElementCollectionFem : public ElementCollection
{ // owns some ElementFem
};
class ElementCollectionIga : public ElementCollection
{ // owns some ElementIga
};

class MeshFem
{ // owns a container of ElementCollectionFem
};
class MeshIga
{ // owns a container of ElementCollectionIga
};

This allows to work on algorithms specific to one interpolation/mesh and make progress without thinking about all other possible interpolations. If it turns out that some algorithms may work on multiple mesh types, we can think about mesh interfaces.

[joergfunger]

What exactly is the reason to introduce this separation. As far as I understand the issue, mesh generation is done before these interpolations are defined. Thus, we only need different mesh generation tools, but why do we need separate containers? This somehow corresponds to the idea of a mesh file #85 that is created during meshing (and is different for fem and Iga) but the final solution data structure could be identical.

[TTitscher]

The reason is mainly the different algorithms I mentioned above. Node selection, for example. We found out that IGA requires different algorithms than FEM. And apparently no one knows how to do it in IGA. So this separation allows implementing it for FEM (and ignore IGA for now) to make some progress. If I only work on ElementInterface, I have no access to required FEM methods like GetLocalCoords(..)

Another thing is PrismsCreate. Gmsh is not capable of generating a prism layer around aggregates in 3D. This is done with PrismCreate() and specific methods of an FEM interpolation.

Another reason is that different meshes require different methods in general. E.g. a Mesh::AddInterpolationType(interpolationQuadLinear...) is only valid for FEM meshes. And Mesh::RefineKnots() only makes sense for IGA.

[joergfunger]

All the examples you are mentioning relate to preprocessing steps, but the data structure is just for the calculation. Wouldn't it be better to have separate meshing (data) structures, but then perform the computation on a unified data structure?

[TTitscher]

Maybe there is a misunderstanding. IMO, you rephrased my idea. :) Perform meshing (+algorithms) on separate data structures (MeshFem, MeshIga, MeshTitscherSpecial, maybe MeshMpi) and perform the calculation on the common interface ElementCollection. So nothing changes in the cell. The cell still works with references of ElementCollection.

[joergfunger]

What I meant were different data structures for mesh generation (preprocessing) and then forget all those information and perform the calculation on a unified structure (that is now very different from the mesh data structure that might include much more geometry and neighborhood information than is actually required for the computation part

[TTitscher]

Can you maybe provide a brief outline of how many classes of what types you mean? I still have troubles understanding the problem/your suggestion.

I propose those mesh classes. MeshFem, MeshIga, ... They can store arbitrary information, temporary 3D octrees, knots, nodes, ... And they own the elements.

Then, there is a container of cells for the calculation, each with a reference to a unified element interface, ElementCollection. So the cell has no clue if it is working on ElementFem or ElementIga.

Is the only difference that you want to forget about the MeshFem and MeshIga at some point? We could move the ownership of the elements from the mesh classes to the cell. This could make things like mesh adaption hard (but YAGNI...). My idea was to keep the element ownership at the meshes and provide references to the cell.

[TTitscher]

Settled those issues in a phone call. We should start implementing it in the proposed way. @joergfunger notes to keep the mesh file issues (85 to 92 or so) in mind. In an MPI context, we could have an interface like:

// single core
MeshFem mesh = MeshFem::ImportGeometryFromGmsh("mesh.msh");
mesh.AddDofInterpolationType(dofType, QUADRATIC); // similar to ConvertToInterpolationType()
mesh.ExportMeshFile("meshMpi.nutomesh", numProcs, partition(?));
// finished.

// other procs
MeshFem mesh = MeshFem::ImportFromMeshFile("meshMpi.nutumesh", rank());
auto cells = Cell::Create(mesh.ElementInterfaces(), integrand, integrationType, ...);
// solve...

This is not thought through and just a quick idea to show that the proposed idea does not clash with mesh files. They seem like a requirement for MPI. But IMO we should have the full functionality (for non MPI stuff) without those.

[joergfunger]

I think this issue is discussed further in the PDE project including the parallelization.