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 28 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( NULL ), srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ),
      my_rank( 0 )
#ifdef MOAB_HAVE_MPI
      ,
      parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
#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 43 of file Intx2Mesh.cpp.

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

Member Function Documentation

clean()

void moab::Intx2Mesh::clean

(

)

Definition at line 836 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 854 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 93 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?
    ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "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
    rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "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;
 
    rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
 
    rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
 
    rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
 
    // for each cell in set 1, determine its neigh in set 1 (could be null too)
    // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
    rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
    rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "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
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
    if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                               MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
        // 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;
            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
        rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
    }
    return MB_SUCCESS;
}

References moab::Range::begin(), countTag, DetermineOrderedNeighbors(), moab::Range::end(), ErrorCode, 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 168 of file Intx2Mesh.cpp.

{
    Range cells;
    ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "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
    rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "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;
        rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "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;
            rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't 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() )
                          << "\n";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)
        rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
    }
    return MB_SUCCESS;
}

References moab::Range::begin(), moab::Interface::contains_entities(), moab::Range::end(), ErrorCode, 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 817 of file Intx2Mesh.cpp.

{
    if( !imaskTag ) return MB_SUCCESS;  // nothing to do
    size_t sz = cells.size();
    std::vector< int > masks( sz );
 
    ErrorCode rval = mb->tag_get_data( imaskTag, cells, &masks[0] );MB_CHK_ERR( rval );
    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(), ErrorCode, imaskTag, moab::Range::insert(), mb, MB_CHK_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 85 of file Intx2Mesh.cpp.

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

References ErrorCode, 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 55 of file Intx2Mesh.cpp.

{
    Range cells;
    ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );
 
    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;
        rval       = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "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(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), mb, MB_CHK_ERR, 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 455 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;
#ifdef VERBOSE
    std::stringstream ffs, fft;
    ffs << "source_rank0" << my_rank << ".vtk";
    rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
    fft << "target_rank0" << my_rank << ".vtk";
    rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
 
#endif
    // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
 
    EntityHandle startSrc = 0, startTgt = 0;
 
    rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
    rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
 
    // filter rs1 and rs2 by mask; remove everything with 0 mask
    // get the mask tag if it exists; if not, leave it uninitialized (NULL)
    mb->tag_get_handle( "GRID_IMASK", imaskTag );
    if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
 
    rval = filterByMask( rs1 );MB_CHK_ERR( rval );
    rval = filterByMask( rs2 );MB_CHK_ERR( rval );
 
    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
 
    // 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;
    rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
 
    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
                      << "\n";
            std::stringstream ffo;
            ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
            rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
        }
#endif
        bool seedFound = false;
        for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
        {
            startTgt  = *it;
            int found = 0;
            // find vertex positions
            const EntityHandle* conn = NULL;
            int nnodes               = 0;
            rval                     = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
            // find leaves close to those positions
            std::vector< double > positions;
            positions.resize( nnodes * 3 );
            rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
            // 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();
                rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
 
                for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
                {
                    Range tmp;
                    rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
 
                    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;
                rval =
                    computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
                if( area > 0 )
                {
                    found     = 1;
                    seedFound = true;
                    break;  // found 2 elements that intersect; these will be the seeds
                }
            }
            if( found )
                break;
            else
            {
#if defined( VERBOSE )
                std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
                          << " not intx with any source\n";
#endif
                rs22.erase( startTgt );
            }
        }
        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 );
 
        Range toResetSrcs;  // will be used to reset src flags for every tgt element processed
 
        /*if (my_rank==0)
          dbg_1 = 1;*/
        unsigned char used = 1;
        // mark the start tgt quad as used, so it will not come back again
        rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
        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 };
            rval                                = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
#ifdef ENABLE_DEBUG
            if( dbg_1 )
            {
                std::cout << "Next: neighbors for current tgt ";
                for( int kk = 0; kk < nsidesTgt; kk++ )
                {
                    if( tgtNeighbors[kk] > 0 )
                        std::cout << 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;
                rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval );                     // status 0 is unused
                if( 1 == status ) tgtNeighbors[j] = 0;  // so will not look anymore on this side of tgt
            }
 
#ifdef ENABLE_DEBUG
            if( dbg_1 )
            {
                std::cout << "reset sources: ";
                for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
                    std::cout << mb->id_from_handle( *itr ) << " ";
                std::cout << std::endl;
            }
