moab::Intx2Mesh Class Reference

#include <Intx2Mesh.hpp>

Inheritance diagram for moab::Intx2Mesh:

Collaboration diagram for moab::Intx2Mesh:

Public Member Functions
	Intx2Mesh (Interface *mbimpl)
virtual	~Intx2Mesh ()
ErrorCode	intersect_meshes (EntityHandle mbs1, EntityHandle mbs2, EntityHandle &outputSet)
ErrorCode	intersect_meshes_kdtree (EntityHandle mbset1, EntityHandle mbset2, EntityHandle &outputSet)
virtual ErrorCode	computeIntersectionBetweenTgtAndSrc (EntityHandle tgt, EntityHandle src, double *P, int &nP, double &area, int markb[MAXEDGES], int markr[MAXEDGES], int &nsidesSrc, int &nsidesTgt, bool check_boxes_first=false)=0
virtual ErrorCode	findNodes (EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double *iP, int nP)=0
virtual double	setup_tgt_cell (EntityHandle tgt, int &nsTgt)=0
virtual ErrorCode	FindMaxEdgesInSet (EntityHandle eset, int &max_edges)
virtual ErrorCode	FindMaxEdges (EntityHandle set1, EntityHandle set2)
virtual ErrorCode	createTags ()
virtual ErrorCode	filterByMask (Range &cells)
ErrorCode	DetermineOrderedNeighbors (EntityHandle inputSet, int max_edges, Tag &neighTag)
void	set_error_tolerance (double eps)
void	clean ()
void	set_box_error (double berror)
ErrorCode	create_departure_mesh_2nd_alg (EntityHandle &euler_set, EntityHandle &covering_lagr_set)
ErrorCode	create_departure_mesh_3rd_alg (EntityHandle &lagr_set, EntityHandle &covering_set)
void	correct_polygon (EntityHandle *foundIds, int &nP)

Protected Attributes
Interface *	mb
EntityHandle	mbs1
EntityHandle	mbs2
Range	rs1
Range	rs2
EntityHandle	outSet
Tag	gid
Tag	TgtFlagTag
Range	TgtEdges
Tag	tgtParentTag
Tag	srcParentTag
Tag	countTag
Tag	srcNeighTag
Tag	tgtNeighTag
Tag	neighTgtEdgeTag
Tag	orgSendProcTag
Tag	imaskTag
	for coverage mesh, will store the original sender More...
const EntityHandle *	tgtConn
const EntityHandle *	srcConn
CartVect	tgtCoords [MAXEDGES]
CartVect	srcCoords [MAXEDGES]
double	tgtCoords2D [MAXEDGES2]
double	srcCoords2D [MAXEDGES2]
std::vector< std::vector< EntityHandle > * >	extraNodesVec
double	epsilon_1
double	epsilon_area
std::vector< double >	allBoxes
double	box_error
EntityHandle	localRoot
Range	localEnts
unsigned int	my_rank
int	max_edges_1
int	max_edges_2
int	counting

Detailed Description

Definition at line 55 of file Intx2Mesh.hpp.

Constructor & Destructor Documentation

Intx2Mesh()

moab::Intx2Mesh::Intx2Mesh

(

Interface *

mbimpl

)

Definition at line 42 of file Intx2Mesh.cpp.

    : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
      countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), imaskTag( 0 ),
      tgtConn( nullptr ), srcConn( nullptr ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ),
      my_rank( 0 )
#ifdef MOAB_HAVE_MPI
      ,
      parcomm( nullptr ), remote_cells( nullptr ), remote_cells_with_tracers( nullptr )
#endif
      ,
      max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
{
    gid = mbimpl->globalId_tag();
}

References gid, and moab::Interface::globalId_tag().

~Intx2Mesh()

moab::Intx2Mesh::~Intx2Mesh ( )

virtual

Definition at line 57 of file Intx2Mesh.cpp.

{
    // TODO Auto-generated destructor stub
#ifdef MOAB_HAVE_MPI
    if( remote_cells )
    {
        delete remote_cells;
        remote_cells = nullptr;
    }
#endif
}

Member Function Documentation

clean()

void moab::Intx2Mesh::clean

(

)

Definition at line 1000 of file Intx2Mesh.cpp.

{
    //
    int indx = 0;
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
    {
        delete extraNodesVec[indx];
    }
    // extraNodesMap.clear();
    extraNodesVec.clear();
    // also, delete some bit tags, used to mark processed tgts and srcs
    mb->tag_delete( TgtFlagTag );
    counting = 0;  // reset counting to original value
}

References moab::Range::begin(), counting, moab::Range::end(), extraNodesVec, mb, moab::Interface::tag_delete(), TgtEdges, and TgtFlagTag.

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

computeIntersectionBetweenTgtAndSrc()

virtual ErrorCode moab::Intx2Mesh::computeIntersectionBetweenTgtAndSrc ( EntityHandle  tgt, EntityHandle  src, double *  P, int &  nP, double &  area, int  markb[MAXEDGES], int  markr[MAXEDGES], int &  nsidesSrc, int &  nsidesTgt, bool  check_boxes_first = false  )

pure virtual

Implemented in moab::IntxRllCssphere, moab::Intx2MeshOnSphere, and moab::Intx2MeshInPlane.

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

correct_polygon()

void moab::Intx2Mesh::correct_polygon	(	EntityHandle *	foundIds,
		int &	nP
	)

Definition at line 1018 of file Intx2Mesh.cpp.

{
    int i = 0;
    while( i < nP )
    {
        int nextIndex = ( i + 1 ) % nP;
        if( nodes[i] == nodes[nextIndex] )
        {
#ifdef ENABLE_DEBUG
            // we need to reduce nP, and collapse nodes
            if( dbg_1 )
            {
                std::cout << " nodes duplicated in list: ";
                for( int j = 0; j < nP; j++ )
                    std::cout << nodes[j] << " ";
                std::cout << "\n";
                std::cout << " node " << nodes[i] << " at index " << i << " is duplicated" << "\n";
            }
#endif
            // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
            // then next thing does nothing
            //  (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
            for( int k = i; k < nP - 1; k++ )
                nodes[k] = nodes[k + 1];
            nP--;  // decrease the number of nodes; also, decrease i, just if we may need to check
                   // again
            i--;
        }
        i++;
    }
    return;
}

Referenced by moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), and moab::IntxRllCssphere::findNodes().

create_departure_mesh_2nd_alg()

ErrorCode moab::Intx2Mesh::create_departure_mesh_2nd_alg	(	EntityHandle &	euler_set,
		EntityHandle &	covering_lagr_set
	)

create_departure_mesh_3rd_alg()

ErrorCode moab::Intx2Mesh::create_departure_mesh_3rd_alg	(	EntityHandle &	lagr_set,
		EntityHandle &	covering_set
	)

