Mesh Optimisation

1. Work around the Mesh

Previously, we use the method createGMSH ( and so GMSH ) to build the mesh associated to our surface reconstrution problem. As shown ahead, it take a lot of time ( \(\simeq 70 s\)). Two possibilties will be explored here :

First, a field, named update, of the createGMSH method can be used to determine how the mesh will be built.
After this, we want to generate manually our mesh, without createGMSH.

We will also add the condition on the heigth mean, for the moment equals to 0, in order to compare directly our results to the real surface data.

1.1. Mesh Enhancement

1.1.1. Update field

The update field on the createGMSH method allow us to determine what we want the mesh built and what we want it let behind. In our case, the only useful information is the adjacency between nodes and elements, in order to parallelize in the future the code. Things, like informations on elements face, aren’t needed here and can be shelve to gain time.

This field is used like this

_uh = createGMSHMesh( _mesh=new Mesh<Hypercube<2>>,
                      _structured=1,
                      _h=doption(_name="pixelsize"),
                      _desc=domain(_name="polymere",
                                   _xmax=(x.cols()-1)*doption(_name="pixelsize"),
                                   _ymax=(x.rows()-1)*doption(_name="pixelsize")),
                      _update=size_type(MESH_UPDATE_ELEMENTS_ADJACENCY|MESH_NO_UPDATE_MEASURES));

After recompilation and execution, we finally obtain the mesh built time

[Mesh built] Time : 41.2557s

We have gain \(\simeq 30 s\) with this mesh beside the initial one.

1.1.2. Manually mesh built

An another way to gain time on mesh built is to manually build it, without using GMSH. We can so direclty build the wanted informations whitout have to read it from a gmsh file. To manage this, we inspire us from methods used to read nodes and elements from GMSH to use these informations into Feel++.

It used like this in the main code :

auto mesh = boost::make_shared<MeshStructured>( x.rows()-1, x.cols()-1,
        doption(_name="pixelsize"),Environment::worldComm() );

mesh->components().reset();
mesh->components().set( size_type(MESH_NO_UPDATE_MEASURES) );
mesh->updateForUse();

Nodes

The main part of the nodes are define like this

int start = (partId == 0) ? 0 : cx[partId]+1;
int stop = (partId == (procSize-1)) ? M_ny : cx[partId+1];
for( int j = start ; j < stop; ++j )
{
    std::cout << "j = " << j << std::endl;
    ghosts.clear();
    if (j == cx[partId] && j != 0)
    ghosts.push_back(partId-1);
    else if (j == cx[partId+1]-1 && j != M_ny-1)
    ghosts.push_back(partId+1);

    for( int i = 0; i < M_nx; ++i )
 {
        // point
        int ptid = (M_ny)*i+j;
        //std::cout << "ptId = " << ptid << std::endl;
        coords[0]=M_pixelsize*j;
        coords[1]=M_pixelsize*i;
        point_type pt( ptid,
                       coords,
                       i == 0 || i == M_nx-1
                       || j%((M_ny-1)/procSize) == 0 );

        //// Is on boundary ?
        pt.setMarker(1);
        //// Actual rank ID
        pt.setProcessId( partId );
        //// NO GHOST POINTS I SAID
        pt.setProcessIdInPartition( partId );
        pt.setNeighborPartitionIds( ghosts );
        this->addPoint( pt );
    }
}

Nodes are composed of several elements we will describe here :

An id, named ptid
A pair of coordinates, coords, which refer to its space position
A ProcessId equals to the id of the processor the node belongs
A ProcessIdInPartition equals to the id of the processor we currently work on

We choose here to define nodes present at internal boundaries, i.e at interfaces between processors, apart, in case of specific actions we could realise on them.

For left interfaces

// First column
if(partId > 0)
{
    int j = cx[partId];
    std::cout << "j = " << j << std::endl;
    ghosts.clear();
    ghosts.push_back(partId-1);
    for( int i = 0; i < M_nx; ++i )
    {
        // point
        int ptid = (M_ny)*i+j;
        //std::cout << "ptId = " << ptid << std::endl;
        coords[0]=M_pixelsize*j;
        coords[1]=M_pixelsize*i;
        point_type pt( ptid,
                       coords,
                       i == 0 || i == M_nx-1
                       || j%((M_ny-1)/procSize) == 0 );
                       // Is on boundary ?

        pt.setMarker(1);
        // Actual rank ID
        pt.setProcessId( -1 );
        pt.setProcessIdInPartition( partId );
        pt.setNeighborPartitionIds( ghosts );
        this->addPoint( pt );
    }
}

