Template Numerical Library version develop:6514e4815
Unstructured meshes tutorial

Introduction

The Mesh class template is a data structure for conforming unstructured homogeneous meshes, which can be used as the fundamental data structure for numerical schemes based on finite volume or finite element methods. The abstract representation supports almost any cell shape which can be described by an entity topology. Currently there are common 2D quadrilateral, 3D hexahedron and arbitrarily dimensional simplex topologies built in the library. The implementation is highly configurable via templates of the C++ language, which allows to avoid the storage of unnecessary dynamic data. The internal memory layout is based on state–of–the–art sparse matrix formats, which are optimized for different hardware architectures in order to provide high performance computations. The DistributedMesh class template is an extended data structure based on Mesh, which allows to represent meshes decomposed into several subdomains for distributed computing using the Message Passing Interface (MPI).

Reading a mesh from a file

The most common way of mesh initialization is by reading a prepared input file created by an external program. TNL provides classes and functions for reading the common VTK, VTU and Netgen file formats.

The main difficulty is mapping the mesh included in the file to the correct C++ type, which can represent the mesh stored in the file. This can be done with the MeshTypeResolver class, which needs to be configured to enable the processing of the specific cell topologies, which we want our program to handle. For example, in the following code we enable loading of 2D triangular and quadrangular meshes:

#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
// Define the tag for the MeshTypeResolver configuration
struct MyConfigTag {};
namespace TNL {
namespace Meshes {
namespace BuildConfigTags {
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle > { enum { enabled = true }; };
} // namespace BuildConfigTags
} // namespace Meshes
} // namespace TNL
The main TNL namespace.
Definition: AtomicOperations.h:22

There are other build config tags which can be used to enable or disable specific types used in the mesh: RealType, GlobalIndexType and LocalIndexType. See the BuildConfigTags namespace for an overview of these tags.

Next, we can define the main task of our program as a templated function, which will be ultimately launched with the correct mesh type based on the input file. We can also use any number of additional parameters, such as the input file name:

// Define the main task/function of the program
template< typename Mesh >
bool task( const Mesh& mesh, const std::string& inputFileName )
{
std::cout << "The file '" << inputFileName << "' contains the following mesh: "
<< TNL::getType<Mesh>() << std::endl;
return true;
}
T endl(T... args)

Of course in practice, the function would be much more complex than this example, where we just print the file name and some textual representation of the mesh to the standard output.

Finally, we define the main function, which sets the input parameters (hard-coded in this example) and calls the resolveAndLoadMesh function to resolve the mesh type and load the mesh from the file into the created object:

int main( int argc, char* argv[] )
{
const std::string inputFileName = "example-triangles.vtu";
auto wrapper = [&] ( auto& reader, auto&& mesh ) -> bool
{
return task( mesh, inputFileName );
};
return ! TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );
}

We need to specify two template parameters when calling resolveAndLoadMesh:

  1. our build config tag (MeshConfigTag in this example),
  2. and the device where the mesh should be allocated.

Then we pass the the function which should be called with the initialized mesh, the input file name, and the input file format ("auto" means auto-detection based on the file name). In order to show the flexibility of passing other parameters to our main task function as defined above, we suggest to implement a wrapper lambda function (called wrapper in the example), which captures the relevant variables and forwards them to the task.

The return value of the resolveAndLoadMesh function is a boolean value representing the success (true) or failure (false) of the whole function call chain. Hence, the return type of both wrapper and task needs to be bool as well.

For completeness, the full example follows:

1
2#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
3
4// Define the tag for the MeshTypeResolver configuration
5struct MyConfigTag {};
6
7namespace TNL {
8namespace Meshes {
9namespace BuildConfigTags {
10
11template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };
12template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle > { enum { enabled = true }; };
13
14} // namespace BuildConfigTags
15} // namespace Meshes
16} // namespace TNL
18
20// Define the main task/function of the program
21template< typename Mesh >
22bool task( const Mesh& mesh, const std::string& inputFileName )
23{
24 std::cout << "The file '" << inputFileName << "' contains the following mesh: "
25 << TNL::getType<Mesh>() << std::endl;
26 return true;
27}
29
31int main( int argc, char* argv[] )
32{
33 const std::string inputFileName = "example-triangles.vtu";
34
35 auto wrapper = [&] ( auto& reader, auto&& mesh ) -> bool
36 {
37 return task( mesh, inputFileName );
38 };
39 return ! TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );
40}