createTags()

ErrorCode moab::Intx2Mesh::createTags ( )

virtual

Definition at line 107 of file Intx2Mesh.cpp.

{
    if( tgtParentTag ) mb->tag_delete( tgtParentTag );
    if( srcParentTag ) mb->tag_delete( srcParentTag );
    if( countTag ) mb->tag_delete( countTag );
 
    unsigned char def_data_bit = 0;  // unused by default
    // maybe the tgt tag is better to be deleted every time, and recreated;
    // or is it easy to set all values to something again? like 0?
    MB_CHK_SET_ERR( mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit ),
                    "can't get tgt flag tag" );
    // create tgt edges if they do not exist yet; so when they are looked upon, they are found
    // this is the only call that is potentially NlogN, in the whole method
    MB_CHK_SET_ERR( mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION ), "can't get adjacent tgt edges" );
 
    // now, create a map from each edge to a list of potential new nodes on a tgt edge
    // this memory has to be cleaned up
    // change it to a vector, and use the index in range of tgt edges
    int indx = 0;
    extraNodesVec.resize( TgtEdges.size() );
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
    {
        std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
        extraNodesVec[indx]             = nv;
    }
 
    int defaultInt = -1;
 
    MB_CHK_SET_ERR( mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create TargetParent tag" );
 
    MB_CHK_SET_ERR( mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create SourceParent tag" );
 
    MB_CHK_SET_ERR( mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create Counting tag" );
 
    // for each cell in set 1, determine its neigh in set 1 (could be nullptr too)
    // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
    MB_CHK_SET_ERR( DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag ),
                    "can't determine neighbors for set 1" );
    MB_CHK_SET_ERR( DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag ),
                    "can't determine neighbors for set 2" );
 
    // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
    // them each time edges were for sure created before (tgtEdges)
    std::vector< EntityHandle > zeroh( max_edges_2, 0 );
    // if we have a tag with this name, it could be of a different size, so delete it if it exists
    if( mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag ) == MB_SUCCESS ) mb->tag_delete( neighTgtEdgeTag );
    MB_CHK_SET_ERR( mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                                        MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] ),
                    "can't create target edge neighbors tag" );
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        MB_CHK_SET_ERR( mb->get_connectivity( tgtCell, tgtConn, num_nodes ), "can't get target connectivity" );
        // account for padded polygons
        while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
            num_nodes--;
 
        int i = 0;
        for( i = 0; i < num_nodes; i++ )
        {
            EntityHandle v[2] = { tgtConn[i],
                                  tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
            std::vector< EntityHandle > adj_entities;
            MB_CHK_SET_ERR( mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT ),
                            "can't get adjacencies" );
            if( !adj_entities.size() ) MB_CHK_SET_ERR( MB_FAILURE, "no adjacencies found" );  // get out , big error
            zeroh[i] = adj_entities[0];  // should be only one edge between 2 nodes
            // also, even if number of edges is less than max_edges_2, they will be ignored, even if
            // the tag is dense
        }
        // zero out the rest
        for( i = num_nodes; i < max_edges_2; i++ )
            zeroh[i] = 0;
        // now set the value of the tag
        MB_CHK_SET_ERR( mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) ),
                        "can't set edge target edge neighbors tag" );
    }
    return MB_SUCCESS;
}

References moab::Range::begin(), countTag, DetermineOrderedNeighbors(), moab::Range::end(), extraNodesVec, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::INTERSECT, max_edges_1, max_edges_2, mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_BIT, MB_TYPE_HANDLE, MB_TYPE_INTEGER, mbs1, mbs2, neighTgtEdgeTag, rs2, moab::Range::size(), srcNeighTag, srcParentTag, moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), tgtConn, TgtEdges, TgtFlagTag, tgtNeighTag, tgtParentTag, and moab::Interface::UNION.

Referenced by intersect_meshes().

DetermineOrderedNeighbors()

ErrorCode moab::Intx2Mesh::DetermineOrderedNeighbors	(	EntityHandle	inputSet,
		int	max_edges,
		Tag &	neighTag
	)

Definition at line 194 of file Intx2Mesh.cpp.

{
    Range cells;
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( inputSet, 2, cells ), "can't get cells in set" );
 
    std::vector< EntityHandle > neighbors( max_edges );
    std::vector< EntityHandle > zeroh( max_edges, 0 );
    // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
    // one on src mesh
    MB_CHK_SET_ERR( mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &zeroh[0] ),
                    "can't create neighbors tag" );
 
    for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
    {
        EntityHandle cell = *cit;
        int nnodes        = 3;
        // will get the nnodes ordered neighbors;
        // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
        const EntityHandle* conn4;
        MB_CHK_SET_ERR( mb->get_connectivity( cell, conn4, nnodes ), "can't get connectivity of a cell" );
        int nsides = nnodes;
        // account for possible padded polygons
        while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
            nsides--;
 
        for( int i = 0; i < nsides; i++ )
        {
            EntityHandle v[2];
            v[0] = conn4[i];
            v[1] = conn4[( i + 1 ) % nsides];
            // get all cells adjacent to these 2 vertices on the edge
            std::vector< EntityHandle > adjcells;
            std::vector< EntityHandle > cellsInSet;
            MB_CHK_SET_ERR( mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT ),
                            "can't get adjacency to 2 verts" );
            // now look for the cells contained in the input set;
            // the input set should be a correct mesh, not overlapping cells, and manifold
            size_t siz = adjcells.size();
            for( size_t j = 0; j < siz; j++ )
                if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
            siz = cellsInSet.size();
 
            if( siz > 2 )
            {
                std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
                          << std::endl;
                MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" );  // non-manifold
            }
            if( siz == 1 )
            {
                // it must be the border of the input mesh;
                neighbors[i] = 0;  // we are guaranteed that ids are !=0; this is marking a border
                // borders do not appear for a sphere in serial, but they do appear for
                // parallel processing anyway
                continue;
            }
            // here siz ==2, it is either the first or second
            if( cell == cellsInSet[0] )
                neighbors[i] = cellsInSet[1];
            else
                neighbors[i] = cellsInSet[0];
        }
        // fill the rest with 0
        for( int i = nsides; i < max_edges; i++ )
            neighbors[i] = 0;
        // now simply set the neighbors tag; the last few positions will not be used, but for
        // simplicity will keep them all (MAXEDGES)
        MB_CHK_SET_ERR( mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] ), "can't set neighbors tag" );
    }
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::Interface::contains_entities(), moab::Range::end(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), moab::Interface::INTERSECT, moab::Interface::list_entities(), mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by createTags().

filterByMask()

ErrorCode moab::Intx2Mesh::filterByMask ( Range &  cells )

