MeshBuilder

[hiptmair]

Imagine the following use case, which is of practical relevance: The cells are affine, but the edge carry a parameterization that provides an exact description of the boundary. In this case just providing vertex coordinates and cell geometries is not enough.

Conversely, vertex positions can be deduced from cell geometries; hence, this piece of information is really redundant.

[craffael]

1. This is something I haven't considered yet. So far I have always assumed that geometrically speaking, a codim=n entity is always bounded by its codim=n+1 entities. But I know, this doesn't have to be the case, I've used such an approach myself in my master thesis.

2. Yes you are of course correct, this could be deduced from the element geometries.

A possible solution to both issues would be to replace the AddPoint() and AddElement() methods with the following method:

/**
 * @brief Insert a new entity explicitly into the mesh
 * @param ref_el The reference element of the entity
 * @param node_indices The zero based node-indices that make up this entity. 
                       If ref_el == RefEl::kPoint, this range should be empty.
 * @param geom  The geometric description of this entity. 
                If it is a nullptr, we try to deduce the geometric shape from entities with lower codimensions that have this entity as a sub-entity.
                If this is not, the entity will not be equipped with a geometry and calling Entity::Geometry() on it will return an exception.
 * @returns The unique, zero-based consecutive index of the entity (after the mesh has been constructed).
 */
size_type AddEntity(RefEl ref_el, ForwardRange<size_type> node_indices, std::unique_ptr<Geometry> geom)` 

[hiptmair]

I think, mesh generation facilities should not be included in the mesh class hierarchy and should be solely provided by MeshBuilders.

Therefore remove the constructor from hybrid2d::Mesh and include it into a MeshBuilder.

This comes at the price of allowing MeshBuilders unrestricted access to the data members of Mesh.

[craffael]

I agree with you that the constructor of the Mesh class should be private. We want the users to build meshes using the MeshBuilder and not the constructor of the hybrid2d mesh class. This means that the hybrid2d MeshBuilder must be a friend of the Mesh class.

From a users point of view I think it doesn't matter if the functionality that is currently stored in the constructor of the Mesh class is moved to the MeshBuilder. From the users perspective, its the package MeshBuilder + Mesh that work together in an expected way. The user shouldn't have to care about how the mesh internals work...

[craffael]

I just realized that we have different ideas about what a MeshBuilder is: For me the term MeshBuilder describes two concepts and I guess we should change the naming if we want to continue using them

1. Concept A: A Mesh Builder is an abstract interface that describes how a hybrid mesh can be constructed. This is currently the class lf::mesh::MeshBuilder. Then there are implementations of this interface such as lf::mesh::hybrid2d::MeshBuilder.

2. Concept B: A Mesh builder is any class that constructs a mesh from a set of parameters and uses the builder pattern. Currently class lf::mesh::TPTriagMeshBuilder would fulfill this concept. I would argue that these classes should not depend on the concrete implementation of the Mesh Manager. I.e. they should take as an input a lf::mesh::MeshBuilder and use it to construct the particular mesh. In this way they would also work with other mesh managers and if internals of the Mesh Manager are changed it is not necessary to change them too.

Maybe we should rename MeshBuilders of type A to MeshFactory...?

[hiptmair]

An idea for a very simple implementation of the generalized constructor for hybrid2d::Mesh:

• First initialize vector of Node entities and do not touch it anymore.
• Use std::map<pair<size_type,size_type>,vector<size_type>> to collect information about edges. First insert the supplied edges and then run through all cells to insert their edges. The mapped_type, a vector of indices, contains the indices of adjacent elements. After this associative container has been filled, run std::unique on all index vectors. Afterwards initialize vector of edges and do not touch it anymore.
• Finally initialize vector(s) of cells (Triangles and Quads) also using the information in the associative container for edges.

This initializes the Entity vectors in strict sequential order so that pointers to entities will always be valid. Moreover, it is a very straightforward implementation that dumps much of the algorithmic complexity onto std::map.

[craffael]

I agree this would be an option. However I'm usually not a big fan of using std::map for a large number of elements unless there is no other option. Although using std::map instead of std::sort gives similar (theoretical) algorithmic complexity, I prefer the latter for the following (performance) reasons:

• Because of the way red-black trees work, all the data that is stored in std::map must be stored on the heap and every insert operation leads to a new allocation. Moreover, the overhead to insert an element into the heap is non-trivial. Also memory locality is not guaranteed at all.
• std::sort relies on quick-sort which has worst-case complexity O(n^2), but on average it is O(n*log(n)) and in practice it's usually even much faster because some elements in the range are usually already correctly sorted.
• For sorting there exist powerful parallelized routines (see e.g. MKL) that scale more or less with the number of cpu cores. std::map is not thread safe and I think it's very hard to come-up with a thread-safe red-black tree implementation ;)

So if we want to write an algorithm that enables us to use entity pointers instead of "pointer to mesh" + index, I suggest to use a "two-stage sort"-like algorithm.

[hiptmair]

Of using std::map might sacrifice efficiency, but this is irrelevant, because the construction of meshes will hardly ever be a performance bottleneck unless incessant adaptive refinement and coarsening is employed in the context of timestepping. There are no plans to do this.