Mesh configuration

The Mesh class template is configurable via its first template parameter, Config. By default, the TNL::Meshes::DefaultConfig template is used. Alternative, user-specified configuration templates can be specified by defining the mesh configuration as the MeshConfig template in the MeshConfigTemplateTag build config tag specialization. For example, here we derive the MeshConfig template from the DefaultConfig template and override the subentityStorage member function to store only those subentity incidence matrices, where the subentity dimension is 0 and the other dimension is at least $D-1$. Hence, only faces and cells will be able to access their subvertices and there will be no other links from entities to their subentities.

namespace TNL {
namespace Meshes {
namespace BuildConfigTags {
// Create a template specialization of the tag specifying the MeshConfig template to use as the Config parameter for the mesh.
template<>
struct MeshConfigTemplateTag< MyConfigTag >
{
template< typename Cell,
int WorldDimension = Cell::dimension,
typename Real = double,
typename GlobalIndex = int,
typename LocalIndex = short int >
struct MeshConfig
: public DefaultConfig< Cell, WorldDimension, Real, GlobalIndex, LocalIndex >
{
template< typename EntityTopology >
static constexpr bool subentityStorage( EntityTopology, int SubentityDimension )
{
return SubentityDimension == 0 && EntityTopology::dimension >= Cell::dimension - 1;
}
};
};
} // namespace BuildConfigTags
} // namespace Meshes
} // namespace TNL

Public interface and basic usage

The whole public interface of the unstructured mesh and its mesh entity class can be found in the reference manual: TNL::Meshes::Mesh, TNL::Meshes::MeshEntity. Here we describe only the basic member functions.

The main purpose of the Mesh class template is to provide access to the mesh entities. Firstly, there is a member function template called getEntitiesCount which returns the number of entities of the dimension specified as the template argument. Given a mesh instance denoted as mesh, it can be used like this:

const int num_vertices = mesh.template getEntitiesCount< 0 >();
const int num_cells = mesh.template getEntitiesCount< Mesh::getMeshDimension() >();

Note that this member function and all other member functions presented below are marked as __cuda_callable__, so they can be called from usual host functions as well as CUDA kernels.

The entity of given dimension and index can be accessed via a member function template called getEntity. Again, the entity dimension is specified as a template argument and the index is specified as a method argument. The getEntity member function does not provide a reference access to an entity stored in the mesh, but each entity is created on demand and contains only a pointer to the mesh and the supplied entity index. Hence, the mesh entity is kind of a proxy object where all member functions call just an appropriate member function via the mesh pointer. The getEntity member function can be used like this:

typename Mesh::Vertex vert = mesh.template getEntity< 0 >( idx );
typename Mesh::Cell elem = mesh.template getEntity< Mesh::getMeshDimension() >( idx2 );

Here we assume that idx < num_vertices and idx2 < num_cells. Note that both Mesh::Vertex and Mesh::Cell are specific instances of the MeshEntity class template.

The information about the subentities and superentities can be accessed via the following member functions:

For example, they can be combined with the getEntity member function to iterate over all subvertices of a specific cell:

typename Mesh::Cell elem = mesh.template getEntity< Mesh::getMeshDimension() >( idx2 );
const int n_subvert = elem.template getSubentitiesCount< 0 >();
for( int v = 0; v < n_subvert; v++ ) {
const int v_idx = elem.template getSubentityIndex< 0 >( v );
typename Mesh::Vertex vert = mesh.template getEntity< 0 >( v_idx );
// [Do some work...]
(void) vert;
}