virtual

Definition at line 981 of file Intx2Mesh.cpp.

{
    if( !imaskTag ) return MB_SUCCESS;  // nothing to do
    size_t sz = cells.size();
    std::vector< int > masks( sz );
 
    MB_CHK_SET_ERR( mb->tag_get_data( imaskTag, cells, &masks[0] ), "can't get mask tag" );
    Range cellsToRemove;
    size_t indx = 0;
    for( Range::iterator eit = cells.begin(); eit != cells.end(); ++eit, ++indx )
    {
        if( masks[indx] ) continue;
        cellsToRemove.insert( *eit );
    }
    cells = subtract( cells, cellsToRemove );
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), imaskTag, moab::Range::insert(), mb, MB_CHK_SET_ERR, MB_SUCCESS, moab::Range::size(), moab::subtract(), and moab::Interface::tag_get_data().

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

FindMaxEdges()

ErrorCode moab::Intx2Mesh::FindMaxEdges ( EntityHandle  set1, EntityHandle  set2  )

virtual

Definition at line 99 of file Intx2Mesh.cpp.

{
    MB_CHK_SET_ERR( FindMaxEdgesInSet( set1, max_edges_1 ), "can't determine max_edges in set 1" );
    MB_CHK_SET_ERR( FindMaxEdgesInSet( set2, max_edges_2 ), "can't determine max_edges in set 2" );
 
    return MB_SUCCESS;
}

References FindMaxEdgesInSet(), max_edges_1, max_edges_2, MB_CHK_SET_ERR, and MB_SUCCESS.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

FindMaxEdgesInSet()

ErrorCode moab::Intx2Mesh::FindMaxEdgesInSet ( EntityHandle  eset, int &  max_edges  )

virtual

Definition at line 69 of file Intx2Mesh.cpp.

{
    Range cells;
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( eset, 2, cells ), "can't get entities by dimension" );
 
    max_edges = 0;  // can be 0 for point clouds
    for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
    {
        EntityHandle cell = *cit;
        const EntityHandle* conn4;
        int nnodes = 3;
        MB_CHK_SET_ERR( mb->get_connectivity( cell, conn4, nnodes ), "can't get connectivity of a cell" );
        if( nnodes > max_edges ) max_edges = nnodes;
    }
    // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
    // global) but it is better to be consistent
#ifdef MOAB_HAVE_MPI
    if( parcomm )
    {
        int local_max_edges = max_edges;
        // now reduce max_edges over all processors
        int mpi_err =
            MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
        if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
    }
#endif
 
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), mb, MB_CHK_SET_ERR, and MB_SUCCESS.

Referenced by FindMaxEdges().

findNodes()

virtual ErrorCode moab::Intx2Mesh::findNodes ( EntityHandle  tgt, int  nsTgt, EntityHandle  src, int  nsSrc, double *  iP, int  nP  )

pure virtual

Implemented in moab::IntxRllCssphere, moab::Intx2MeshOnSphere, and moab::Intx2MeshInPlane.

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

intersect_meshes()

ErrorCode moab::Intx2Mesh::intersect_meshes	(	EntityHandle	mbs1,
		EntityHandle	mbs2,
		EntityHandle &	outputSet
	)

Definition at line 555 of file Intx2Mesh.cpp.

