Feel++

Using computational meshes

Introduction

Feel++ provides some tools to manipulate meshes. Here is a basic example that shows how to generate a mesh for a square geometry

Excerpt from codes/mymesh.cpp
Unresolved directive in README.adoc - include::../../../../codes/03-mymesh.cpp[tag=mesh]

As always, we initialize the {feelpp} environment (see section [FirstApp] ). The unitSquare() will generate a mesh for a square geometry. Feel++ provides several functions to automate the GMSH mesh generation for different topologies. These functions will create a geometry file .geo and a mesh file .msh. We can visualize them in {uri-gmsh-www}[Gmsh].

$ gmsh <entity_name>.msh

Finally we use the exporter() (see \ref Exporter) function to export the mesh for post processing. It will create by default a Paraview format file .sos and an Ensight format file .case.

$ paraview <app_name>.case

In this section, we present some of the mesh definition and manipulation tools provided by Feel++. For more information you can also read {uri-gmsh-manual}[the Gmsh manual].

Basic Meshes

There is a list of basic geometries you can automatically generate with Feel++ library.

Table 1. Table of function generation basic primitive mesh

Feel++ function

Dim

Description

unitSegment()

1

Build a mesh of the unit segment \$[0,1\$]

unitSquare()

2

Build a mesh of the unit square \$[0,1\$^2] using triangles

unitCircle()

2

Build a mesh of the unit circle using triangles

unitHypercube()

3

Build a mesh of the unit hypercube \$[0,1\$^3] using tetrahedrons

unitSphere()

3

Build a mesh of the unit sphere using tetrahedrons

Example
auto mesh = unitSquare();

Loading Meshes

loadMesh

You can use this function to load mesh from different formats

in the case of a .geo file, {feelpp} automatically generate a mesh data structure on this geometrical structure.
Interface
mesh_ptrtype loadMesh(_mesh, _filename, _refine, _update, physical_are_elementary_regions);

Required Parameters:

  • _mesh a mesh data structure.

Optional Parameters:

Name Type Description Default Value Option

_hsize

double

characteristic size of the mesh. This option will edit the .geo file and change the variable h if defined

0.1

gmsh.hsize

_geo_variables_list

string

Set a list of variable that may be defined in a .geo file

""

gmsh.geo-variables-list

_filename

string

filename with extension

feel.geo

gmsh.filename

_depends

string

list of files (separated by , or ;) on which gmsh.filename depends

""

gmsh.depends

_refine

boolean

optionally refine with \p refine levels the mesh.

0.

gmsh.refine

_update

integer

update the mesh data structure (build internal faces and edges)

true

_physical_are_elementary_regions

boolean

to load specific meshes formats.

false.

gmsh.physical_are_elementary_regions

_straighten

boolean

in case of curvilinear elements, straighten the elements which are not touching with a face the boundary of the domain

true

gmsh.straighten

_partitioner

The file you want to load has to be in an appropriate repository. {feelpp} looks for the files in the following directories (in this order):

  • current path

  • paths that went through changeRepository(), it means that we look for example into the path from which the executable was run

  • localGeoRepository() which is usually "$HOME/feel/geo" (Environment )

  • systemGeoRepository() which is usually "$FEELPP_DIR/share/feel/geo" (Environment)

Examples

Load a mesh data structure from the file $HOME/feel/mymesh.msh.

auto mesh = loadMesh(_mesh=new mesh_type,
                     _filename="mymesh.msh");

Load a geometric structure from the file ./mygeo.geo and automatically create a mesh data structure.

auto mesh = loadMesh(_mesh=new mesh_type,
                     _filename="mygeo.geo");

Create a mesh data structure from the file ./feel.geo.

auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );

loadGMSHMesh

In order to load only .msh file, you can also use the loadGMSHMesh.

Interface:

mesh_ptrtype loadGMSHMesh(_mesh, _filename, _refine, _update, _physical_are_elementary_regions);

Required Parameters:

  • _mesh a mesh data structure.

  • _filename filename with extension.

