wrap_intx.cpp

Go to the documentation of this file.

/*
 * wrap_intx.cpp
 *  Will implement the intersection method that will be callable from fortran too
 *  will be added to the sub-library IntxMesh
 *
 *
 *  Created on: Dec 14, 2013
 *      Author: iulian
 */
#include "iMesh.h"
#include "iMeshP.h"
#include "MBiMesh.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/Intx2MeshOnSphere.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/ParallelComm.hpp"
#include "MBTagConventions.hpp"
#include "moab/ParallelMergeMesh.hpp"
#include "moab/IntxUtils.hpp"

#include <sstream>
#include <mpi.h>

using namespace moab;
static double radius = 1.;
static double gtol   = 1.e-9;
static bool debug    = false;

// this mapping to coordinates will keep an index into the coords and dep_coords array
// more exactly, the fine vertices are in a Range fineVerts;
// then the vertex index i in fineVerts , fineVerts[i] will have 3 coordinates at
// index  mapping_to_coords[i] * 3
// more exactly, x = coords[ mapping_to_coords[i] * 3    ]
//               y = coords[ mapping_to_coords[i] * 3 + 1]
//               z = coords[ mapping_to_coords[i] * 3 + 2]

// in the same way, departure position for vertex depVerts[i] will be at index
//   mapping_to_coords[i] * 2 in dep_coords array (it has just latitude and longitude)
//
static int* mapping_to_coords = NULL;
static int numVertices        = 0;
// this will be instantiated at create mesh step
// should be cleaned up at the end
static Intx2MeshOnSphere* pworker = NULL;

// should get rid of this; instead of using array[NC+1][NC+1], use  row based indexing (C-style):
// parray =  new int[ (NC+1)*(NC+1) ] , and instead of array[i][j], use parray[ i*(NC+1) + j ]
#define NC 3

/*
 *  methods defined here:
 *  void update_tracer(iMesh_Instance instance,
 iBase_EntitySetHandle imesh_euler_set, int * ierr);

 void create_mesh(iMesh_Instance instance,
 iBase_EntitySetHandle * imesh_euler_set, double * coords, int * corners,
 int nc, int nelem, MPI_Fint comm, int * ierr) ;

 void intersection_at_level(iMesh_Instance instance,
 iBase_EntitySetHandle fine_set, iBase_EntitySetHandle * intx_set, double * dep_coords, double
 radius, int nc, int nelem, MPI_Fint comm, int * ierr)

 */
#ifdef __cplusplus
extern "C" {
#endif

void initialize_area_and_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, double* area_vals,
                                 int* ierr )
{

    EntityHandle eul_set = (EntityHandle)imesh_euler_set;

    moab::Interface* mb = MOABI;
    *ierr               = 1;

    ErrorCode rval;

    Range eulQuads;
    rval = mb->get_entities_by_type( eul_set, MBQUAD, eulQuads );
    ERRORV( rval, "can't get eulerian quads" );

    // area of the euler element is fixed, store it; it is used to recompute the averages at each
    // time step
    Tag tagArea = 0;
    std::string tag_name4( "Area" );
    rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get area tag" );

    std::cout << " num quads = " << eulQuads.size() << "\n";
    ERRORV( rval, "can't set area for each quad" );

    rval = mb->tag_set_data( tagArea, eulQuads, &area_vals[0] );

    *ierr = 0;
    return;
}