{
    mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
    mbs2   = mbset2;
    outSet = outputSet;
#ifdef VERBOSE
    std::stringstream ffs, fft;
    ffs << "source_rank0" << my_rank << ".vtk";
    MB_CHK_SET_ERR( mb->write_mesh( ffs.str().c_str(), &mbset1, 1 ), "can't write source mesh" );
    fft << "target_rank0" << my_rank << ".vtk";
    MB_CHK_SET_ERR( mb->write_mesh( fft.str().c_str(), &mbset2, 1 ), "can't write target mesh" );
 
#endif
    // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
 
    EntityHandle startSrc = 0, startTgt = 0;
 
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs1, 2, rs1 ), "can't get source entities by dimension" );
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs2, 2, rs2 ), "can't get target entities by dimension" );
 
    // filter rs1 and rs2 by mask; remove everything with 0 mask
    // get the mask tag if it exists; if not, leave it uninitialized (nullptr)
    ErrorCode rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
    if( imaskTag != nullptr && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
 
    MB_CHK_SET_ERR( filterByMask( rs1 ), "can't filter source by mask" );
    MB_CHK_SET_ERR( filterByMask( rs2 ), "can't filter target by mask" );
 
    createTags();  // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
 
    Range rs22 = rs2;  // a copy of the initial range; we will remove from it elements as we
                       // advance ; rs2 is needed for marking the polygon to the tgt parent
 
    // create the local kdd tree with source elements; will use it to search
    // more efficiently for the seeds in advancing front;
    // some of the target cells will not be covered by source cells, and they need to be eliminated
    // early from contention
    FileOptions kdOpts( "PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;" );
    AdaptiveKDTree kd( mb );
    kd.parse_options( kdOpts );
    EntityHandle tree_root = 0;
 
    // build a kd tree with the rs1 (source) cells
    MB_CHK_SET_ERR( kd.build_tree( rs1, &tree_root ), "can't build kd tree on source cells" );
#if defined( ENABLE_DEBUG )
    for( auto it = rs22.begin(); it != rs22.end(); ++it )
    {
        EntityHandle cell        = *it;
        const EntityHandle* conn = nullptr;
        int nnodes               = 0;
        rval                     = mb->get_connectivity( cell, conn, nnodes );MB_CHK_ERR( rval );
        std::cout << " cell: \t" << " ht:" << mb->id_from_handle( cell ) << " nodes: " << nnodes
                  << " gt:" << global_id_ent( mb, cell, gid ) << " stat: " << char_stat_ent( mb, cell, TgtFlagTag )
                  << "\n";
    }
#endif
 
    while( !rs22.empty() )
    {
#if defined( ENABLE_DEBUG ) || defined( VERBOSE )
        if( rs22.size() < rs2.size() )
        {
            std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
                      << " rs22.size():" << rs22.size() << "\n";
            std::stringstream ffo;
            ffo << "file0" << counting << "rank0" << my_rank << ".h5m";
            MB_CHK_SET_ERR( mb->write_mesh( ffo.str().c_str(), &outSet, 1 ), "can't write output mesh" );
            for( auto it = rs22.begin(); it != rs22.end(); ++it )
            {
                EntityHandle cell        = *it;
                const EntityHandle* conn = nullptr;
                int nnodes               = 0;
                rval                     = mb->get_connectivity( cell, conn, nnodes );MB_CHK_ERR( rval );
                std::cout << " cell: \t" << " ht:" << mb->id_from_handle( cell ) << " nodes: " << nnodes
                          << " g:" << global_id_ent( mb, cell, gid )
                          << " stat: " << char_stat_ent( mb, cell, TgtFlagTag ) << "\n";
            }
        }
#endif
        bool seedFound = false;
        Range verified_seeds;
        for( Range::reverse_iterator it = rs22.rbegin(); it != rs22.rend(); ++it )
        {
            startTgt             = *it;
            unsigned char status = 0;
            rval                 = mb->tag_get_data( TgtFlagTag, &startTgt, 1, &status );MB_CHK_ERR( rval );
            if( 1 == status )
            {
                verified_seeds.insert( startTgt );
                continue;
            }
            int found = 0;
            // find vertex positions
            const EntityHandle* conn = nullptr;
            int nnodes               = 0;
            MB_CHK_SET_ERR( mb->get_connectivity( startTgt, conn, nnodes ), "can't get target connectivity" );
            // find leaves close to those positions
            std::vector< double > positions;
            positions.resize( nnodes * 3 );
            MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &positions[0] ), "can't get target coordinates" );
            // find leaves within a distance from each vertex of target
            // in those leaves, collect all cells; we will try for an intx in there, instead of
            // looping over all rs1 cells, as before
            Range close_source_cells;
            std::vector< EntityHandle > leaves;
            for( int i = 0; i < nnodes; i++ )
            {
                leaves.clear();
                MB_CHK_SET_ERR( kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 ),
                                "can't search for leaves" );
 
                for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
                {
                    Range tmp;
                    MB_CHK_SET_ERR( mb->get_entities_by_dimension( *j, 2, tmp ), "can't get entities by dimension" );
 
                    close_source_cells.merge( tmp.begin(), tmp.end() );
                }
            }
 
            for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
            {
                startSrc    = *it2;
                double area = 0;
                // if area is > 0 , we have intersections
                double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
                //
                int nP = 0;
                int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
                int nsTgt, nsSrc;
                MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc,
                                                                     nsTgt, true ),
                                "can't compute intersection between target and source" );
                if( area > 0 )
                {
                    found     = 1;
                    seedFound = true;
                    break;  // found 2 elements that intersect; these will be the seeds
                }
            }
            if( found )
                break;
            else
            {
#ifdef ENABLE_DEBUG
                std::cout << " on rank " << my_rank << " target cell " << " ht:" << mb->id_from_handle( startTgt )
                          << " g:" << global_id_ent( mb, startTgt, gid ) << " not intx with any source\n";
#endif
                verified_seeds.insert( startTgt );
            }
        }
        rs22 = subtract( rs22, verified_seeds );
        if( !seedFound ) continue;  // continue while(!rs22.empty())
 
        std::queue< EntityHandle > srcQueue;  // these are corresponding to Ta,
        srcQueue.push( startSrc );
        std::queue< EntityHandle > tgtQueue;
        tgtQueue.push( startTgt );
 
        unsigned char used = 1;
        // mark the start tgt quad as used, so it will not come back again
        MB_CHK_SET_ERR( mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used ), "can't set target flag" );
        while( !tgtQueue.empty() )
        {
            // flags for the side : 0 means a src cell not found on side
            // a paired src not found yet for the neighbors of tgt
            Range nextSrc[MAXEDGES];  // there are new ranges of possible next src cells for
                                      // seeding the side j of tgt cell
 
            EntityHandle currentTgt = tgtQueue.front();
            tgtQueue.pop();
            int nsidesTgt;                                                   // will be initialized now
            double areaTgtCell   = setup_tgt_cell( currentTgt, nsidesTgt );  // this is the area in the gnomonic plane
            double recoveredArea = 0;
            // get the neighbors of tgt, and if they are solved already, do not bother with that
            // side of tgt
            EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
            MB_CHK_SET_ERR( mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors ),
                            "can't get target neighbors" );
#ifdef ENABLE_DEBUG
            if( dbg_1 )
            {
                std::cout << "Next: neighbors for current tgt nsidesTgt: " << nsidesTgt << " ";
                for( int kk = 0; kk < nsidesTgt; kk++ )
                {
                    if( tgtNeighbors[kk] > 0 )
                        std::cout << " ht:" << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
                    else
                        std::cout << 0 << " ";
                }
                std::cout << std::endl;
            }
#endif
            // now get the status of neighbors; if already solved, make them 0, so not to bother
            // anymore on that side of tgt
            for( int j = 0; j < nsidesTgt; j++ )
            {
                EntityHandle tgtNeigh = tgtNeighbors[j];
                unsigned char status  = 1;
                if( tgtNeigh == 0 ) continue;
                MB_CHK_SET_ERR( mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status ),
                                "can't get target flag" );  // status 0 is unused
                if( 1 == status ) tgtNeighbors[j] = 0;      // so will not look anymore on this side of tgt
            }
 
            EntityHandle currentSrc = srcQueue.front();
            // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
            // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
            // pairs;
            //  at every moment, the queue contains pairs of cells that intersect, and they form the
            //  "advancing front"
            srcQueue.pop();
 
            Range localSrc;
            Range localSrcAlreadyTested;
            localSrc.insert( currentSrc );
#ifdef VERBOSE
            int countingStart = counting;