Optional Parameters:

  • _refine optionally refine with \p refine levels the mesh. - Default =0

  • _update update the mesh data structure (build internal faces and edges).

    • Default =true

  • _physical_are_elementary_regions to load specific meshes formats.

    • Default = false

The file you want to load has to be in an appropriate repository. See LoadMesh.

Examples

From doc/manual/heatns.cpp

 mesh_ptrtype mesh = loadGMSHMesh( _mesh=new mesh_type,
                                   _filename="piece.msh",
                                   _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER );

From applications/check/check.cpp

mesh = loadGMSHMesh( _mesh=new mesh_type,
                     _filename=soption("filename"),
                     _rebuild_partitions=(Environment::worldComm().size() > 1),
                     _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );

Creating Meshes

createGMSHMesh

Interface
mesh_ptrtype createGMSHMesh(_mesh, _desc, _h, _order, _parametricnodes, _refine, _update, _force_rebuild, _physical_are_elementary_regions);

Required Parameters:

  • _mesh mesh data structure.

  • _desc descprition. See further.

Optional Parameters:

  • _h characteristic size.

    • Default = 0.1

  • _order order.

    • Default = 1

  • _parametricnodes

    • Default = 0

  • _refine optionally refine with \p refine levels the mesh.

    • Default =0

  • _update update the mesh data structure (build internal faces and edges).

    • Default =true

  • _force_rebuild rebuild mesh if already exists.

    • Default = false

  • _physical_are_elementary_regions to load specific meshes formats.

    • Default = false

To generate your mesh you need a description parameter. This one can be create by one the two following function.

geo

Use this function to create a description from a .geo file.

Interface
gmsh_ptrtype geo(_filename, _h, _dim, _order, _files_path);

Required Parameters:

  • filename: file to load.

Optional Parameters:

  • _h characteristic size of the mesh.

    • Default = 0.1.

  • _dim dimension.

    • Default = 3.

  • _order order.

    • Default = 1.

  • _files_path path to the file.

    • Default = localGeoRepository().

The file you want to load has to be in an appropriate repository. See LoadMesh.

Example

From doc/manual/heat/ground.cpp

mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=geo( _filename="ground.geo",
                                  _dim=2,
                                  _order=1,
                                  _h=meshSize ) );
Excerpt from doc/manual/fd/penalisation.cpp
mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=geo( _filename=File_Mesh,
                                  _dim=Dim,
                                  _h=doption(_name="gmsh.hsize"),
                                  _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER );

domain

Use this function to generate a simple geometrical domain from parameters.

Interface
gmsh_ptrtype domain(_name, _shape, _h, _dim, _order, _convex, \
                    _addmidpoint, _xmin, _xmax, _ymin, _ymax, _zmin, _zmax);

Required Parameters:

  • _name name of the file that will ge generated without extension.

  • _shape shape of the domain to be generated (simplex or hypercube).

Optional Parameters:

  • _h characteristic size of the mesh.

    • Default = 0.1

  • _dim dimension of the domain.

    • Default = 2

  • _order order of the geometry.

    • Default = 1

  • _convex type of convex used to mesh the domain.

    • Default = simplex

  • _addmidpoint add middle point.

    • Default = true

  • _xmin minimum x coordinate.

    • Default = 0

  • _xmax maximum x coordinate.

    • Default = 1

  • _ymin minimum y coordinate.

    • Default = 0

  • _ymax maximum y coordinate.

    • Default = 1.

  • _zmin minimum z coordinate.

    • Default = 0

  • _zmax maximum z coordinate.

    • Default = 1

Example

From doc/manual/laplacian/laplacian.ccp

mesh_ptrtype mesh = createGMSHMesh( _mesh=new mesh_type,
                                    _desc=domain( _name=( boost::format( "%1%-%2%" ) % shape % Dim ).str() ,
                                                  _usenames=true,
                                                  _shape=shape,
                                                  _h=meshSize,
                                                  _xmin=-1,
                                                  _ymin=-1 ) );

From doc/manual/stokes/stokes.cpp

mesh = createGMSHMesh( _mesh=new mesh_type,
                       _desc=domain( _name=(boost::format("%1%-%2%-%3%")%"hypercube"%convex_type().dimension()%1).str() ,
                                     _shape="hypercube",
                                     _dim=convex_type().dimension(),
                                     _h=meshSize ) );

