Unstructured meshes

Table of Contents

• Introduction
• Reading a mesh from a file
• Mesh configuration
• Public interface and basic usage
• Parallel iteration over mesh entities
• Writing a mesh and data to a file
• Example: Game of Life

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::Meshes::BuildConfigTags {
 
template<>
struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle >
{
   static constexpr bool enabled = true;
};
 
template<>
struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle >
{
   static constexpr bool enabled = true;
};
 
}  // namespace TNL::Meshes::BuildConfigTags

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;
}

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 );
   };
   const bool result = TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );
   return static_cast< int >( ! result );
}

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:

 
#include <TNL/Meshes/TypeResolver/resolveMeshType.h>
 
// Define the tag for the MeshTypeResolver configuration
struct MyConfigTag
{};
 
namespace TNL::Meshes::BuildConfigTags {
 
template<>
struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle >
{
   static constexpr bool enabled = true;
};
 
template<>
struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle >
{
   static constexpr bool enabled = true;
};
 
}  // namespace TNL::Meshes::BuildConfigTags
 
// 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;
}
 
int
main( int argc, char* argv[] )
{
   const std::string inputFileName = "example-triangles.vtu";
 
   auto wrapper = [ & ]( auto& reader, auto&& mesh ) -> bool
   {
      return task( mesh, inputFileName );
   };
   const bool result = TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );
   return static_cast< int >( ! result );
}

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::Meshes::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 SpaceDimension = Cell::dimension,
             typename Real = double,
             typename GlobalIndex = int,
             typename LocalIndex = short int >
   struct MeshConfig : public DefaultConfig< Cell, SpaceDimension, Real, GlobalIndex, LocalIndex >
   {
      static constexpr bool
      subentityStorage( int entityDimension, int subentityDimension )
      {
         return subentityDimension == 0 && entityDimension >= Cell::dimension - 1;
      }
   };
};
 
}  // namespace TNL::Meshes::BuildConfigTags

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:

• getSubentitiesCount
• getSubentityIndex
• getSuperentitiesCount
• getSuperentityIndex

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
   using DeviceMesh = TNL::Meshes::Mesh< typename HostMesh::Config, TNL::Devices::Cuda >;
   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 );

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 );
   VectorType 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" );
   };

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;

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::Meshes::BuildConfigTags {
 
// clang-format off
// disable all grids
template< int Dimension, typename Real, typename Device, typename Index >
struct GridTag< MyConfigTag, Grid< Dimension, Real, Device, Index > >
{ static constexpr bool enabled = false; };
 
// Meshes are enabled only for topologies explicitly listed below.
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Edge > { static constexpr bool enabled = true; };
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { static constexpr bool enabled = true; };
template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle > { static constexpr bool enabled = true; };
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Tetrahedron > { static constexpr bool enabled = true; };
//template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Hexahedron > { static constexpr bool enabled = true; };
 
// Meshes are enabled only for the space dimension equal to the cell dimension.
template< typename CellTopology, int SpaceDimension >
struct MeshSpaceDimensionTag< MyConfigTag, CellTopology, SpaceDimension >
{ enum { enabled = ( SpaceDimension == CellTopology::dimension ) }; };
 
// Meshes are enabled only for types explicitly listed below.
template<> struct MeshRealTag< MyConfigTag, float > { static constexpr bool enabled = false; };
template<> struct MeshRealTag< MyConfigTag, double > { static constexpr bool enabled = true; };
template<> struct MeshGlobalIndexTag< MyConfigTag, int > { static constexpr bool enabled = true; };
template<> struct MeshGlobalIndexTag< MyConfigTag, long int > { static constexpr bool enabled = false; };
template<> struct MeshLocalIndexTag< MyConfigTag, short int > { static constexpr bool enabled = true; };
// clang-format on
 
// Config tag specifying the MeshConfig template to use.
template<>
struct MeshConfigTemplateTag< MyConfigTag >
{
   template< typename Cell,
             int SpaceDimension = Cell::dimension,
             typename Real = double,
             typename GlobalIndex = int,
             typename LocalIndex = short int >
   struct MeshConfig : public DefaultConfig< Cell, SpaceDimension, Real, GlobalIndex, LocalIndex >
   {
      static constexpr bool
      subentityStorage( int entityDimension, int SubentityDimension )
      {
         return SubentityDimension == 0 && entityDimension >= Cell::dimension - 1;
      }
 
      static constexpr bool
      superentityStorage( int entityDimension, int SuperentityDimension )
      {
         //         return false;
         return ( entityDimension == 0 || entityDimension == Cell::dimension - 1 ) && SuperentityDimension == Cell::dimension;
      }
 
      static constexpr bool
      entityTagsStorage( int entityDimension )
      {
         //         return false;
         return entityDimension == 0 || entityDimension >= Cell::dimension - 1;
      }
 
      static constexpr bool
      dualGraphStorage()
      {
         return true;
      }
 
      static constexpr int dualGraphMinCommonVertices = 1;
   };
};
 
}  // namespace TNL::Meshes::BuildConfigTags
 
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 );
   VectorType f_out( cellsCount );
   f_in.setValue( 0 );
 
#if 1
   // generate random initial condition
   std::random_device dev;
   std::mt19937 rng( dev() );
   std::uniform_int_distribution<> dist( 0, 1 );
   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 );
 
   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 ) );
   };
   const bool result = Meshes::resolveAndLoadMesh< MyConfigTag, Devices::Host >( wrapper, inputFileName, inputFileFormat );
   return static_cast< int >( ! result );
}