#endif
            // will advance-front search in the neighborhood of tgt cell, until we finish processing
            // all
            //   possible src cells; localSrc set will contain all possible src cells that cover
            //   the current tgt cell
            while( !localSrc.empty() )
            {
                //
                EntityHandle srcT = localSrc.pop_front();  // also remove from local range
                double P[10 * MAXEDGES], area;             //
                int nP           = 0;
                int nb[MAXEDGES] = { 0 };
                int nr[MAXEDGES] = { 0 };
 
                int nsidesSrc;
                // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
                // intersection points could include the vertices of initial elements
                // nb [j] = 0 means no intersection on the side j for element src (markers)
                // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
                // tgt poly.  A potential next poly in the tgt queue is the tgt poly that is
                // adjacent to this side
                MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr,
                                                                     nsidesSrc, nsidesTgt ),
                                "can't compute intersection between target and source" );
                localSrcAlreadyTested.insert( srcT );
 
                if( nP > 0 )
                {
#ifdef ENABLE_DEBUG
                    std::cout << " srcT: " << " hs:" << mb->id_from_handle( srcT )
                              << " g:" << global_id_ent( mb, srcT, gid ) << " nsidesSrc:" << nsidesSrc << "\n";
                    std::cout << " currentTgt: " << " ht:" << mb->id_from_handle( currentTgt )
                              << " g:" << global_id_ent( mb, currentTgt, gid )
                              << " stat:" << char_stat_ent( mb, currentTgt, TgtFlagTag ) << " nsidesTgt:" << nsidesTgt
                              << "\n";
                    unsigned char status = 1;
                    rval                 = mb->tag_set_data( TgtFlagTag, &currentTgt, 1, &status );MB_CHK_ERR( rval );
                    if( dbg_1 )
                    {
                        for( int k = 0; k < nsidesSrc; k++ )
                            std::cout << " nb[" << k << "]=" << nb[k];
                        std::cout << "\n";
                        for( int k = 0; k < nsidesTgt; k++ )
                            std::cout << " nr[" << k << "]=" << nr[k];
                        std::cout << "\n";
                    }
#endif
 
                    // intersection found: output P and original triangles if nP > 2
                    EntityHandle neighbors[MAXEDGES] = { 0 };
 
                    MB_CHK_SET_ERR( mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors ),
                                    "failed to get the neighbors for source element " << mb->id_from_handle( srcT ) );
 
                    Range newPotentialSrc;
                    // add neighbors to the localSrc queue, if they are not marked
                    for( int nn = 0; nn < nsidesSrc; nn++ )
                    {
                        EntityHandle neighbor = neighbors[nn];
                        if( nb[nn] > 0 )  // advance across src boundary nn
                        {
                            if( neighbor > 0 )
                            {
                                if( localSrcAlreadyTested.index( neighbor ) < 0 )  // -1
                                {
                                    localSrc.insert( neighbor );
#ifdef ENABLE_DEBUG
                                    std::cout << " local src elem " << " hs:" << mb->id_from_handle( neighbor )
                                              << " for tgt:" << "ht:" << mb->id_from_handle( currentTgt ) << "\n";
#endif
                                }
                            }
                            else  // if it is on the boundary it is a special case, maybe we need to advance more on that side,
                            // because the boundary is non-convex
                            {
                                // find the ends of edge nn, and add adjacent sources to those vertices (if they are in rs1)
                                int NumNodesSrc = 0;
                                const EntityHandle* connS;
                                rval = mb->get_connectivity( srcT, connS, NumNodesSrc );MB_CHK_SET_ERR( rval, "can't get connectivity" );
                                EntityHandle v1 = connS[nn];
                                EntityHandle v2 = connS[( nn + 1 ) % NumNodesSrc];
                                Range adjacentSourceCells1, adjacentSourceCells2;
                                rval = mb->get_adjacencies( &v1, 1, 2, false, adjacentSourceCells1 );MB_CHK_SET_ERR( rval, "can't get adjacent cells" );
                                adjacentSourceCells1 = intersect( adjacentSourceCells1, rs1 );
                                rval                 = mb->get_adjacencies( &v2, 1, 2, false, adjacentSourceCells2 );MB_CHK_SET_ERR( rval, "can't get adjacent cells" );
                                adjacentSourceCells2 = intersect( adjacentSourceCells2, rs1 );
                                adjacentSourceCells1.merge( adjacentSourceCells2 );
                                Range potentialSrc = subtract( adjacentSourceCells1, localSrcAlreadyTested );
                                newPotentialSrc.merge( potentialSrc );
                            }
                        }
                    }
                    // these might come from non-convex boundary
                    if( !newPotentialSrc.empty() )
                    {
                        localSrc.merge( newPotentialSrc );
                    }
                    // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
                    for( int nn = 0; nn < nsidesTgt; nn++ )
                    {
                        if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
                            nextSrc[nn].insert( srcT );  // potential src cell that can intersect
                                                         // the tgt neighbor nn
                    }
                    if( nP > 1 )
                    {  // this will also construct triangles/polygons in the new mesh, if needed
                        MB_CHK_SET_ERR( findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP ),
                                        "can't find nodes" );
#ifdef ENABLE_DEBUG
                        std::cout << " intersect: " << " ht:" << mb->id_from_handle( currentTgt ) << " "
                                  << " hs:" << mb->id_from_handle( srcT )
                                  << " g:" << global_id_ent( mb, currentTgt, gid ) << " "
                                  << " g:" << global_id_ent( mb, srcT, gid ) << " counting: " << counting << "\n";
#endif
                    }
 
                    recoveredArea += area;
                }
#ifdef ENABLE_DEBUG
                else if( dbg_1 )
                {
                    std::cout << " tgt, src, do not intersect: " << "ht:" << mb->id_from_handle( currentTgt ) << " "
                              << "hs:" << mb->id_from_handle( srcT ) << "\n";
                }
#endif
            }  // end while (!localSrc.empty())
            recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fraction
#if defined( ENABLE_DEBUG ) || defined( VERBOSE )
            if( fabs( recoveredArea ) > epsilon_1 )
            {
#ifdef VERBOSE
                std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
                          << " fraction error recovery:" << recoveredArea
                          << " tgtID: " << "ht:" << mb->id_from_handle( currentTgt )
                          << " countingStart:" << countingStart << "\n";
#endif
            }
#endif
            // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
            // mesh)
            rs22.erase( currentTgt );
#ifdef ENABLE_DEBUG
            std::cout << " remove cell: " << "ht:" << mb->id_from_handle( currentTgt )
                      << " g:" << global_id_ent( mb, currentTgt, gid ) << "  rs22.size():" << rs22.size() << "\n";
#endif
            // also, look at its neighbors, and add to the seeds a next one
 
            //int savedGnoPlane = plane;
            for( int j = 0; j < nsidesTgt; j++ )
            {
                EntityHandle tgtNeigh = tgtNeighbors[j];
                if( tgtNeigh == 0 || nextSrc[j].size() == 0 )  // if tgt is bigger than src, there could be no src
                                                               // to advance on that side
                    continue;
                int nsidesTgt2 = 0;
                // this also changes gnomonic plane //plane//
                setup_tgt_cell( tgtNeigh, nsidesTgt2 );  // find possible intersection with src cell from nextSrc
                for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
                {
                    EntityHandle nextB = *nit;
                    // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
                    // src quad
                    double P[10 * MAXEDGES], area;  //
                    int nP           = 0;
                    int nb[MAXEDGES] = { 0 };
                    int nr[MAXEDGES] = { 0 };
 
                    int nsidesSrc;  ///
                    MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc(
                                        /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 ),
                                    "can't compute intersection between target and source" );
                    if( area > 0 )
                    {
                        unsigned char is_used = 0;
                        rval                  = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &is_used );MB_CHK_ERR( rval );
                        if( 0 == is_used )
                        {
                            tgtQueue.push( tgtNeigh );
                            srcQueue.push( nextB );
#ifdef ENABLE_DEBUG
                            if( dbg_1 )
                                std::cout << "new polys pushed: src, tgt:" << " ht:" << mb->id_from_handle( tgtNeigh )
                                          << " hs:" << mb->id_from_handle( nextB ) << " counting: " << counting
                                          << std::endl;
#endif
                            MB_CHK_SET_ERR( mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used ),
                                            "can't set target flag" );
                        }
                        break;  // so we are done with this side of tgt, we have found a proper next
                                // seed
                    }
                }
            }
 
        }  // end while (!tgtQueue.empty())
    }
 
    // before cleaning up , we need to settle the position of the intersection points
    // on the boundary edges
    // this needs to be collective, so we should maybe wait something