From doc/manual/solid/beam.cpp

mesh_ptrtype mesh = createGMSHMesh( _mesh=new mesh_type,
                                    _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK,
                                    _desc=domain( _name=( boost::format( "beam-%1%" ) % nDim ).str(),
                                                  _shape="hypercube",
                                                  _xmin=0., _xmax=0.351,
                                                  _ymin=0., _ymax=0.02,
                                                  _zmin=0., _zmax=0.02,
                                                  _h=meshSize ) );

Mesh iterators

Feel++ mesh data structure allows to iterate over its entities: elements, faces, edges and points.

The following table describes free-functions that allow to define mesh region over which to operate. MeshType denote the type of mesh passed to the free functions in the table.

MeshType can be a pointer, a shared_pointer or a reference to a mesh type.

For example :

auto mesh = loadMesh( _mesh=Mesh<Simplex<2>>);
auto r1 = elements(mesh); // OK
auto r2 = elements(*mesh); // OK
Table 2. Table of mesh iterators
Type Function Description

elements_t<MeshType>

elements(mesh)

All the elements of a mesh

markedelements_t<MeshType>

markedelements(mesh, id)

All the elements marked by marked id

boundaryelements_t<MeshType>

boundaryelements(mesh)

All the elements of the mesh which share a face with the boundary of the mesh.

internalelements_t<MeshType>

internalelements(mesh)

All the elements of the mesh which share a face with the boundary of the mesh.

pid_faces_t<MeshType>

faces(mesh)

All the faces of the mesh.

markedfaces_t<MeshType>

markedfaces(mesh)

All the faces of the mesh which are marked.

boundaryfaces_t<MeshType>

boundaryfaces(mesh)

All elements that own a topological dimension one below the mesh. For example, if you mesh is a 2D one, boundaryfaces(mesh) will return all the lines (because of dimension 2-1=1). These elements which have one dimension less, are corresponding to the boundary faces.

internalfaces_t<MeshType>

internalelements(mesh)

All the elements of the mesh which are stricly within the domain that is to say they do not share a face with the boundary.

edges_t<MeshType>

edges(mesh)

All the edges of the mesh.

boundaryedges_t<MeshType>

boundaryedges(mesh)

All boundary edges of the mesh.

points_t<MeshType>

points(mesh)

All the points of the mesh.

markedpoints_t<MeshType>

markedpoints(mesh,id)

All the points marked id of mesh.

boundarypoints_t<MeshType>

boundarypoints(mesh)

All boundary points of the mesh.

internalpoints_t<MeshType>

internalpoints(mesh)

All internal points of the mesh(not on the boundary)

Here are some examples on how to use these functionSpace

auto mesh = ...;

