## 1. Using computational meshes

### 1.1. 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 include directive in modules/reference/pages/Mesh/README.adoc - include::../../../../codes/03-mymesh.cpp[]``

As always, we initialize the Feel++ 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].

### 1.2. Basic Meshes

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

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

### 1.3. Loading Meshes

#### 1.3.1. `loadMesh`

You can use this function to load mesh from different formats

 in the case of a `.geo` file, Feel++ 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. Feel++ 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 > > );``

#### 1.3.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.

#### 1.3.3. 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 );``````

### 1.4. Creating Meshes

#### 1.4.1. 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.

#### 1.4.2. 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 );``````

#### 1.4.3. 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 ) );``````

### 1.5. 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
...
}``````

#### 1.5.1. 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` `elements(mesh,entity_process_t)` all the elements of mesh associated to entity_process_t. `ext_elements_t` `markedelements(mesh, id, entity_process_t)` all the elements marked id of mesh associated to entity_process_t. `ext_faces_t` `faces(mesh,entity_process_t)` all the faces of mesh associated to entity_process_t. `ext_faces_t` `markedfaces(mesh, id, entity_process_t)` all the faces marked id of mesh associated to entity_process_t. `ext_edges_t` `edges(mesh,entity_process_t)` all the edges of mesh associated to entity_process_t. `ext_edges_t` `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` 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() )
{...}
}``````

#### 1.5.2. 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;``````

#### 1.5.3. 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;``````

#### 1.5.4. 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.

 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

#### 1.5.5. 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);``

### 1.6. 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

### 1.7. Mesh Operations

#### 1.7.1. `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.

#### 1.7.2. `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