#ifdef MOAB_HAVE_MPI
    MB_CHK_SET_ERR( resolve_intersection_sharing(), "can't resolve intersection sharing" );
#endif
 
    this->clean();
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), counting, createTags(), moab::AdaptiveKDTree::distance_search(), moab::Range::empty(), moab::Range::end(), epsilon_1, moab::Range::erase(), ErrorCode, filterByMask(), findNodes(), moab::Range::front(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), gid, moab::Interface::id_from_handle(), imaskTag, moab::Range::index(), moab::Range::insert(), moab::intersect(), MAXEDGES, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, mbs1, mbs2, moab::Range::merge(), my_rank, outSet, moab::AdaptiveKDTree::parse_options(), moab::Range::pop_front(), moab::Range::rbegin(), moab::Range::rend(), rs1, rs2, setup_tgt_cell(), moab::Range::size(), srcNeighTag, moab::subtract(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), TgtFlagTag, tgtNeighTag, and moab::Interface::write_mesh().

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().

intersect_meshes_kdtree()

ErrorCode moab::Intx2Mesh::intersect_meshes_kdtree	(	EntityHandle	mbset1,
		EntityHandle	mbset2,
		EntityHandle &	outputSet
	)

Slow KD-tree-based mesh intersection routine (no advancing-front).

Overview:

• Builds a KD-tree over source (mbs1) faces and, for each target (mbs2) face, queries nearby source leaves using a distance-based search around target vertices.
• For the candidate source faces gathered from nearby KD-tree leaves, computes exact polygonal intersections in a gnomonic plane, accumulates overlap area, and creates intersection polygons/nodes in outSet via findNodes.
• This path is intentionally simpler and potentially more expensive than the advancing-front algorithm used by intersect_meshes.

Inputs/Assumptions:

• mbset1 (source) fully covers mbset2 (target) on the sphere.
• Both sets contain 2D elements (triangles, quads, or generic convex polygons).
• On-sphere intersection math is performed using gnomonic projection; tolerances are derived from maximum edge lengths on the source mesh.

High-level Steps: 1) Cache 2D entities of source (rs1) and target (rs2). Optionally filter by GRID_IMASK tag to exclude masked-out elements. 2) Precompute and tag target-edge adjacency (__tgtEdgeNeighbors) for quick access when locating/creating intersection points on target boundaries. 3) Estimate tolerances: compute maximum source edge length to derive KD-tree search tolerance and box overlap epsilon; reduce across ranks under MPI. 4) Build an AdaptiveKDTree on the source faces with spherical options (PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;). 5) For each target face:

• Gather its vertex coordinates; compute an average edge length av_len.
• For each target vertex, perform kd.distance_search within radius av_len to collect nearby KD-tree leaves; accumulate their contained 2D source faces into close_source_cells.
• For each candidate source face in that range, call computeIntersectionBetweenTgtAndSrc to compute polygon intersection points and area; if area > 0, call findNodes to create nodes/polygons in outSet.
• Track recovered area vs. the target cell area (diagnostic). 6) Under MPI, reconcile shared intersection points across process boundaries via resolve_intersection_sharing. 7) Cleanup transient state and return.

Complexity Notes:

• Building the KD-tree is roughly O(N log N). For each target face, the search radius heuristic (av_len) aims to limit candidates; worst-case behavior can still approach quadratic if meshes overlap densely.

Key Data/Tags:

• tgtParentTag, srcParentTag, countTag maintain provenance and counters for created intersection entities; they are (re)created in this routine.
• __tgtEdgeNeighbors stores per-target-face edge handles to speed boundary ops.

Error handling:

• Uses MB_CHK_ERR/MB_CHK_SET_ERR macros for MOAB ErrorCode propagation.
• Cleans up tags and temporary state before returning on success.

Definition at line 322 of file Intx2Mesh.cpp.

{
    ErrorCode rval;
    mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
    mbs2   = mbset2;
    outSet = outputSet;
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs1, 2, rs1 ), "can't get source entities by dimension" );
    MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs2, 2, rs2 ), "can't get target entities by dimension" );
    // from create tags, copy relevant ones
    if( tgtParentTag ) mb->tag_delete( tgtParentTag );
    if( srcParentTag ) mb->tag_delete( srcParentTag );
    if( countTag ) mb->tag_delete( countTag );
 
    // filter rs1 and rs2 by mask; remove everything with 0 mask
    // get the mask tag if it exists; if not, leave it uninitialized (nullptr)
    rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
    if( imaskTag != nullptr && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
    MB_CHK_SET_ERR( filterByMask( rs1 ), "can't filter source by mask" );
    MB_CHK_SET_ERR( filterByMask( rs2 ), "can't filter target by mask" );
    // create tgt edges if they do not exist yet; so when they are looked upon, they are found
    // this is the only call that is potentially NlogN, in the whole method
    MB_CHK_SET_ERR( mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION ),
                    "can't get adjacent target edges" );
 
    int index = 0;
    extraNodesVec.resize( TgtEdges.size() );
    for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, index++ )
    {
        std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
        extraNodesVec[index]            = nv;
    }
 
    int defaultInt = -1;
    // Now let us create the association tags to source and target parent, along with internal counters
    MB_CHK_SET_ERR( mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create target parent tag" );
    MB_CHK_SET_ERR( mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create source parent tag" );
    MB_CHK_SET_ERR( mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                        &defaultInt ),
                    "can't create Counting tag" );
 
    // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
    // them each time edges were for sure created before (tgtEdges)
    // if we have a tag with this name, it could be of a different size, so delete it
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
    if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
    std::vector< EntityHandle > zeroh( max_edges_2, 0 );
    MB_CHK_SET_ERR( mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                                        MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] ),
                    "can't create tgt edge neighbors tag" );
 
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        MB_CHK_SET_ERR( mb->get_connectivity( tgtCell, tgtConn, num_nodes ), "can't get target connectivity" );
        // account for padded polygons
        while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
            num_nodes--;
 
        for( int i = 0; i < num_nodes; i++ )
        {
            EntityHandle v[2] = { tgtConn[i],
                                  tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
            std::vector< EntityHandle > adj_entities;
            rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
            if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval;  // get out , big error
            zeroh[i] = adj_entities[0];                                       // should be only one edge between 2 nodes
            // also, even if number of edges is less than max_edges_2, they will be ignored, even if
            // the tag is dense
        }
        // now set the value of the tag
        MB_CHK_SET_ERR( mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) ),
                        "can't set edge target edge neighbors tag" );
    }
 
    // find out max edge on source mesh;
    double max_length = 0;
    {
        std::vector< double > coords( 3 * max_edges_1, 0.0 );
        for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
        {
            const EntityHandle* conn = nullptr;
            int nnodes;
            MB_CHK_SET_ERR( mb->get_connectivity( *it, conn, nnodes ), "can't get source connectivity" );
            while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
                nnodes--;
            MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &coords[0] ), "can't get source coordinates" );
            for( int j = 0; j < nnodes; j++ )
            {
                int next = ( j + 1 ) % nnodes;
                double edge_length =
                    ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
                    ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
                    ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
                if( edge_length > max_length ) max_length = edge_length;
            }
        }
        max_length = std::sqrt( max_length );
    }
 
    // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
    double tolerance = 1.e-15;
    if( max_length < 1. )
    {
        // basically, the sag for an arc of length max_length on a circle of radius 1
        tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
        if( box_error < tolerance ) box_error = tolerance;
        tolerance = 3 * tolerance;  // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
        // be more generous, use 1.5 ~= sqrt(2.)
 
        if( !my_rank )
        {
            std::cout << " max edge length: " << max_length << "  tolerance for kd tree: " << tolerance << "\n";
            std::cout << " box overlap tolerance: " << box_error << "\n";
        }
    }