The iteration over superentities adjacent to an entity is very similar and left as an exercise for the reader.

Note that the implementations of all templated member functions providing access to subentities and superentities contain a static_assert expression which checks if the requested subentities or superentities are stored in the mesh according to its configuration.

Parallel iteration over mesh entities

The Mesh class template provides a simple interface for the parallel iteration over mesh entities of a specific dimension. There are several member functions:

  • forAll – iterates over all mesh entities of a specific dimension
  • forBoundary – iterates over boundary mesh entities of a specific dimension
  • forInterior – iterates over interior (i.e., not boundary) mesh entities of a specific dimension

For distributed meshes there are two additional member functions:

  • forGhost – iterates over ghost mesh entities of a specific dimension
  • forLocal – iterates over local (i.e., not ghost) mesh entities of a specific dimension

All of these member functions have the same interface: they take one parameter, which should be a functor, such as a lambda expression $f$ that is called as $f(i)$, where $i$ is the mesh entity index in the current iteration. Remember that the iteration is performed in parallel, so all calls to the functor must be independent since they can be executed in any order.

Note that only the mesh entity index is passed to the functor, it does not get the mesh entity object or even (a reference to) the mesh. All additional information needed by the functor must be handled manually, e.g. via a lambda capture.

For example, the iteration over cells on a mesh allocated on the host can be done as follows:

auto kernel = [&mesh] ( typename Mesh::GlobalIndexType i ) mutable
{
typename Mesh::Cell elem = mesh.template getEntity< Mesh::getMeshDimension() >( i );
// [Do some work with the current cell `elem`...]
(void) elem;
};
mesh.template forAll< Mesh::getMeshDimension() >( kernel );

The parallel iteration is more complicated for meshes allocated on a GPU, since the lambda expression needs to capture a pointer to the copy of the mesh, which is allocated on the right device. This can be achieved with a smart pointer as follows:

// Copy the mesh from host to the device
DeviceMesh deviceMesh = hostMesh;
// Create a device pointer to the device mesh
TNL::Pointers::DevicePointer< const DeviceMesh > meshDevicePointer( deviceMesh );
const DeviceMesh* meshPointer = &meshDevicePointer.template getData< typename DeviceMesh::DeviceType >();
// Define and execute the kernel on the device
auto kernel = [meshPointer] __cuda_callable__ ( typename DeviceMesh::GlobalIndexType i ) mutable
{
typename DeviceMesh::Cell elem = meshPointer->template getEntity< DeviceMesh::getMeshDimension() >( i );
// [Do some work with the current cell `elem`...]
(void) elem;
};
deviceMesh.template forAll< DeviceMesh::getMeshDimension() >( kernel );
Definition: Mesh.h:70
The DevicePointer is like SharedPointer, except it takes an existing host object - there is no call t...
Definition: DevicePointer.h:51

Alternatively, you can use a SharedPointer instead of a DevicePointer to allocate the mesh, but it does not allow to bind to an object which has already been created outside of the SharedPointer.

Writing a mesh and data to a file

Numerical simulations typically produce results which can be interpreted as mesh functions or fields. In C++ they can be stored simply as arrays or vectors with the appropriate size. For example, here we create two arrays f_in and f_out, which represent the input and output state of an iterative algorithm (f_in and f_out will be swapped after each iteration):

using Index = typename Mesh::GlobalIndexType;
const Index cellsCount = mesh.template getEntitiesCount< Mesh::getMeshDimension() >();
using VectorType = Containers::Vector< std::uint8_t, typename Mesh::DeviceType, Index >;
VectorType f_in( cellsCount ), f_out( cellsCount );
f_in.setValue( 0 );

Note that here we used std::uint8_t as the value type. The following value types are supported for the output into VTK file formats: std::int8_t, std::uint8_t, std::int16_t, std::uint16_t, std::int32_t, std::uint32_t, std::int64_t, std::uint64_t, float, double.

The output into a specific file format can be done with an appropriate writer class, see TNL::Meshes::Writers. For example, using VTUWriter for the .vtu file format:

auto make_snapshot = [&] ( Index iteration )
{
const std::string filePath = "GoL." + std::to_string(iteration) + ".vtu";
std::ofstream file( filePath );
using Writer = Meshes::Writers::VTUWriter< Mesh >;
Writer writer( file );
writer.writeMetadata( iteration, iteration );
writer.template writeEntities< Mesh::getMeshDimension() >( mesh );
writer.writeCellData( f_in, "function values" );
};
T to_string(T... args)

Note that this writer supports writing metadata (iteration index and time level value), then we call writeEntities to write the mesh cells and writeCellData to write the mesh function values. The writeCellData call can be repeated multiple times for different mesh functions that should be included in the snapshot.

Then we can take the snapshot of the initial state,

// write initial state
make_snapshot( 0 );

and similarly use make_snapshot in the iteration loop.

Example: Game of Life

In this example we will show how to implement the Conway's Game of Life using the Mesh class template. Although the game is usually implemented on structured grids rather than unstructured meshes, it will nicely illustrate how the building blocks for a numerical simulation are connected together.

The kernel of the Game of Life can be implemented as follows:

auto kernel = [f_in_view, f_out_view, meshPointer] __cuda_callable__ ( Index i ) mutable
{
// sum values of the function on the neighbor cells
typename VectorType::RealType sum = 0;
for( Index n = 0; n < meshPointer->getCellNeighborsCount( i ); n++ ) {
const Index neighbor = meshPointer->getCellNeighborIndex( i, n );
sum += f_in_view[ neighbor ];
}
const bool live = f_in_view[ i ];
// Conway's rules for square grid
if( live ) {
// any live cell with less than two live neighbors dies
if( sum < 2 )
f_out_view[ i ] = 0;
// any live cell with two or three live neighbors survives
else if( sum < 4 )
f_out_view[ i ] = 1;
// any live cell with more than three live neighbors dies
else
f_out_view[ i ] = 0;
}
else {
// any dead cell with exactly three live neighbors becomes a live cell
if( sum == 3 )
f_out_view[ i ] = 1;
// any other dead cell remains dead
else
f_out_view[ i ] = 0;
}
};

The kernel function takes f_in_view (the input state of the current iteration) and for the $i$-th cell sums the values of the neighbor cells, which are accessed using the dual graph – see TNL::Meshes::Mesh::getCellNeighborsCount and TNL::Meshes::Mesh::getCellNeighborIndex. Then it writes the resulting state of the $i$-th cell into f_out_view according to Conway's rules for a square grid:

  • any live cell with less than two live neighbors dies,
  • any live cell with two or three live neighbors survives,
  • any live cell with more than three live neighbors dies,
  • any dead cell with exactly three live neighbors becomes a live cell,
  • and any other dead cell remains dead.

The kernel function is evaluated for all cells in the mesh, followed by swapping f_in and f_out (including their views), writing the output into a VTU file and checking if this was the last iteration:

// iterate over all cells
mesh.template forAll< Mesh::getMeshDimension() >( kernel );
// swap input and output arrays
f_in.swap( f_out );
// remember to update the views!
f_in_view.bind( f_in.getView() );
f_out_view.bind( f_out.getView() );
// write output
make_snapshot( iteration );
// check if finished
all_done = max( f_in ) == 0 || iteration > max_iter || f_in == f_out;
T max(T... args)

The remaining pieces needed for the implementation have either been already presented on this page, or they are left as an exercise to the reader. For the sake of completeness, we include the full example below.

