Documentation: writing mesh tutorial

)

set( CPP_EXAMPLES

   ReadMeshExample

   MeshConfigurationExample

   MeshIterationExample

   ${COMMON_EXAMPLES}

)

set( CUDA_EXAMPLES

#include <TNL/Meshes/TypeResolver/resolveMeshType.h>

// Define the tag for the MeshTypeResolver configuration

struct MyConfigTag {};

//! [Configuration example]

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 = GlobalIndex >

   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

//! [Configuration example]

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

template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };

} // namespace BuildConfigTags

} // namespace Meshes

} // namespace TNL

int main( int argc, char* argv[] )

{

   const std::string inputFileName = "example-triangles.vtu";

   auto wrapper = [] ( auto& reader, auto&& mesh ) -> bool

   {

      return true;

   };

   return ! TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );

}

#include <TNL/Meshes/TypeResolver/resolveMeshType.h>

// Define the tag for the MeshTypeResolver configuration

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

template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };

} // namespace BuildConfigTags

} // namespace Meshes

} // namespace TNL

// Define the main task/function of the program

template< typename Mesh >

bool task( const Mesh& mesh )

{

   //! [getEntitiesCount]

   const int num_vertices = mesh.template getEntitiesCount< 0 >();

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

   //! [getEntitiesCount]

   // shut up warnings about unused variables

   (void) num_vertices;

   (void) num_cells;

   const int idx = num_vertices - 1;

   const int idx2 = num_cells - 1;

   //! [getEntity]

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

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

   //! [getEntity]

   // shut up warnings about unused variables

   (void) vert;

   (void) elem;

{

   //! [Iteration over subentities]

   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;

   }

   //! [Iteration over subentities]

}

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

   };

   return ! TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );

}

//! [config]

#include <TNL/Meshes/TypeResolver/resolveMeshType.h>

// Define the tag for the MeshTypeResolver configuration

struct MeshConfigTag {};

struct MyConfigTag {};

namespace TNL {

namespace Meshes {

namespace BuildConfigTags {

template<> struct MeshCellTopologyTag< MeshConfigTag, Topologies::Triangle > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MeshConfigTag, Topologies::Quadrangle > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Triangle > { enum { enabled = true }; };

template<> struct MeshCellTopologyTag< MyConfigTag, Topologies::Quadrangle > { enum { enabled = true }; };

} // namespace BuildConfigTags

} // namespace Meshes

} // namespace TNL

//! [config]

//! [task]

// Define the main task/function of the program

template< typename Mesh >

bool task( const Mesh& mesh, const std::string& inputFileName )

             << TNL::getType<Mesh>() << std::endl;

   return true;

}

//! [task]

//! [main]

int main( int argc, char* argv[] )

{

   const std::string inputFileName = "example-triangles.vtu";

   {

      return task( mesh, inputFileName );

   };

   return ! TNL::Meshes::resolveAndLoadMesh< MeshConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );

   return ! TNL::Meshes::resolveAndLoadMesh< MyConfigTag, TNL::Devices::Host >( wrapper, inputFileName, "auto" );

}

//! [main]

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.

\dontinclude ReadMeshExample.cpp

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](@ref TNL::Meshes::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:

\skip #include

\until // namespace TNL

\snippet ReadMeshExample.cpp config

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](@ref TNL::Meshes::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:

\skip Define the main task

\until }

\snippet ReadMeshExample.cpp task

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](@ref TNL::Meshes::resolveAndLoadMesh) function to resolve the mesh type and load the mesh from the file into the created object:

\skip int main

\until }

\until }

\snippet ReadMeshExample.cpp main

We need to specify two template parameters when calling `resolveAndLoadMesh`:

For completeness, the full example follows:

\includelineno ReadMeshExample.cpp

## Mesh configuration

## Mesh configuration   {#configuration}

The [Mesh](@ref TNL::Meshes::Mesh) class template is configurable via its first template parameter, `Config`.

By default, the \ref TNL::Meshes::DefaultConfig template is used.

For example, here we derive the `MeshConfig` template from the [DefaultConfig](@ref TNL::Meshes::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.

```cpp

// Define the tag for the MeshTypeResolver configuration

struct MyConfigTag {};

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 = GlobalIndex >

   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

```

\snippet MeshConfigurationExample.cpp Configuration example

## 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: \ref TNL::Meshes::Mesh, \ref TNL::Meshes::MeshEntity.

Here we describe only the basic member functions.

The main purpose of the [Mesh](@ref TNL::Meshes::Mesh) class template is to provide access to the mesh entities.

Firstly, there is a member function template called [getEntitiesCount](@ref TNL::Meshes::Mesh::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:

\snippet MeshIterationExample.cpp getEntitiesCount

Note that this member function and all other member functions presented below are marked as [\_\_cuda\_callable\_\_](tutorial_GeneralConcepts.html), 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](@ref TNL::Meshes::Mesh::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:

\snippet MeshIterationExample.cpp getEntity

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](@ref TNL::Meshes::MeshEntity) class template.

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

- [getSubentitiesCount](@ref TNL::Meshes::MeshEntity::getSubentitiesCount)

- [getSubentityIndex](@ref TNL::Meshes::MeshEntity::getSubentityIndex)

- [getSuperentitiesCount](@ref TNL::Meshes::MeshEntity::getSuperentitiesCount)

- [getSuperentityIndex](@ref TNL::Meshes::MeshEntity::getSuperentityIndex)

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

\snippet MeshIterationExample.cpp Iteration over subentities

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](#configuration).

## Parallel iteration over mesh entities

The [Mesh](@ref TNL::Meshes::Mesh) class template provides a simple interface for the parallel iteration over mesh entities of a specific dimension.

There are several member functions:

- [forAll](@ref TNL::Meshes::Mesh::forAll) -- iterates over all mesh entities of a specific dimension

- [forBoundary](@ref TNL::Meshes::Mesh::forBoundary) -- iterates over boundary mesh entities of a specific dimension

- [forInterior](@ref TNL::Meshes::Mesh::forInterior) -- iterates over interior (i.e., not boundary) mesh entities of a specific dimension

For distributed meshes there are two additional member functions:

- [forGhost](@ref TNL::Meshes::Mesh::forGhost) -- iterates over ghost mesh entities of a specific dimension

- [forLocal](@ref TNL::Meshes::Mesh::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:

```cpp

TODO

```

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](tutorial_Pointers.html) as follows:

```cpp

TODO

```

Alternatively, you can use a [SharedPointer](@ref TNL::Pointers::SharedPointer) instead of a [DevicePointer](@ref TNL::Pointers::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

## Example: Game of Life