#ifdef MOAB_HAVE_MPI
    // reduce box tolerance on every task, if needed
    double min_box_eps = box_error;
    if( nullptr != parcomm ) MPI_Allreduce( &box_error, &min_box_eps, 1, MPI_DOUBLE, MPI_MIN, parcomm->comm() );
    box_error = min_box_eps;
#endif
 
    // create the kd tree on source cells, and intersect all targets in an expensive loop
    // build a kd tree with the rs1 (source) cells
    FileOptions kdOpts( "PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;" );
    AdaptiveKDTree kd( mb );
    kd.parse_options( kdOpts );
    EntityHandle tree_root = 0;
    MB_CHK_SET_ERR( kd.build_tree( rs1, &tree_root ), "can't build Kd-tree on source cells" );
 
    for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
    {
        EntityHandle tcell = *it;
        // find vertex positions
        const EntityHandle* conn = nullptr;
        int nnodes               = 0;
        MB_CHK_SET_ERR( mb->get_connectivity( tcell, conn, nnodes ), "can't get target connectivity" );
        // find leaves close to those positions
        double areaTgtCell   = setup_tgt_cell( tcell, nnodes );  // this is the area in the gnomonic plane
        double recoveredArea = 0;
        std::vector< double > positions;
        positions.resize( nnodes * 3 );
        MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &positions[0] ), "can't get target coordinates" );
 
        // distance to search will be based on average edge length
        double av_len = 0;
        for( int k = 0; k < nnodes; k++ )
        {
            int ik      = ( k + 1 ) % nnodes;
            double len1 = 0;
            for( int j = 0; j < 3; j++ )
            {
                double len2 = positions[3 * k + j] - positions[3 * ik + j];
                len1 += len2 * len2;
            }
            av_len += sqrt( len1 );
        }
        if( nnodes > 0 ) av_len /= nnodes;
        // find leaves within a distance from each vertex of target
        // in those leaves, collect all cells; we will try for an intx in there
        Range close_source_cells;
        std::vector< EntityHandle > leaves;
        for( int i = 0; i < nnodes; i++ )
        {
            leaves.clear();
            MB_CHK_SET_ERR( kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 ),
                            "can't search for leaves" );
 
            for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
            {
                Range tmp;
                MB_CHK_SET_ERR( mb->get_entities_by_dimension( *j, 2, tmp ), "can't get entities by dimension" );
 
                close_source_cells.merge( tmp.begin(), tmp.end() );
            }
        }
#ifdef VERBOSE
        if( close_source_cells.empty() )
        {
            std::cout << " there are no close source cells to target cell " << tcell
                      << " ht:" << mb->id_from_handle( tcell ) << "\n";
        }
#endif
        for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
        {
            EntityHandle startSrc = *it2;
            double area           = 0;
            // if area is > 0 , we have intersections
            double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
            //
            int nP = 0;
            int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
            int nsTgt, nsSrc;
            MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt,
                                                                 true ),
                            "can't compute intersection between target and source" );
            if( area > 0 )
            {
                if( nP > 1 )
                {  // this will also construct triangles/polygons in the new mesh, if needed
                    MB_CHK_SET_ERR( findNodes( tcell, nnodes, startSrc, nsSrc, P, nP ), "can't find nodes" );
#ifdef ENABLE_DEBUG
                    std::cout << " intersect: " << " ht:" << mb->id_from_handle( tcell ) << " "
                              << " hs:" << mb->id_from_handle( startSrc ) << " g:" << global_id_ent( mb, tcell, gid )
                              << " g:" << global_id_ent( mb, startSrc, gid ) << " counting: " << counting << "\n";
#endif
                }
                recoveredArea += area;
            }
        }
        recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fract
    }
    // before cleaning up , we need to settle the position of the intersection points
    // on the boundary edges
    // this needs to be collective, so we should maybe wait something
#ifdef MOAB_HAVE_MPI
    if( nullptr != parcomm )
    {
        MB_CHK_SET_ERR( resolve_intersection_sharing(), "can't resolve intersection sharing (correct position)" );
    }
#endif
 
    this->clean();
    return MB_SUCCESS;
}

References moab::Range::begin(), box_error, moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), counting, countTag, moab::AdaptiveKDTree::distance_search(), edge_length(), moab::Range::empty(), moab::Range::end(), epsilon_1, ErrorCode, extraNodesVec, filterByMask(), findNodes(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), gid, moab::Interface::id_from_handle(), imaskTag, moab::index, moab::Interface::INTERSECT, max_edges_1, max_edges_2, MAXEDGES, mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, MB_TYPE_INTEGER, mbs1, mbs2, moab::Range::merge(), my_rank, neighTgtEdgeTag, outSet, moab::AdaptiveKDTree::parse_options(), rs1, rs2, setup_tgt_cell(), moab::Range::size(), srcParentTag, moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), tgtConn, TgtEdges, tgtParentTag, moab::tolerance, and moab::Interface::UNION.

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().

set_box_error()

void moab::Intx2Mesh::set_box_error ( double  berror )

inline

Definition at line 138 of file Intx2Mesh.hpp.

    {
        box_error = berror;
    }

References box_error.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

set_error_tolerance()

void moab::Intx2Mesh::set_error_tolerance ( double  eps )

inline

Definition at line 121 of file Intx2Mesh.hpp.

    {
        epsilon_1    = eps;
        epsilon_area = eps * sqrt( eps );
    }

