Uniform Mesh

1. Mesh Presentation

For this project, we work with data obtain from taken images of a surface. That’s why we choose to work with an uniform mesh composed with square elements, where an element equals one pixel of the original image.

2. Work around Mesh creation

At first, we use the the method createGMSH ( and so GMSH ) to build the mesh associated to our surface reconstrution problem. However, it takes too much time ( \simeq 70 s ). That’s why we decide to explore other methods and optimize its construction :

First, a field, named update, of the createGMSH method can be used to determine how the mesh will be built.
In a second place, we try to generate manually our mesh, without createGMSH.

An another concern here is to add the condition on the reconstruct surface heigth mean. For the moment, we choose arbitrarily a mean equals to 0, it will also allow us to compare directly our results to the real surface data.

2.1. Mesh Enhancement

2.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.

2.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

// active points
    for( int j = startColPtId[partId] ; j < startColPtId[partId]+nPtByCol[partId]; ++j )
    {
        for( int i = 0 ; i < M_nx; ++i )
        {
            int ptid = (M_ny)*i+j;
            coords[0]=M_pixelsize*j;
            coords[1]=M_pixelsize*i;
            point_type pt( ptid,coords );
            pt.setProcessId( partId );
            pt.setProcessIdInPartition( partId );
            this->addPoint( pt );
        }
    }

    // ghost points (not belong to active partition)
    std::vector<size_type> ghostPointColDesc;
    if ( partId > 0 )
        ghostPointColDesc.push_back( startColPtId[partId]-1 );
    if ( partId < (nProc-1) )
        ghostPointColDesc.push_back( startColPtId[partId]+nPtByCol[partId] );
    for ( int j : ghostPointColDesc )
    {
        for( int i = 0 ; i < M_nx; ++i )
        {
            int ptid = (M_ny)*i+j;
            coords[0]=M_pixelsize*j;
            coords[1]=M_pixelsize*i;
            point_type pt( ptid,coords );
            pt.setProcessId( invalid_rank_type_value );
            pt.setProcessIdInPartition( partId  );
            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 can see two types of nodes : the normal ones and the ones named ghost, who will be useful for parallelization.

One of the main difference is that we set ProcessId to a invalid number for the ghost one as we know they doesn’t really belong to our partition.

Elements

In the same way as nodes, the elements are defined like :

// active elements
    size_type startColPtIdForElt = startColPtId[partId];
    size_type stopColPtIdForElt = startColPtId[partId]+(nPtByCol[partId]-1);
    for( int j = startColPtIdForElt ; j < stopColPtIdForElt; ++j )
    {
        rank_type neighborProcessId = invalid_rank_type_value;
        if ( partId > 0 )
            if ( j == startColPtIdForElt )
                neighborProcessId=partId-1;
        if ( partId < (nProc-1) )
            if ( j == (stopColPtIdForElt-1) )
                neighborProcessId=partId+1;

        for( int i = 0 ; i < M_nx-1; ++i )
        {
            size_type eid = (M_ny-1)*i+j; // StructuredMesh Id
            element_type e;
            e.setProcessIdInPartition( partId );
            e.setProcessId( partId );
            size_type 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( uint16_type k = 0; k < 4; ++k )
                e.setPoint( k, this->point( ptid[k]  ) );
            if ( neighborProcessId != invalid_rank_type_value )
                e.addNeighborPartitionId( neighborProcessId );

            auto const& eltInserted = this->addElement( e, true ); // e.id() is defined by Feel++
            idStructuredMeshToFeelMesh.insert( std::make_pair(eid,eltInserted.id()));
        }
    }

    // ghost elements (touch active partition with a point)
    std::vector<std::pair<size_type,rank_type> > ghostEltColDesc;
    if ( partId > 0 )
        ghostEltColDesc.push_back( std::make_pair( startColPtId[partId]-1, partId-1) );
    if ( partId < (nProc-1) )
        ghostEltColDesc.push_back( std::make_pair( startColPtId[partId]+nPtByCol[partId]-1, partId+1) );
    for ( auto const& ghostEltCol : ghostEltColDesc )
    {
        int j = ghostEltCol.first;
        rank_type partIdGhost = ghostEltCol.second;
        for( int i = 0 ; i < M_nx-1; ++i )
        {
            size_type eid = (M_ny-1)*i+j; // StructuredMesh Id
            element_type e;
            e.setProcessIdInPartition( partId );
            e.setProcessId( partIdGhost );

            size_type 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( uint16_type k = 0; k < 4; ++k )
            {
                CHECK( this->hasPoint( ptid[k] ) ) << "mesh doesnt have this point id : " << ptid[k];
                e.setPoint( k, this->point( ptid[k]  ) );
            }

            auto const& eltInserted = this->addElement( e, true ); // e.id() is defined by Feel++

            idStructuredMeshToFeelMesh.insert( std::make_pair(eid,eltInserted.id()));
            mapGhostElt.insert( std::make_pair( eid,boost::make_tuple( idStructuredMeshToFeelMesh[eid], partIdGhost ) ) );
        }
    }

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 from 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, idStructuredMeshToFeelMesh, that will make the connection between GMSH Id and Feel Id.

Time 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.