#endif
            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();
            toResetSrcs.clear();  // empty the range of used srcs, will have to be set unused again,
            // at the end of tgt element processing
            toResetSrcs.insert( currentSrc );
            // mb2->set_tag_data
            std::queue< EntityHandle > localSrc;
            localSrc.push( 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 queue will contain all possible src cells that cover
            //   the current tgt cell
            while( !localSrc.empty() )
            {
                //
                EntityHandle srcT = localSrc.front();
                localSrc.pop();
                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
                rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
                                                            nsidesTgt );MB_CHK_ERR( rval );
                if( nP > 0 )
                {
#ifdef ENABLE_DEBUG
                    if( dbg_1 )
                    {
                        for( int k = 0; k < 3; k++ )
                        {
                            std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
                        }
                    }
#endif
 
                    // intersection found: output P and original triangles if nP > 2
                    EntityHandle neighbors[MAXEDGES] = { 0 };
                    rval                             = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
                    if( rval != MB_SUCCESS )
                    {
                        std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
                        return MB_FAILURE;
                    }
 
                    // add neighbors to the localSrc queue, if they are not marked
                    for( int nn = 0; nn < nsidesSrc; nn++ )
                    {
                        EntityHandle neighbor = neighbors[nn];
                        if( neighbor > 0 && nb[nn] > 0 )  // advance across src boundary nn
                        {
                            if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
                            {
                                localSrc.push( neighbor );
#ifdef ENABLE_DEBUG
                                if( dbg_1 )
                                {
                                    std::cout << " local src elem " << mb->id_from_handle( neighbor )
                                              << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
                                    mb->list_entities( &neighbor, 1 );
                                }
#endif
                                toResetSrcs.insert( neighbor );
                            }
                        }
                    }
                    // 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
                        rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
                    }
 
                    recoveredArea += area;
                }
#ifdef ENABLE_DEBUG
                else if( dbg_1 )
                {
                    std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
                              << 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: " << 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 );
            // also, look at its neighbors, and add to the seeds a next one
 
            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;
                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;  ///
                    rval = computeIntersectionBetweenTgtAndSrc(
                        /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
                    if( area > 0 )
                    {
                        tgtQueue.push( tgtNeigh );
                        srcQueue.push( nextB );
#ifdef ENABLE_DEBUG
                        if( dbg_1 )
                            std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
                                      << mb->id_from_handle( nextB ) << std::endl;
#endif
                        rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
                        break;  // so we are done with this side of tgt, we have found a proper next
                                // seed
                    }
                }
            }
 
        }  // end while (!tgtQueue.empty())
    }
#ifdef ENABLE_DEBUG
    if( dbg_1 )
    {
        for( int k = 0; k < 6; k++ )
            mout_1[k].close();
    }
#endif
    // 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
    rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
#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(), moab::Range::find(), findNodes(), moab::Range::front(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::id_from_handle(), moab::ID_FROM_HANDLE(), imaskTag, moab::Range::insert(), moab::Interface::list_entities(), MAXEDGES, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, mbs1, mbs2, moab::Range::merge(), my_rank, nr, outSet, moab::AdaptiveKDTree::parse_options(), rs1, rs2, setup_tgt_cell(), size, moab::Range::size(), srcNeighTag, 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
	)

Definition at line 239 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;
    rval   = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
    rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
    // 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 (NULL)
    rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
    if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
    rval = filterByMask( rs1 );MB_CHK_ERR( rval );
    rval = filterByMask( rs2 );MB_CHK_ERR( rval );
    // 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
    rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt 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;
    rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
 
    rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
                               &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
 
    rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "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 );
    rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
                               MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
 
    for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
    {
        EntityHandle tgtCell = *rit;
        int num_nodes        = 0;
        rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
        // 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
        rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt 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 = NULL;
            int nnodes;
            rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
            while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
                nnodes--;
            rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get 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;
    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;
    rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
 
    for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
    {
        EntityHandle tcell = *it;
        // find vertex positions
        const EntityHandle* conn = NULL;
        int nnodes               = 0;
        rval                     = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
        // 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 );
        rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
 
        // 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();
            rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
 
            for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
            {
                Range tmp;
                rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
 
                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 << " id from handle "
                      << 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;
            rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
            if( area > 0 )
            {
                if( nP > 1 )
                {  // this will also construct triangles/polygons in the new mesh, if needed
                    rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
                }
                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
    rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
#endif
 
    this->clean();
    return MB_SUCCESS;
}

References moab::Range::begin(), box_error, moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), 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(), moab::Interface::id_from_handle(), imaskTag, moab::Interface::INTERSECT, max_edges_1, max_edges_2, MAXEDGES, mb, MB_CHK_ERR, 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, nr, 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().

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 253 of file Intx2Mesh.hpp.

box_error

double moab::Intx2Mesh::box_error

protected

Definition at line 254 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 268 of file Intx2Mesh.hpp.

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

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 250 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 251 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 248 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::Intx2MeshOnSphere::findNodes(), 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 257 of file Intx2Mesh.hpp.

localRoot

EntityHandle moab::Intx2Mesh::localRoot

protected

Definition at line 256 of file Intx2Mesh.hpp.

max_edges_1

int moab::Intx2Mesh::max_edges_1

protected

Definition at line 266 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 267 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(), filterByMask(), FindMaxEdgesInSet(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), intersect_meshes(), intersect_meshes_kdtree(), 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(), intersect_meshes(), and intersect_meshes_kdtree().

my_rank

unsigned int moab::Intx2Mesh::my_rank

protected

Definition at line 258 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::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