void update_tracer_test( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set,
                         iBase_EntitySetHandle imesh_output_set, int numTracers, double* tracer_vals, int* ierr )
{

    EntityHandle eul_set    = (EntityHandle)imesh_euler_set;
    EntityHandle output_set = (EntityHandle)imesh_output_set;

    moab::Interface* mb = MOABI;
    *ierr               = 1;

    ErrorCode rval;

    Range eulQuads;
    rval = mb->get_entities_by_type( eul_set, MBQUAD, eulQuads );
    ERRORV( rval, "can't get eulerian quads" );

    // tagElem is the average computed at each element, from nodal values
    Tag tagElem = 0;
    std::string tag_name2( "TracerAverage" );
    rval = mb->tag_get_handle( tag_name2.c_str(), numTracers, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get tracer tag " );

    // area of the euler element is fixed, store it; it is used to recompute the averages at each
    // time step
    Tag tagArea = 0;
    std::string tag_name4( "Area" );
    rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get area tag" );

    std::cout << " num quads = " << eulQuads.size() << "\n";

    rval = mb->tag_set_data( tagElem, eulQuads, &tracer_vals[0] );
    ERRORV( rval, "can't set tracer data" );

    rval = pworker->update_tracer_data( output_set, tagElem, tagArea );
    ERRORV( rval, "can't update tracer " );

    rval = mb->tag_get_data( tagElem, eulQuads, &tracer_vals[0] );
    ERRORV( rval, "can't get tracer data" );

    *ierr = 0;
    return;
}

void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, int* ierr )
{
    using namespace moab;
    const double radius = 1.;
    const double gtol   = 1.e-9;
    const bool debug    = false;

    Range ents;
    moab::Interface* mb = MOABI;
    *ierr               = 1;

    EntityHandle euler_set = (EntityHandle)imesh_euler_set;

    Intx2MeshOnSphere worker( mb );
    worker.SetRadius( radius );

    worker.SetErrorTolerance( gtol );

    EntityHandle covering_lagr_set;

    ErrorCode rval = mb->create_meshset( MESHSET_SET, covering_lagr_set );
    ERRORV( rval, "can't create covering set " );

    // we need to update the correlation tag and remote tuples
    rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );
    ERRORV( rval, "can't populate covering set " );

    if( debug )
    {
        rval = mb->write_file( "lagr.h5m", 0, 0, &covering_lagr_set, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    //
    rval = enforce_convexity( mb, covering_lagr_set );
    ERRORV( rval, "can't write covering set " );

    EntityHandle outputSet;
    rval = mb->create_meshset( MESHSET_SET, outputSet );
    ERRORV( rval, "can't create output set " );

    rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );
    ERRORV( rval, "can't intersect " );

    if( debug )
    {
        rval = mb->write_file( "output.vtk", 0, 0, &outputSet, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    // tagElem is the average computed at each element, from nodal values
    Tag tagElem = 0;
    std::string tag_name2( "TracerAverage" );
    rval = mb->tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get tracer tag " );

    // area of the euler element is fixed, store it; it is used to recompute the averages at each
    // time step
    Tag tagArea = 0;
    std::string tag_name4( "Area" );
    rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
    ERRORV( rval, "can't get area tag" );

    rval = worker.update_tracer_data( outputSet, tagElem, tagArea );
    ERRORV( rval, "can't update tracer " );

    // everything can be deleted now from intx data; polygons, etc.

    *ierr = 0;
    return;
}

ErrorCode create_coarse_mesh( Interface* mb, ParallelComm* pcomm, EntityHandle coarseSet, double* coords, int* corners,
                              int nc, int nelem, EntityHandle& start_vert, int& totalNumVertices,
                              int& numCornerVertices, std::vector< double* >& coordv )
{

    int rank = pcomm->proc_config().proc_rank();

    Tag tagid;  // global id at corners
    ErrorCode rval;
    tagid = mb->globalId_tag();

    int dum_id = -1;
    Tag partitionTag;
    mb->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, partitionTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                        &dum_id );

    // there are nelem*4 corners, and 3*(nc2+1)*(nc2+1)*nelem coordinates
    // create first a coarse mesh,
    int size_corners = 4 * nelem;
    // int size_coords=3*(nc+1)*(nc+1)*nelem;
    // order first the corners array, and eliminate duplicates
    std::vector< int > corn1( size_corners );
    std::copy( corners, corners + size_corners, corn1.begin() );
    std::sort( corn1.begin(), corn1.end() );
    corn1.erase( std::unique( corn1.begin(), corn1.end() ), corn1.end() );

    numCornerVertices = (int)corn1.size();

    // these corners are actually dofs in spectral mesh, so they will vary from 0 to
    // np*np*nel_global (at most) they are shared on edges, corners of spectral elements, so number
    // of dofs is lower than this estimate still, the number of coarse vertices is close to number
    // of cells + 2

    std::map< int, int > index_map;
    for( size_t i = 0; i < corn1.size(); i++ )
        index_map[corn1[i]] = i;

    // estimate for the vertices that will get created, from euler formula
    // v-e+f = 2 for whole sphere
    // for one task (distributed) v-e+f = 1, for one connected region
    // if multiple connex regions (out of HFSC distribution from homme), v-e+f = k, k is number of
    // connectivity regs so in any case, e = v+f -k, where k is at least 1 (1, 2, 3, etc) so in any
    // case number of coarse edges is at most e_max = v+f-1

    // this will give an estimate for the number of "fine" vertices
    int e_max = nelem + numCornerVertices - 1;
    // total number of extra vertices will be
    int numVertsOnEdges = ( nc - 1 ) * e_max;
    int numExtraVerts   = numVertsOnEdges + ( nc - 1 ) * ( nc - 1 ) * nelem;  // internal fine vertices

    totalNumVertices = numCornerVertices + numExtraVerts;  // this could be overestimated, because we are not sure
    // about the number of edges

    // used to determine the if the nodes are matching at corners of elements
    std::vector< int > index_used( numCornerVertices, 0 );

    ReadUtilIface* read_iface;
    rval = mb->query_interface( read_iface );ERRORR( rval, "can't get query intf " );

    EntityHandle start_elem, *connect;
    // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
    // over quads later
    rval = read_iface->get_node_coords( 3, totalNumVertices, 0, start_vert, coordv );ERRORR( rval, "can't get node coords " );
    // fill it up
    int stride = ( nc + 1 ) * ( nc + 1 );  // 16, for nc ==3
    // int order_in_coord[] = {0, 3, 15, 12}; // for nc=3
    /*

     *  first j, then i, so this is the order of the points in coords array, now:
     *
     *   nc(nc+1),   ...      (nc+1)*(nc+1)-1
     *
     *   2(nc+1),
     *   nc+1,   nc+2,            2*(nc+1)-1
     *   0,      1 ,    2, ..........., nc
     */
    int order_in_coord[] = { 0, nc, ( nc + 1 ) * ( nc + 1 ) - 1, nc * ( nc + 1 ) };
    for( int i = 0; i < size_corners; i++ )
    {
        int id    = corners[i];
        int index = index_map[id];

        int j         = i % 4;
        int ind_coord = i / 4;
        ind_coord     = ( ind_coord * stride + order_in_coord[j] ) * 3;
        if( index_used[index] )
        {
            if( fabs( coordv[0][index] - coords[ind_coord] ) > 0.000001 )
                std::cout << " id:" << corners[i] << " i:" << i << " j:" << j << " " << coordv[0][index] << " "
                          << coords[ind_coord] << "\n";
            continue;
        }
        coordv[0][index]     = coords[ind_coord];
        coordv[1][index]     = coords[ind_coord + 1];
        coordv[2][index]     = coords[ind_coord + 2];
        index_used[index]    = 1;
        EntityHandle vertexh = start_vert + index;
        rval                 = mb->tag_set_data( tagid, &vertexh, 1, (void*)&id );ERRORR( rval, "can't set tag id on vertex " );
    }
    // create quads; one quad for each spectral element; later, we will create fine quads
    rval = read_iface->get_element_connect( nelem, 4, MBQUAD, 0, start_elem, connect );ERRORR( rval, "can't create elements " );

    for( int i = 0; i < nelem; i++ )
    {
        for( int j = 0; j < 4; j++ )
        {
            int index_v        = index_map[corners[i * 4 + j]];
            connect[i * 4 + j] = start_vert + index_v;
        }
    }

    Range quads( start_elem, start_elem + nelem - 1 );

    mb->add_entities( coarseSet, quads );

    // create all edges adjacent to quads
    Range coarseEdges;
    rval = mb->get_adjacencies( quads, 1, true, coarseEdges, Interface::UNION );ERRORR( rval, "can't create edges " );

    mb->add_entities( coarseSet, coarseEdges );

    // see how much we overestimated the number e_max
    std::cout << " on rank " << rank << " e_max is " << e_max << " actual number of edges: " << coarseEdges.size()
              << "\n";
    rval = pcomm->resolve_shared_ents( coarseSet, 2, 1 );  // resolve vertices and edges
    ERRORR( rval, "can't resolve shared vertices and edges " );

    mb->tag_set_data( partitionTag, &coarseSet, 1, &rank );

    /* // do we have any edges?
     Range edges;
     mb->get_entities_by_dimension(0, 1, edges);
     if (rank == 0)
     std::cout << " number of edges after resolve share ents:" << edges.size()
     << "\n";*/

    rval = mb->write_file( "coarse.h5m", 0, "PARALLEL=WRITE_PART", &coarseSet, 1 );ERRORR( rval, "can't write in parallel coarse mesh" );
    // coarse mesh, from corners

    return rval;
}

// start_v and coordv refer to all vertices, including the coarse ones
ErrorCode fill_coord_on_edges( Interface* mb, std::vector< double* >& coordv, double* coords, Range& edges,
                               EntityHandle start_v, Range& coarseQuads, int nc, int numCornerVertices,
                               Tag& fineVertOnEdgeTag )
{
    ErrorCode rval = MB_SUCCESS;

    double* coordv2[3];
    for( int k = 0; k < 3; k++ )
        coordv2[k] = coordv[k] + numCornerVertices;  // they will start later

    assert( NC == nc );

    int edges_index[4][NC - 1];  // indices in the coords array for vertices, oriented positively
    /*
     *  first j, then i, so this is the order of the points in coords array, now:
     *
     *   nc(nc+1),   ...      (nc+1)*(nc+1)-1
     *
     *   2(nc+1),
     *   nc+1,   nc+2,            2*(nc+1)-1
     *   0,      1 ,    2, ..........., nc
     */
    // for first edge, oriented positive, the coords indices (*3) are 1, 2, (nc-1)
    // second edge                                                  2*(nc+1)-1, 3*(nc+1) -1, ...,
    // nc*(nc+1)-1 third edge:                                                  (nc+1) * (nc+1) -2,
    // ..., nc*(nc+1) +1 fourth edge:                                                 (nc-1)*(nc+1),
    // ..., (nc+1)
    for( int j = 1; j <= nc - 1; j++ )
    {
        edges_index[0][j - 1] = j;                                // for nc = 3: 1, 2
        edges_index[1][j - 1] = ( j + 1 ) * ( nc + 1 ) - 1;       //             7, 11
        edges_index[2][j - 1] = ( nc + 1 ) * ( nc + 1 ) - j - 1;  //             14, 13
        edges_index[3][j - 1] = ( nc - j ) * ( nc + 1 );          //             8, 4
    }
    // int num_quads=(int)coarseQuads.size();
    int stride = ( nc + 1 ) * ( nc + 1 );
    int indexv = 0;
    for( Range::iterator eit = edges.begin(); eit != edges.end(); ++eit )
    {
        EntityHandle edge = *eit;
        std::vector< EntityHandle > faces;
        rval = mb->get_adjacencies( &edge, 1, 2, false, faces );ERRORR( rval, "can't get adjacent faces." );
        if( faces.size() < 1 ) return MB_FAILURE;
        int sense = 0, side_number = -1, offset = -1;
        EntityHandle quad = faces[0];  // just consider first quad
        rval              = mb->side_number( quad, edge, side_number, sense, offset );ERRORR( rval, "can't get side number" );
        int indexq = coarseQuads.index( quad );

        if( indexq == -1 ) return MB_FAILURE;

        EntityHandle firstRefinedV = start_v + numCornerVertices + indexv;
        rval                       = mb->tag_set_data( fineVertOnEdgeTag, &edge, 1, &firstRefinedV );ERRORR( rval, "can't set refined vertex tag" );
        // copy the right coordinates from the coords array to coordv2 array

        double* start_quad = &coords[3 * stride * indexq];
        if( sense > 0 )
        {
            for( int k = 1; k <= nc - 1; k++ )
            {
                int index_in_quad  = edges_index[side_number][k - 1] * 3;
                coordv2[0][indexv] = start_quad[index_in_quad];
                coordv2[1][indexv] = start_quad[index_in_quad + 1];
                coordv2[2][indexv] = start_quad[index_in_quad + 2];
                indexv++;
            }
        }
        else
        {
            // sense < 0, so we will traverse the edge in inverse sense
            for( int k = 1; k <= nc - 1; k++ )
            {
                int index_in_quad  = edges_index[side_number][nc - 1 - k] * 3;
                coordv2[0][indexv] = start_quad[index_in_quad];
                coordv2[1][indexv] = start_quad[index_in_quad + 1];
                coordv2[2][indexv] = start_quad[index_in_quad + 2];
                indexv++;
            }
        }
    }
    return rval;
}
/*
 ErrorCode resolve_interior_verts_on_bound_edges(Interface * mb, ParallelComm * pcomm,
 Range & edges)
 {
 // edges are coarse edges;
 ErrorCode rval;
 int rank=pcomm->proc_config().proc_rank();
 Range sharedCoarseEdges=edges;// filter the non shared ones
 rval = pcomm->filter_pstatus(sharedCoarseEdges, PSTATUS_SHARED, PSTATUS_AND);ERRORR(rval, "can't filter coarse edges
 "); ParallelMergeMesh pmerge(pcomm, 0.0001); ErrorCode rval = pmerge.merge();

 return rval;
 }*/
ErrorCode create_fine_mesh( Interface* mb, ParallelComm* pcomm, EntityHandle coarseSet, EntityHandle fine_set,
                            double* coords, int nc, int nelem, EntityHandle start_vert, int numCornerVertices,
                            std::vector< double* >& coordv )
{
    int rank   = pcomm->proc_config().proc_rank();
    int stride = ( nc + 1 ) * ( nc + 1 );  // 16, for nc ==3
    // there are stride*3 coordinates for each coarse quad, representing the fine mesh

    Tag fineVertTag;
    // will store the first refined vertex on an edge in a entity handle tag
    EntityHandle def = 0;
    ErrorCode rval =
        mb->tag_get_handle( "__firstvertex", 1, MB_TYPE_HANDLE, fineVertTag, MB_TAG_DENSE | MB_TAG_CREAT, &def );ERRORR( rval, "can't get tag on first vertex" );

    Range coarseQuads;
    Range edges;
    rval = mb->get_entities_by_dimension( coarseSet, 2, coarseQuads );ERRORR( rval, "can't get coarse quads  " );

    rval = mb->get_entities_by_dimension( coarseSet, 1, edges );ERRORR( rval, "can't get coarse edges  " );

    Range verts;
    rval = mb->get_connectivity( coarseQuads, verts );ERRORR( rval, "can't get coarse vertices " );

    /*std::cout <<" local coarse mesh on rank " << rank << "  "<< coarseQuads.size() << " quads, "
     << edges.size() << " edges, " << verts.size() <<  " vertices.\n";*/

    int dum_id = -1;
    Tag partitionTag;
    mb->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, partitionTag, MB_TAG_SPARSE | MB_TAG_CREAT,
                        &dum_id );
    // fine mesh, with all coordinates
    // std::vector<double *> coordv2;

    ReadUtilIface* read_iface;
    rval = mb->query_interface( read_iface );ERRORR( rval, "can't get query intf " );

    ;
    // create verts, (nc+1)*(nc+1)*nelem - verts.size()
    //
    /* int numVertsOnEdges = (nc-1)*(int)edges.size();
     int numExtraVerts = numVertsOnEdges + (nc-1)*(nc-1)*nelem; // internal fine vertices
     rval = read_iface->get_node_coords(3, numExtraVerts, 0,
     start_vert, coordv2);*/
    // ERRORR(rval, "can't get coords fine mesh");
    // fill coordinates for vertices on the edges, then vertices in the interior of coarse quads
    // we know that all quads are in order, their index corresponds to index in coords array
    rval =
        fill_coord_on_edges( mb, coordv, coords, edges, start_vert, coarseQuads, nc, numCornerVertices, fineVertTag );ERRORR( rval, "can't fill edges vertex coords on fine mesh" );

    EntityHandle start_elem, *connect;
    rval = read_iface->get_element_connect( nelem * nc * nc, 4, MBQUAD, 0, start_elem, connect );ERRORR( rval, "can't create elements fine mesh" );

    int iv     = ( nc - 1 ) * edges.size() + numCornerVertices;  // iv is in the coordv array indices
    start_vert = start_vert + numCornerVertices + ( nc - 1 ) * edges.size();

    /* // add a child to the mesh set, with the ordered vertices, as they come out in the list of
     coordinates EntityHandle vertSet; rval = mb->create_meshset(MESHSET_ORDERED, vertSet);ERRORR(rval, "can't create
     vertex set ");

     rval = mb->add_parent_child(fine_set, vertSet);ERRORR(rval, "can't create parent child relation between fine set
     and vertSet ");*/

    std::vector< EntityHandle > vertList;
    vertList.reserve( nelem * ( nc + 1 ) * ( nc + 1 ) );  // will have a list of vertices, in order

    // now fill coordinates on interior nodes; also mark the start for each interior vertex
    //  in a coarse quad
    for( int ie = 0; ie < nelem; ie++ )
    {
        // just fill coordinates for an array of (nc-1)*(nc-1) vertices
        EntityHandle firstVert = start_vert + ( nc - 1 ) * ( nc - 1 ) * ie;
        EntityHandle eh        = coarseQuads[ie];
        rval                   = mb->tag_set_data( fineVertTag, &eh, 1, &firstVert );ERRORR( rval, "can't set refined vertex tag" );

        int index_coords = stride * ie;
        for( int j = 1; j <= ( nc - 1 ); j++ )
        {
            for( int i = 1; i <= ( nc - 1 ); i++ )
            {
                int indx2     = 3 * ( index_coords + ( nc + 1 ) * j + i );
                coordv[0][iv] = coords[indx2];
                coordv[1][iv] = coords[indx2 + 1];
                coordv[2][iv] = coords[indx2 + 2];
                iv++;
            }
        }
    }

    Range quads3( start_elem, start_elem + nelem * nc * nc - 1 );

    int ic = 0;
    for( int ie = 0; ie < nelem; ie++ )
    {
        // just fill coordinates for an array of (nc-1)*(nc-1) vertices
        EntityHandle arr2[NC + 1][NC + 1];  //
        /*
         *     (nc,0)         (nc,nc)
         *
         *
         *     (1,0)           (1,nc)
         *     (0,0) (0,1)     (0,nc)
         */
        EntityHandle coarseQ      = coarseQuads[ie];
        const EntityHandle* conn4 = NULL;
        int nnodes                = 0;
        rval                      = mb->get_connectivity( coarseQ, conn4, nnodes );ERRORR( rval, "can't get conn of coarse quad" );
        if( nnodes != 4 ) return MB_FAILURE;

        arr2[0][0]   = conn4[0];
        arr2[nc][0]  = conn4[1];
        arr2[nc][nc] = conn4[2];
        arr2[0][nc]  = conn4[3];

        // get the coarse edges
        std::vector< EntityHandle > aedges;
        rval = mb->get_adjacencies( &coarseQ, 1, 1, false, aedges );ERRORR( rval, "can't get adje edges of coarse quad" );
        assert( (int)aedges.size() == 4 );

        for( int k = 0; k < 4; k++ )
        {
            EntityHandle edh = aedges[k];
            /*
             edges_index[0][j-1] = j;                        // for nc = 3: 1, 2
             edges_index[1][j-1] = (j+1)*(nc+1) - 1;         //             7, 11
             edges_index[2][j-1] = (nc+1)*(nc+1) - j - 1;    //             14, 13
             edges_index[3][j-1] = (nc-j) * (nc+1) ;         //             8, 4
             */
            int sense = 0, side_number = -1, offset = -1;
            rval = mb->side_number( coarseQ, edh, side_number, sense, offset );ERRORR( rval, "can't get side number" );
            EntityHandle firstV;  // first vertex on edge, if edge oriented positively
            rval = mb->tag_get_data( fineVertTag, &edh, 1, &firstV );ERRORR( rval, "can't get first vertex tag on edge" );
            if( sense > 0 )
            {
                if( 0 == side_number )
                {
                    for( int i = 1; i <= nc - 1; i++ )
                        arr2[i][0] = firstV + i - 1;
                }
                else if( 1 == side_number )
                {
                    for( int j = 1; j <= nc - 1; j++ )
                        arr2[nc][j] = firstV + j - 1;
                }
                else if( 2 == side_number )
                {
                    for( int i = nc - 1; i >= 1; i-- )
                        arr2[i][nc] = firstV + nc - 1 - i;
                }
                else if( 3 == side_number )
                {
                    for( int j = nc - 1; j >= 1; j-- )
                        arr2[0][j] = firstV + nc - 1 - j;
                }
            }
            else  // if (sense<0)
            {
                if( 0 == side_number )
                {
                    for( int i = 1; i <= nc - 1; i++ )
                        arr2[i][0] = firstV + nc - 1 - i;
                }
                else if( 1 == side_number )
                {
                    for( int j = 1; j <= nc - 1; j++ )
                        arr2[nc][j] = firstV + nc - 1 - j;
                }
                else if( 2 == side_number )
                {
                    for( int i = nc - 1; i >= 1; i-- )
                        arr2[i][nc] = firstV + i - 1;
                }
                else if( 3 == side_number )
                {
                    for( int j = nc - 1; j >= 1; j-- )
                        arr2[0][j] = firstV + j - 1;
                }
            }
        }
        // fill the interior points matrix
        EntityHandle firstV;  // first vertex on interior of coarse quad
        rval = mb->tag_get_data( fineVertTag, &coarseQ, 1, &firstV );ERRORR( rval, "can't get first vertex tag on coarse tag" );
        int inc = 0;
        for( int j = 1; j <= nc - 1; j++ )
        {
            for( int i = 1; i <= nc - 1; i++ )
            {
                arr2[i][j] = firstV + inc;
                inc++;
            }
        }
        // fill now the matrix of quads; vertices are from 0 to (nc+1)*(nc+1)-1
        for( int j1 = 0; j1 < nc; j1++ )
        {
            for( int i1 = 0; i1 < nc; i1++ )
            {
                connect[ic++] = arr2[i1][j1];          // first one
                connect[ic++] = arr2[i1 + 1][j1];      //
                connect[ic++] = arr2[i1 + 1][j1 + 1];  // opp diagonal
                connect[ic++] = arr2[i1][j1 + 1];      //
            }
        }

        for( int j1 = 0; j1 <= nc; j1++ )
        {
            for( int i1 = 0; i1 <= nc; i1++ )
            {
                vertList.push_back( arr2[i1][j1] );
            }
        }
    }
    /*
     rval = mb->add_entities(vertSet, &vertList[0], (int)vertList.size());ERRORR(rval,"can't add to the vert set the
     list of ordered vertices");*/

    mb->add_entities( fine_set, quads3 );
    // notify MOAB of the new elements
    rval = read_iface->update_adjacencies( start_elem, nelem * nc * nc, 4, connect );ERRORR( rval, "can't update adjacencies on fine quads" );

    rval = mb->tag_set_data( partitionTag, &fine_set, 1, &rank );ERRORR( rval, "can't set partition tag on fine set" );

    /*
     // delete the coarse mesh, except vertices
     mb->delete_entities(coarseQuads);
     mb->delete_entities(edges);
     */

    // the vertices on the boundary edges of the partition need to be shared and resolved
    ParallelMergeMesh pmerge( pcomm, 0.0001 );
    rval = pmerge.merge();ERRORR( rval, "can't resolve vertices on interior of boundary edges " );

    rval = mb->get_connectivity( quads3, verts );ERRORR( rval, "can't get vertices " );

    Range owned_verts = verts;
    rval              = pcomm->filter_pstatus( owned_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT );ERRORR( rval, "can't filter for owned vertices only" );

    Range entities[4];
    entities[0] = owned_verts;
    entities[2] = quads3;
    // assign new ids only for owned entities
    // will eliminate gaps in global id space for vertices
    rval = pcomm->assign_global_ids( entities, 2, 1, true, false );ERRORR( rval, "can't assign global ids for vertices " );

    /*ErrorCode ParallelComm::assign_global_ids( Range entities[],
     const int dimension,
     const int start_id,
     const bool parallel,
     const bool owned_only) */

    std::stringstream fff;
    fff << "fine0" << pcomm->proc_config().proc_rank() << ".h5m";
    mb->write_mesh( fff.str().c_str(), &fine_set, 1 );

    rval = mb->write_file( "fine.h5m", 0, "PARALLEL=WRITE_PART", &fine_set, 1 );ERRORR( rval, "can't write set 3, fine " );

    // we need to keep a mapping index, from the coords array to the vertex handles
    // so, for a given vertex entity handle, at what index in the coords array the vertex
    // coordinates are?
    // in the coords array, vertices are repeated 2 times if they are interior to a coarse edge, and
    // repeated 3 or 4 times if they are a corner vertex in a coarse quad

    numVertices       = (int)verts.size();
    mapping_to_coords = new int[numVertices];
    for( int k = 0; k < numVertices; k++ )
        mapping_to_coords[k] = -1;  // it means it was not located yet in vertList
    // vertList is parallel to the coords and dep_coords array

    // now loop over vertsList, and see where
    // vertList has nelem * (nc+1)*(nc+1) vertices; loop over them, and see where are they located

    for( int kk = 0; kk < (int)vertList.size(); kk++ )
    {
        EntityHandle v = vertList[kk];
        int index      = verts.index( v );
        if( -1 == index )
        {
            std::cout << " can't locate vertex " << v << " in vertex Range \n";
            return MB_FAILURE;
        }
        if( mapping_to_coords[index] == -1 )  // it means the vertex v was not yet encountered in the vertList
        {
            mapping_to_coords[index] = kk;
        }
    }
    // check that every mapping has an index different from -1
    for( int k = 0; k < numVertices; k++ )
    {
        if( mapping_to_coords[k] == -1 )
        {
            {
                std::cout << " vertex at index " << k << " in vertex Range " << verts[k] << " is not mapped \n";
                return MB_FAILURE;  //
            }
        }
    }

    return rval;
}