References epsilon_1, and epsilon_area.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

setup_tgt_cell()

virtual double moab::Intx2Mesh::setup_tgt_cell ( EntityHandle  tgt, int &  nsTgt  )

pure virtual

Implemented in moab::IntxRllCssphere, moab::Intx2MeshOnSphere, and moab::Intx2MeshInPlane.

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

Member Data Documentation

allBoxes

std::vector< double > moab::Intx2Mesh::allBoxes

protected

Definition at line 249 of file Intx2Mesh.hpp.

box_error

double moab::Intx2Mesh::box_error

protected

Definition at line 250 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), intersect_meshes_kdtree(), and set_box_error().

counting

int moab::Intx2Mesh::counting

protected

Definition at line 265 of file Intx2Mesh.hpp.

Referenced by clean(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), and intersect_meshes_kdtree().

countTag

Tag moab::Intx2Mesh::countTag

protected

Definition at line 216 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and intersect_meshes_kdtree().

epsilon_1

double moab::Intx2Mesh::epsilon_1

protected

Definition at line 246 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), and set_error_tolerance().

epsilon_area

double moab::Intx2Mesh::epsilon_area

protected

Definition at line 247 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), and set_error_tolerance().

extraNodesVec

std::vector< std::vector< EntityHandle >* > moab::Intx2Mesh::extraNodesVec

protected

Definition at line 244 of file Intx2Mesh.hpp.

Referenced by clean(), createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and intersect_meshes_kdtree().

gid

Tag moab::Intx2Mesh::gid

protected

Definition at line 205 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshEdges::EdgeSplits(), moab::Intx2MeshOnSphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), Intx2Mesh(), and moab::Intx2MeshOnSphere::update_tracer_data().

imaskTag

Tag moab::Intx2Mesh::imaskTag

protected

for coverage mesh, will store the original sender

Definition at line 226 of file Intx2Mesh.hpp.

Referenced by filterByMask(), intersect_meshes(), and intersect_meshes_kdtree().

localEnts

Range moab::Intx2Mesh::localEnts

protected

Definition at line 253 of file Intx2Mesh.hpp.

localRoot

EntityHandle moab::Intx2Mesh::localRoot

protected

Definition at line 252 of file Intx2Mesh.hpp.

max_edges_1

int moab::Intx2Mesh::max_edges_1

protected

Definition at line 263 of file Intx2Mesh.hpp.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), createTags(), FindMaxEdges(), and intersect_meshes_kdtree().

max_edges_2

int moab::Intx2Mesh::max_edges_2

protected

Definition at line 264 of file Intx2Mesh.hpp.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), createTags(), FindMaxEdges(), and intersect_meshes_kdtree().

mb

Interface* moab::Intx2Mesh::mb

protected

Definition at line 197 of file Intx2Mesh.hpp.

Referenced by clean(), moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), createTags(), DetermineOrderedNeighbors(), moab::Intx2MeshEdges::EdgeSplits(), filterByMask(), FindMaxEdgesInSet(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), moab::Intx2MeshEdges::orderSubEdges(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), and moab::Intx2MeshOnSphere::update_tracer_data().

mbs1

EntityHandle moab::Intx2Mesh::mbs1

protected

Definition at line 199 of file Intx2Mesh.hpp.

Referenced by createTags(), intersect_meshes(), and intersect_meshes_kdtree().

mbs2

EntityHandle moab::Intx2Mesh::mbs2

protected

Definition at line 200 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshEdges::EdgeSplits(), intersect_meshes(), and intersect_meshes_kdtree().

my_rank

unsigned int moab::Intx2Mesh::my_rank

protected

Definition at line 254 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshOnSphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), and moab::Intx2MeshOnSphere::update_tracer_data().

neighTgtEdgeTag

Tag moab::Intx2Mesh::neighTgtEdgeTag

protected

Definition at line 223 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and intersect_meshes_kdtree().

orgSendProcTag

Tag moab::Intx2Mesh::orgSendProcTag

protected

Definition at line 225 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshOnSphere::findNodes().

outSet

EntityHandle moab::Intx2Mesh::outSet

protected

Definition at line 204 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshEdges::EdgeSplits(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), and intersect_meshes_kdtree().

rs1

Range moab::Intx2Mesh::rs1

protected

Definition at line 201 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), and moab::Intx2MeshOnSphere::update_tracer_data().

rs2

Range moab::Intx2Mesh::rs2

protected

Definition at line 202 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), and moab::Intx2MeshOnSphere::update_tracer_data().

srcConn

const EntityHandle* moab::Intx2Mesh::srcConn

protected

Definition at line 231 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), and moab::IntxRllCssphere::findNodes().

srcCoords

CartVect moab::Intx2Mesh::srcCoords[MAXEDGES]

protected

Definition at line 233 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), and moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

srcCoords2D

double moab::Intx2Mesh::srcCoords2D[MAXEDGES2]

protected

Definition at line 235 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), and moab::IntxRllCssphere::findNodes().

srcNeighTag

Tag moab::Intx2Mesh::srcNeighTag

protected

Definition at line 218 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

srcParentTag

Tag moab::Intx2Mesh::srcParentTag

protected

Definition at line 215 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes_kdtree(), and moab::Intx2MeshOnSphere::update_tracer_data().

tgtConn

const EntityHandle* moab::Intx2Mesh::tgtConn

protected

Definition at line 230 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes_kdtree(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), and moab::IntxRllCssphere::setup_tgt_cell().

tgtCoords

CartVect moab::Intx2Mesh::tgtCoords[MAXEDGES]

protected

Definition at line 232 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), and moab::IntxRllCssphere::setup_tgt_cell().

tgtCoords2D

double moab::Intx2Mesh::tgtCoords2D[MAXEDGES2]

protected

Definition at line 234 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), and moab::IntxRllCssphere::setup_tgt_cell().

TgtEdges

Range moab::Intx2Mesh::TgtEdges

protected

Definition at line 210 of file Intx2Mesh.hpp.

Referenced by clean(), createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and intersect_meshes_kdtree().

TgtFlagTag

Tag moab::Intx2Mesh::TgtFlagTag

protected

Definition at line 208 of file Intx2Mesh.hpp.

Referenced by clean(), createTags(), and intersect_meshes().

tgtNeighTag

Tag moab::Intx2Mesh::tgtNeighTag

protected

Definition at line 220 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

tgtParentTag

Tag moab::Intx2Mesh::tgtParentTag

protected

Definition at line 214 of file Intx2Mesh.hpp.

Referenced by createTags(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes_kdtree(), and moab::Intx2MeshOnSphere::update_tracer_data().

---

The documentation for this class was generated from the following files:

• Intx2Mesh.hpp
• Intx2Mesh.cpp