auto r1 = elements(mesh);
// iterate over the set of elements local to the process(no ghost cell selected, see next section)
for ( auto const&  e : r2 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r2 = markedelements(mesh,"iron");
// iterate over the set of elements marked iron in the mesh
for ( auto const&  e : r2 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r3 = boundaryfaces(mesh);
// iterate over the set of faces on the boundary of the mesh
for ( auto const&  e : r3 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

auto r4 = markededges(mesh,"line");
// iterate over the set of edges marked "line" in the mesh
for ( auto const&  e : r4 )
{
  auto const& elt = unwrap_ref( e );
  // work with element elt
  ...
}

Extended set of entities

Feel++ allows also to select an extended sets of entities from the mesh, you can extract entities which belongs to the local process but also ghost entities which satisfy the same property as the local ones.

Actually you can select both or one one of them thanks to the enum data structure entity_process_t which provides the following options

entity_process_t

Description

LOCAL_ONLY

only local entities

GHOST_ONLY

only ghost entities

ALL

both local and ghost entities

Type

Function

Description

ext_elements_t<MeshType>

elements(mesh,entity_process_t)

all the elements of mesh associated to entity_process_t.

ext_elements_t<MeshType>

markedelements(mesh, id, entity_process_t)

all the elements marked id of mesh associated to entity_process_t.

ext_faces_t<MeshType>

faces(mesh,entity_process_t)

all the faces of mesh associated to entity_process_t.

ext_faces_t<MeshType>

markedfaces(mesh, id, entity_process_t)

all the faces marked id of mesh associated to entity_process_t.

ext_edges_t<MeshType>

edges(mesh,entity_process_t)

all the edges of mesh associated to entity_process_t.

ext_edges_t<MeshType>

markededges(mesh, id, entity_process_t)

all the edges marked id of mesh associated to entity_process_t.

The type of the object returned for an entity is always the same, for elements it is ext_elements_t<MeshType> whether the elements are marked or not. The reason is that in fact we have to create a temporary data structure embedded in the range object that stores a reference to the elements which are selected.

Here is how to select both local and ghost elements from a Mesh

auto mesh =...;
auto r = elements(mesh,entity_process_t::ALL);
for (auto const& e : r )
{
  // do something on the local and ghost element
  ...
  // do something special on ghost cells
  if ( unwrap_ref(e).isGhostCell() )
  {...}
}

Concatenate sets of entities

Denote \$\mathcal{E}_{1}, \ldots ,\mathcal{E}_{n}\$ \$n\$ disjoints sets of the same type of entities (eg elements, faces,edges or points), \$\cup_{i=1}^{n} \mathcal{E}_i\$ with \$\cap_{i=0}^{n} \mathcal{E}_i = \emptyset\$.

We wish to concatenate these \$n\$ sets. To this end, we use concatenate which takes an arbitrary number of disjoints sets.

#include <feel/feelmesh/concatenate.hpp>
...
auto E_1 = internalfaces(mesh);
auto E_2 = markedfaces(mesh,"Gamma_1");
auto E_3 = markedfaces(mesh,"Gamma_2");
auto newset = concatenate( E_1, E_2, E_3 );
cout << "measure of newset = " << integrate(_range=newset, _expr=cst(1.)).evaluate() << std::endl;

Compute the complement of a set of entities

Denote \$\mathcal{E}\$ a set of entities, eg. the set of all faces (both internal and boundary faces). Denote \$\mathcal{E}_\Gamma\$ a set of entities marked by \$\Gamma\$. We wish to build \${\Gamma}^c=\mathcal{E}\backslash\Gamma\$. To compute the complement, Feel++ provides a complement template function that requires \$\mathcal{E}\$ and a predicate that return true if an entity of \$\mathcal{E}\$ belongs to \$\Gamma\$, false otherwise. The function returns mesh iterators over \$\Gamma^c\$.

#include <feel/feelmesh/complement.hpp>
...
auto E = faces(mesh);
// build set of boundary faces, equivalent to boundaryfaces(mesh)
auto bdyfaces = complement(E,[](auto const& e){return e.isOnBoundary()});
cout << "measure of bdyfaces = " << integrate(_range=bdyfaces, _expr=cst(1.)).evaluate() << std::endl;
// should be the same as above
cout << "measure of boundaryfaces = " << integrate(_range=boundaryfaces(mesh), _expr=cst(1.)).evaluate() << std::endl;

Helper function on entities set

Feel++ provides some helper functions to apply on set of entities. We denote by range_t the type of the entities set.

Table 3. Utility functions

Type

Function

Description

size_type

nelements(range_t,bool)

returns the local number of elements in entities set range_t of bool is false, other the global number which requires communication (default: global number)

WorldComm

worldComm(range_t)

returns the WorldComm associated to the entities set

Create a new range

A range can be also build directly by the user. This customized range is stored in a std container which contains the c++ references of entity object. We use boost::reference_wrapper for take c++ references and avoid copy of mesh data. All entities enumerated in the range must have same type (elements,faces,edges,points). Below we have an example which select all active elements in mesh for the current partition (i.e. identical to elements(mesh)).

auto mesh = ...;
// define reference entity type
typedef boost::reference_wrapper<typename mesh_type::element_type const> element_ref_type;
// store entities in a vector
typedef std::vector<element_ref_type> cont_range_type;
boost::shared_ptr<cont_range_type> myelts( new cont_range_type );
for (auto const& elt : elements(mesh) )
{
    myelts->push_back(boost::cref(elt));
}
// generate a range object usable in feel++
auto myrange = boost::make_tuple( mpl::size_t<MESH_ELEMENTS>(),
                                  myelts->begin(),myelts->end(),myelts );

Next, this range can be used in feel++ language.

double eval = integrate(_range=myrange,_expr=cst(1.)).evaluate()(0,0);

Mesh Markers

Elements and their associated sub-entities can be marked.

A marker is an integer specifying for example a material id, a boundary condition id or some other property associated with the entity.

A dictionary can map string to marker ids.

The dictionary is stored in the Mesh data structures and provides the set of correspondances between strings and ids.

To access a marker, it is necessary to verify that it exists as follows

for( auto const& ewrap : elements(mesh))
{
  auto const& e = unwrap_ref( ewrap );
  if ( e.hasMarker() ) (1)
  {
    std::cout << "Element " << e.id() << " has marker " << e.marker() << std::endl;
  }
  if ( e.hasMarker(5) ) (2)
  {
    std::cout << "Element " << e.id() << " has marker 5 " << e.marker(5) << std::endl;
  }
}
1 check if marker 1 (the default marker) exists, if yes then print it
2 check if marker 5 exists, if yes then print it

Mesh Operations

straightenMesh

One of the optimisations that allows to have a huge gain in computational effort is to straighten all the high order elements except for the boundary faces of the computational mesh. This is achieved by moving all the nodes associated to the high order transformation to the position these nodes would have if a first order geometrical transformation were applied. This procedure can be formalized in the following operator

\$\mathbf{\eta}^{\mathrm{straightening}}_K(\mathbf{\varphi}^N_{K}(\mathbf{x}^*)) = \left(\mathbf{\varphi}^1_{K}(\mathbf{x}^*)-\mathbf{\varphi}^N_{K}(\mathbf{x}^*)\right) - {\left( \mathbf{\varphi}^1_{K \cap \Gamma}(\mathbf{x}^*)-\mathbf{\varphi}^N_{K \cap \Gamma}(\mathbf{x}^*)\right)}\$

where \$\mathbf{x}^*\$ is any point in \$K^*\$ and \$\mathbf{\varphi}^1_{K}(\mathbf{x}^*)\$ and \$\mathbf{\varphi}^N_{K}(\mathbf{x}^*)\$ its images by the geometrical transformation of order one and order \$N\$, respectively. On one hand, the first two terms ensure that for all \$K\$ not intersecting \$\Gamma\$, the order one and \$N\$ transformations produce the same image. On the other hand, the last two terms are 0 unless the image of \$\mathbf{x}^*\$ in on \$\Gamma\$ and, in this case, we don’t move the high order image of \$\mathbf{x}^*\$. This allows to have straight internal elements and elements touching the boundary to remain high order. When applying numerical integration, specific quadratures are considered when dealing with internal elements or elements sharing a face with the boundary. The performances, thanks to this transformation, are similar to the ones obtained with first order meshes. However, it needs to be used with care as it can generate folded meshes.

createSubmesh

In multiphysics applications or using advanced numerical methods e.g. involving Lagrange multipliers, it is often required to define Function Spaces on different meshes. Theses meshes are often related. Consider for example a heat transfer problem on a domain \$\Omega\$ coupled with fluid flow problem on a domain \$\Omega_f \subset \Omega\$ as in the {uri-benchmarks}heat[Heat Transfer benchmarks]. createSubmesh allows to extract \$\Omega_f\$ out of \$\Omega\$ while keeping information on the relation between the two meshes to be able to transfer data between these meshes very efficiently.

auto mesh=loadMesh(_mesh=new Mesh<Simplex<d>>); (1)
auto fluid_mesh = createSubmesh( mesh, markedelements(mesh,"Air") ); (2)
auto face_mesh = createSubmesh( mesh, faces(mesh) ); (3)
1 create a mesh of simplices of dimension \$d\$
2 extract a mesh subregion \$\Omega_f\$ marked Air from a mesh \$\Omega\$
3 extract a \$d-1\$ mesh made of all the faces of the \$d\$ mesh