// this is called from Fortran 90
void create_mesh( iMesh_Instance instance, iBase_EntitySetHandle* imesh_euler_set,
                  iBase_EntitySetHandle* imesh_departure_set, iBase_EntitySetHandle* imesh_intx_set, double* coords,
                  int* corners, int nc, int nelem, MPI_Fint comm, int* ierr )
{
    /* double * coords=(double*) icoords;
     int * corners = (int*) icorners;*/
    *ierr            = 1;
    Interface* mb    = MOABI;
    MPI_Comm mpicomm = MPI_Comm_f2c( comm );
    // instantiate parallel comm now or not?
    ParallelComm* pcomm = new ParallelComm( mb, mpicomm );

    EntityHandle coarseSet;
    ErrorCode rval = mb->create_meshset( MESHSET_SET, coarseSet );
    ERRORV( rval, "can't create coarse set " );

    EntityHandle start_vert;
    int totalNumVertices;
    int numCornerVertices;
    std::vector< double* > coordv;
    rval = create_coarse_mesh( mb, pcomm, coarseSet, coords, corners, nc, nelem, start_vert, totalNumVertices,
                               numCornerVertices, coordv );
    ERRORV( rval, "can't create coarse set " );

    EntityHandle fine_set;
    rval = mb->create_meshset( MESHSET_SET, fine_set );
    ERRORV( rval, "can't create fine set " );

    rval = create_fine_mesh( mb, pcomm, coarseSet, fine_set, coords, nc, nelem, start_vert, numCornerVertices, coordv );
    ERRORV( rval, "can't create fine mesh set " );

    *imesh_departure_set = (iBase_EntitySetHandle)fine_set;

    EntityHandle euler_set;
    rval = mb->create_meshset( MESHSET_SET, euler_set );
    ERRORV( rval, "can't create moab euler set " );
    *imesh_euler_set = (iBase_EntitySetHandle)euler_set;

    // call in Intx utils
    // it will copy the second set from the first set
    rval = deep_copy_set_with_quads( mb, fine_set, euler_set );
    ERRORV( rval, "can't populate lagrange set " );

    EntityHandle intx_set;
    rval = mb->create_meshset( MESHSET_SET, intx_set );
    ERRORV( rval, "can't create output set " );

    *imesh_intx_set = (iBase_EntitySetHandle)intx_set;

    pworker = new Intx2MeshOnSphere( mb );

    pworker->set_box_error( 100 * gtol );
    Range local_verts;
    rval = pworker->build_processor_euler_boxes( euler_set,
                                                 local_verts );  // output also the local_verts
    ERRORV( rval, "can't compute euler boxes " );
    pworker->SetErrorTolerance( gtol );

    *ierr = 0;
    return;
}
ErrorCode set_departure_points_position( Interface* mb, EntityHandle lagrSet, double* dep_coords, double radius2 )
{

    // the departure quads are created in the same order as the fine quads
    // for each coarse element, there are nc*nc fine quads
    // their vertices are in order
    Range lagr_quads;
    ErrorCode rval = mb->get_entities_by_type( lagrSet, MBQUAD, lagr_quads );ERRORR( rval, "can't get lagrange quads" );

    // get all vertices from lagr_quads
    Range lVerts;
    rval = mb->get_connectivity( lagr_quads, lVerts );ERRORR( rval, "can't get lagrangian vertices (departure)" );

    // they are parallel to the verts Array, they must have the same number of vertices
    assert( numVertices == (int)lVerts.size() );

    for( int i = 0; i < numVertices; i++ )
    {
        EntityHandle v = lVerts[i];
        int index      = mapping_to_coords[i];
        assert( -1 != index );

        SphereCoords sph;
        sph.R   = radius2;
        sph.lat = dep_coords[2 * index];
        sph.lon = dep_coords[2 * index + 1];

        CartVect depPoint = spherical_to_cart( sph );
        rval              = mb->set_coords( &v, 1, (double*)depPoint.array() );ERRORR( rval, "can't set position of vertex" );
    }

    return MB_SUCCESS;
}
void intersection_at_level( iMesh_Instance instance, iBase_EntitySetHandle fine_set, iBase_EntitySetHandle lagr_set,
                            iBase_EntitySetHandle intx_set, double* dep_coords, double radius2, int* ierr )
{
    *ierr         = 1;
    Interface* mb = MOABI;
    // MPI_Comm mpicomm = MPI_Comm_f2c(comm);
    // instantiate parallel comm now or not?

    EntityHandle lagrMeshSet = (EntityHandle)lagr_set;

    ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
    if( NULL == pcomm ) return;  // error is 1

    // set the departure tag on the fine mesh vertices
    ErrorCode rval = set_departure_points_position( mb, lagrMeshSet, dep_coords, radius2 );
    ERRORV( rval, "can't set departure tag" );
    if( debug )
    {
        std::stringstream fff;
        fff << "lagr0" << pcomm->proc_config().proc_rank() << ".vtk";
        rval = mb->write_mesh( fff.str().c_str(), &lagrMeshSet, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    // it should be done earlier
    pworker->SetRadius( radius );

    EntityHandle covering_set;
    rval = pworker->create_departure_mesh_3rd_alg( lagrMeshSet, covering_set );
    ERRORV( rval, "can't compute covering set " );

    if( debug )
    {
        std::stringstream fff;
        fff << "cover" << pcomm->proc_config().proc_rank() << ".vtk";
        rval = mb->write_mesh( fff.str().c_str(), &covering_set, 1 );

        ERRORV( rval, "can't write covering set " );
    }
    EntityHandle intxSet = (EntityHandle)intx_set;
    rval                 = pworker->intersect_meshes( covering_set, (EntityHandle)fine_set, intxSet );
    ERRORV( rval, "can't intersect " );

    if( debug )
    {
        std::stringstream fff;
        fff << "intx0" << pcomm->proc_config().proc_rank() << ".vtk";
        rval = mb->write_mesh( fff.str().c_str(), &intxSet, 1 );
        ERRORV( rval, "can't write covering set " );
    }

    return;
}
void cleanup_after_intersection( iMesh_Instance instance, iBase_EntitySetHandle fine_set,
                                 iBase_EntitySetHandle lagr_set, iBase_EntitySetHandle intx_set, int* ierr )
{
    *ierr         = 1;
    Interface* mb = MOABI;
    // delete elements
    // delete now the polygons and the elements of out_set
    // also, all verts that are not in euler set or lagr_set
    Range allVerts;
    ErrorCode rval = mb->get_entities_by_dimension( 0, 0, allVerts );
    ERRORV( rval, "can't get all vertices" );

    Range allElems;
    rval = mb->get_entities_by_dimension( 0, 2, allElems );
    ERRORV( rval, "can't get all elems" );

    Range polys;  // first euler polys
    rval = mb->get_entities_by_dimension( (EntityHandle)fine_set, 2, polys );
    ERRORV( rval, "can't get all polys from fine set" );

    // add to polys range the lagr polys
    rval = mb->get_entities_by_dimension( (EntityHandle)lagr_set, 2,
                                          polys );  // do not delete lagr set either, with its vertices
    ERRORV( rval, "can't get all polys from lagr set" );
    // add to the connecVerts range all verts, from all initial polys
    Range vertsToStay;
    rval = mb->get_connectivity( polys, vertsToStay );
    ERRORV( rval, "get verts that stay" );

    Range todeleteVerts = subtract( allVerts, vertsToStay );

    Range todeleteElem = subtract( allElems, polys );  // this is coarse mesh too (if still here)

    // empty the out mesh set
    EntityHandle out_set = (EntityHandle)intx_set;
    rval                 = mb->clear_meshset( &out_set, 1 );
    ERRORV( rval, "clear mesh set" );

    rval = mb->delete_entities( todeleteElem );
    ERRORV( rval, "delete intx elements" );
    rval = mb->delete_entities( todeleteVerts );
    ERRORV( rval, "failed to delete intx vertices" );

    *ierr = 0;
    return;
}

void cleanup_after_simulation( int* ierr )
{
    delete[] mapping_to_coords;
    numVertices = 0;
    *ierr       = 0;
    return;
}
#ifdef __cplusplus
}  // extern "C"
#endif