#include <random>
#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
#include <TNL/Meshes/Geometry/getEntityCenter.h>
#include <TNL/Meshes/Writers/VTUWriter.h>
using namespace TNL;
struct MyConfigTag {};
namespace TNL {
namespace Meshes {
namespace BuildConfigTags {
// disable all grids
template< int Dimension, typename Real, typename Device, typename Index >
struct GridTag< MyConfigTag, Grid< Dimension, Real, Device, Index > >
{ enum { enabled = false }; };
// Meshes are enabled only for topologies explicitly listed below.
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Edge > { enum { enabled = true }; };
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle > { enum { enabled = true }; };
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Tetrahedron > { enum { enabled = true }; };
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Hexahedron > { enum { enabled = true }; };
// Meshes are enabled only for the world dimension equal to the cell dimension.
template< typename CellTopology, int WorldDimension >
struct MeshWorldDimensionTag< MyConfigTag, CellTopology, WorldDimension >
{ enum { enabled = ( WorldDimension == CellTopology::dimension ) }; };
// Meshes are enabled only for types explicitly listed below.
template<> struct MeshRealTag< MyConfigTag, float > { enum { enabled = false }; };
template<> struct MeshRealTag< MyConfigTag, double > { enum { enabled = true }; };
template<> struct MeshGlobalIndexTag< MyConfigTag, int > { enum { enabled = true }; };
template<> struct MeshGlobalIndexTag< MyConfigTag, long int > { enum { enabled = false }; };
template<> struct MeshLocalIndexTag< MyConfigTag, short int > { enum { enabled = true }; };
// Config tag specifying the MeshConfig template to use.
template<>
struct MeshConfigTemplateTag< MyConfigTag >
{
template< typename Cell,
int WorldDimension = Cell::dimension,
typename Real = double,
typename GlobalIndex = int,
typename LocalIndex = short int >
struct MeshConfig
: public DefaultConfig< Cell, WorldDimension, Real, GlobalIndex, LocalIndex >
{
template< typename EntityTopology >
static constexpr bool subentityStorage( EntityTopology, int SubentityDimension )
{
return SubentityDimension == 0 && EntityTopology::dimension >= Cell::dimension - 1;
}
template< typename EntityTopology >
static constexpr bool superentityStorage( EntityTopology, int SuperentityDimension )
{
// return false;
return (EntityTopology::dimension == 0 || EntityTopology::dimension == Cell::dimension - 1) && SuperentityDimension == Cell::dimension;
}
template< typename EntityTopology >
static constexpr bool entityTagsStorage( EntityTopology )
{
// return false;
return EntityTopology::dimension == 0 || EntityTopology::dimension >= Cell::dimension - 1;
}
static constexpr bool dualGraphStorage()
{
return true;
}
static constexpr int dualGraphMinCommonVertices = 1;
};
};
} // namespace BuildConfigTags
} // namespace Meshes
} // namespace TNL
template< typename Mesh >
bool runGameOfLife( const Mesh& mesh )
{
using Index = typename Mesh::GlobalIndexType;
const Index cellsCount = mesh.template getEntitiesCount< Mesh::getMeshDimension() >();
using VectorType = Containers::Vector< std::uint8_t, typename Mesh::DeviceType, Index >;
VectorType f_in( cellsCount ), f_out( cellsCount );
f_in.setValue( 0 );
#if 1
// generate random initial condition
std::mt19937 rng(dev());
for( Index i = 0; i < cellsCount; i++ )
f_in[ i ] = dist(rng);
const Index max_iter = 100;
#else
// find a reference cell - the one closest to the point
typename Mesh::PointType p = {0.5, 0.5};
typename Mesh::RealType dist = 1e5;
Index reference_cell = 0;
for( Index i = 0; i < cellsCount; i++ ) {
const auto cell = mesh.template getEntity< Mesh::getMeshDimension() >( i );
const auto c = getEntityCenter( mesh, cell );
const auto d = TNL::l2Norm( c - p );
if( d < dist ) {
reference_cell = i;
dist = d;
}
}
// R-pentomino (stabilizes after 1103 iterations)
const Index max_iter = 1103;
f_in[reference_cell] = 1;
Index n1 = mesh.getCellNeighborIndex(reference_cell,1); // bottom
Index n2 = mesh.getCellNeighborIndex(reference_cell,2); // left
Index n3 = mesh.getCellNeighborIndex(reference_cell,5); // top
Index n4 = mesh.getCellNeighborIndex(reference_cell,6); // top-right
f_in[n1] = 1;
f_in[n2] = 1;
f_in[n3] = 1;
f_in[n4] = 1;
/*
// Acorn (stabilizes after 5206 iterations)
const Index max_iter = 5206;
f_in[reference_cell] = 1;
Index n1 = mesh.getCellNeighborIndex(reference_cell,4);
f_in[n1] = 1;
Index s1 = mesh.getCellNeighborIndex(n1,4);
Index s2 = mesh.getCellNeighborIndex(s1,4);
Index n2 = mesh.getCellNeighborIndex(s2,4);
f_in[n2] = 1;
Index n3 = mesh.getCellNeighborIndex(n2,4);
f_in[n3] = 1;
Index n4 = mesh.getCellNeighborIndex(n3,4);
f_in[n4] = 1;
f_in[mesh.getCellNeighborIndex(s2,5)] = 1;
f_in[mesh.getCellNeighborIndex(mesh.getCellNeighborIndex(n1,5),5)] = 1;
*/
#endif
auto make_snapshot = [&] ( Index iteration )
{
const std::string filePath = "GoL." + std::to_string(iteration) + ".vtu";
std::ofstream file( filePath );
using Writer = Meshes::Writers::VTUWriter< Mesh >;
Writer writer( file );
writer.writeMetadata( iteration, iteration );
writer.template writeEntities< Mesh::getMeshDimension() >( mesh );
writer.writeCellData( f_in, "function values" );
};
// write initial state
make_snapshot( 0 );
// captures for the iteration kernel
auto f_in_view = f_in.getConstView();
auto f_out_view = f_out.getView();
Pointers::DevicePointer< const Mesh > meshDevicePointer( mesh );
const Mesh* meshPointer = &meshDevicePointer.template getData< typename Mesh::DeviceType >();
bool all_done = false;
Index iteration = 0;
do {
iteration++;
std::cout << "Computing iteration " << iteration << "..." << std::endl;
auto kernel = [f_in_view, f_out_view, meshPointer] __cuda_callable__ ( Index i ) mutable
{
// sum values of the function on the neighbor cells
typename VectorType::RealType sum = 0;
for( Index n = 0; n < meshPointer->getCellNeighborsCount( i ); n++ ) {
const Index neighbor = meshPointer->getCellNeighborIndex( i, n );
sum += f_in_view[ neighbor ];
}
const bool live = f_in_view[ i ];
// Conway's rules for square grid
if( live ) {
// any live cell with less than two live neighbors dies
if( sum < 2 )
f_out_view[ i ] = 0;
// any live cell with two or three live neighbors survives
else if( sum < 4 )
f_out_view[ i ] = 1;
// any live cell with more than three live neighbors dies
else
f_out_view[ i ] = 0;
}
else {
// any dead cell with exactly three live neighbors becomes a live cell
if( sum == 3 )
f_out_view[ i ] = 1;
// any other dead cell remains dead
else
f_out_view[ i ] = 0;
}
};
// iterate over all cells
mesh.template forAll< Mesh::getMeshDimension() >( kernel );
// swap input and output arrays
f_in.swap( f_out );
// remember to update the views!
f_in_view.bind( f_in.getView() );
f_out_view.bind( f_out.getView() );
// write output
make_snapshot( iteration );
// check if finished
all_done = max( f_in ) == 0 || iteration > max_iter || f_in == f_out;
}
while( all_done == false );
return true;
}
int main( int argc, char* argv[] )
{
const std::string inputFileName = "grid-100x100.vtu";
const std::string inputFileFormat = "auto";
auto wrapper = [&] ( auto& reader, auto&& mesh ) -> bool
{
using MeshType = std::decay_t< decltype(mesh) >;
return runGameOfLife( std::forward<MeshType>(mesh) );
};
return ! Meshes::resolveAndLoadMesh< MyConfigTag, Devices::Host >( wrapper, inputFileName, inputFileFormat );
}