and for the right ones

// Last columnif(partId < (procSize-1))
{
    int j = cx[partId+1];
    //std::cout << "j = " << j << std::endl;
     ghosts.clear();
    ghosts.push_back(partId+1);
    for( int i = 0; i < M_nx; ++i )
    {
        // point
        int ptid = (M_ny)*i+j;
        //std::cout << "ptId = " << ptid << std::endl;
        coords[0]=M_pixelsize*j;
        coords[1]=M_pixelsize*i;
        point_type pt( ptid,
                       coords,
                       i == 0 || i == M_nx-1
                        || j%((M_ny-1)/procSize) == 0 );
        // Is on boundary ?
        pt.setMarker(1);
        // Actual rank ID
        pt.setProcessId( partId );
        pt.setProcessIdInPartition( partId );
        pt.setNeighborPartitionIds( ghosts );
        this->addPoint( pt );
    }
}

One of the main difference is that we set NeighborPartitionIds as we know we are at the interface with an another processor and we know its id.

Elements

The elements are then defined

/*
* Generates elements whithout ghosts
*/
tic();
for( int j = cx[partId]; j < cx[partId+1]; ++j )
{
    ghosts.clear();
    if (j == cx[partId] && j != 0)
        ghosts.push_back(partId-1);
    else if (j == cx[partId+1]-1 && j != M_ny-2)
        {
            ghosts.push_back(partId+1);
        }
    for( int i = 0; i < M_nx-1; ++i )
    {
        element_type e;
        int eid = (M_ny-1)*i+j; // GMSH Id

        e.setId( eid ); // eid == gmsh id
        e.setProcessIdInPartition( partId );
        e.setProcessId( partId );
        e.setMarker(1);
        e.setNeighborPartitionIds( ghosts );


        // points definition
        int ptid[4] = { (M_ny)*(i+1)+j, // 0
                        (M_ny)*(i+1)+j+1, // 1
                        (M_ny)*i+j+1, // 2
                        (M_ny)*i+j // 3
                      };

        for( int k = 0; k < 4; ++k )
        {
            e.setPoint( k, this->point( ptid[k] ) );
        }
        e.setOnBoundary( i==0 || i==(M_nx-2)
                         || j%((M_ny-2)/procSize)==0 );
        // true -> id = global elt id
        // false -> id = local elt id
        this->addElement( e, true );
        // e.id() is modified to Feel++ Id
        __idGmshToFeel.insert( std::make_pair(eid,e.id()));
        // e.id() == Feel++ Id
    }
}

First at all we define the possible processors id that can be neighborhood to the ones we want to build. We can obtain this information with the value of j, our horizontal unknown.

After that we build the several part of an element :

An id, named eid
A ProcessId equals to the id of the processor the node belongs
A ProcessIdInPartition equals to the id of the processor we currently work on
A set of points, ptid, which refer to nodes of our element. In this case, we work with square ( hypercube in 2D ) so we have four points.

We also fill a map, __idGmshToFeel, that will make the connection between GMSH Id and Feel Id.

Timing Result

The mesh built time obtain after this is equal to

[Mesh built] Time : 11.2078s

This time is decomposed as

Timer                                                       Count    Total(s)      Max(s)      Min(s)     Mean(s)   StdDev(s)
Mesh built                                                      1    1.12e+01    1.12e+01    1.12e+01    1.12e+01    0.00e+00
Mesh::updateForUse update geomap cache                        1    7.00e-02    7.00e-02    7.00e-02    7.00e-02    0.00e+00
Mesh::updateForUse update setMesh in elements and faces       1    2.85e-01    2.85e-01    2.85e-01    2.85e-01    0.00e+00
MeshStructured Elements                                       1    6.60e+00    6.60e+00    6.60e+00    6.60e+00    0.00e+00
MeshStructured Nodes                                          1    3.52e+00    3.52e+00    3.52e+00    3.52e+00    0.00e+00

We gain \(\simeq 30 s\) from using only the update field and \(\simeq 60 s\) from the original resolution.