moab::Intx2MeshOnSphere Class Reference

#include <Intx2MeshOnSphere.hpp>

Inheritance diagram for moab::Intx2MeshOnSphere:

Collaboration diagram for moab::Intx2MeshOnSphere:

Public Member Functions
	Intx2MeshOnSphere (Interface *mbimpl, IntxAreaUtils::AreaMethod amethod=IntxAreaUtils::lHuiller)
virtual	~Intx2MeshOnSphere ()
void	set_radius_source_mesh (double radius)
void	set_radius_destination_mesh (double radius)
double	setup_tgt_cell (EntityHandle tgt, int &nsTgt)
ErrorCode	computeIntersectionBetweenTgtAndSrc (EntityHandle tgt, EntityHandle src, double *P, int &nP, double &area, int markb[MAXEDGES], int markr[MAXEDGES], int &nsSrc, int &nsTgt, bool check_boxes_first=false)
ErrorCode	findNodes (EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double *iP, int nP)
ErrorCode	update_tracer_data (EntityHandle out_set, Tag &tagElem, Tag &tagArea)
Public Member Functions inherited from moab::Intx2Mesh
	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	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)

Public Attributes
const IntxAreaUtils::AreaMethod	areaMethod

Private Attributes
int	plane
double	Rsrc
double	Rdest

Additional Inherited Members
Protected Attributes inherited from moab::Intx2Mesh
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 16 of file Intx2MeshOnSphere.hpp.

Constructor & Destructor Documentation

Intx2MeshOnSphere()

moab::Intx2MeshOnSphere::Intx2MeshOnSphere	(	Interface *	mbimpl,
		IntxAreaUtils::AreaMethod	amethod = IntxAreaUtils::lHuiller
	)

Definition at line 28 of file Intx2MeshOnSphere.cpp.

    : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
{
}

~Intx2MeshOnSphere()

moab::Intx2MeshOnSphere::~Intx2MeshOnSphere ( )

virtual

Definition at line 33 of file Intx2MeshOnSphere.cpp.

{}

Member Function Documentation

computeIntersectionBetweenTgtAndSrc()

ErrorCode moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc ( EntityHandle  tgt, EntityHandle  src, double *  P, int &  nP, double &  area, int  markb[MAXEDGES], int  markr[MAXEDGES], int &  nsSrc, int &  nsTgt, bool  check_boxes_first = false  )

virtual

Implements moab::Intx2Mesh.

Definition at line 78 of file Intx2MeshOnSphere.cpp.

{
    // the area will be used from now on, to see how well we fill the target cell with polygons
    // the points will be at most 40; they will describe a convex patch, after the points will be
    // ordered and collapsed (eliminate doubles)
 
    // CartVect srccoords[4];
    int num_nodes  = 0;
    ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
    nsBlue = num_nodes;
    // account for possible padded polygons
    while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
        nsBlue--;
    rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
 
    area = 0.;
    nP   = 0;  // number of intersection points we are marking the boundary of src!
    if( check_boxes_first )
    {
        // look at the boxes formed with vertices; if they are far away, return false early
        // make sure the target is setup already
        setup_tgt_cell( tgt, nsTgt );  // we do not need area here
        // use here gnomonic plane (plane) to see where source is
        bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
        int planeb;
        CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
        IntxUtils::decide_gnomonic_plane( mid3, planeb );
        if( !overlap3d && ( plane != planeb ) )  // plane was set at setup_tgt_cell
            return MB_SUCCESS;                   // no error, but no intersection, decide early to get out
        // if same plane, still check for gnomonic plane in 2d
        // if no overlap in 2d, get out
        if( !overlap3d && plane == planeb )  // CHECK 2D too
        {
            for( int j = 0; j < nsBlue; j++ )
            {
                rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j],
                                                       srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
            }
            bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
            if( !overlap2d ) return MB_SUCCESS;  // we are sure they are not overlapping in 2d , either
        }
    }
#ifdef ENABLE_DEBUG
    if( dbg_1 )
    {
        std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
        for( int j = 0; j < nsTgt; j++ )
        {
            std::cout << tgtCoords[j] << "\n";
        }
        std::cout << "src " << mb->id_from_handle( src ) << "\n";
        for( int j = 0; j < nsBlue; j++ )
        {
            std::cout << srcCoords[j] << "\n";
        }
        mb->list_entities( &tgt, 1 );
        mb->list_entities( &src, 1 );
    }
#endif
 
    for( int j = 0; j < nsBlue; j++ )
    {
        rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
    }
 
#ifdef ENABLE_DEBUG
    if( dbg_1 )
    {
        std::cout << "gnomonic plane: " << plane << "\n";
        std::cout << " target                                src\n";
        for( int j = 0; j < nsTgt; j++ )
        {
            std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
        }
        for( int j = 0; j < nsBlue; j++ )
        {
            std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
        }
    }
#endif
 
    rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
 
    int side[MAXEDGES] = { 0 };  // this refers to what side? source or tgt?
    int extraPoints =
        IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
    if( extraPoints >= 1 )
    {
        for( int k = 0; k < nsBlue; k++ )
        {
            if( side[k] )
            {
                // this means that vertex k of source is inside convex tgt; mark edges k-1 and k in
                // src,
                //   as being "intersected" by tgt; (even though they might not be intersected by
                //   other edges, the fact that their apex is inside, is good enough)
                markb[k] = 1;
                markb[( k + nsBlue - 1 ) % nsBlue] =
                    1;  // it is the previous edge, actually, but instead of doing -1, it is
                // better to do modulo +3 (modulo 4)
                // null side b for next call
                side[k] = 0;
            }
        }
    }
    nP += extraPoints;
 
    extraPoints =
        IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
    if( extraPoints >= 1 )
    {
        for( int k = 0; k < nsTgt; k++ )
        {
            if( side[k] )
            {
                // this is to mark that target edges k-1 and k are intersecting src
                markr[k] = 1;
                markr[( k + nsTgt - 1 ) % nsTgt] =
                    1;  // it is the previous edge, actually, but instead of doing -1, it is
                // better to do modulo +3 (modulo 4)
                // null side b for next call
            }
        }
    }
    nP += extraPoints;
 
    // now sort and orient the points in P, such that they are forming a convex polygon
    // this will be the foundation of our new mesh
    // this works if the polygons are convex
    IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 );  // nP should be at most 8 in the end ?
    // if there are more than 3 points, some area will be positive
 
    if( nP >= 3 )
    {
        for( int k = 1; k < nP - 1; k++ )
            area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
#ifdef CHECK_CONVEXITY
        // each edge should be large enough that we can compute angles between edges
        for( int k = 0; k < nP; k++ )
        {
            int k1              = ( k + 1 ) % nP;
            int k2              = ( k1 + 1 ) % nP;
            double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
            if( orientedArea < 0 )
            {
                std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
                          << " " << src << " \n";
            }
        }
#endif
    }
 
    return MB_SUCCESS;  // no error
}

References moab::IntxUtils::area2D(), moab::IntxUtils::borderPointsOfXinY2(), moab::GeomUtil::bounding_boxes_overlap(), moab::GeomUtil::bounding_boxes_overlap_2d(), moab::Intx2Mesh::box_error, moab::IntxUtils::decide_gnomonic_plane(), moab::IntxUtils::EdgeIntersections2(), moab::Intx2Mesh::epsilon_1, moab::Intx2Mesh::epsilon_area, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::IntxUtils::gnomonic_projection(), moab::Interface::id_from_handle(), moab::Interface::list_entities(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_ERR, MB_SUCCESS, plane, Rsrc, setup_tgt_cell(), moab::IntxUtils::SortAndRemoveDoubles2(), moab::Intx2Mesh::srcConn, moab::Intx2Mesh::srcCoords, moab::Intx2Mesh::srcCoords2D, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.

findNodes()

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

virtual

Implements moab::Intx2Mesh.

Definition at line 247 of file Intx2MeshOnSphere.cpp.

{
#ifdef ENABLE_DEBUG
    // first of all, check against target and source vertices
    //
    if( dbg_1 )
    {
        std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
                  << "\n";
        for( int n = 0; n < nP; n++ )
            std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
    }
#endif
 
    // get the edges for the target triangle; the extra points will be on those edges, saved as
    // lists (unordered)
 
    // first get the list of edges adjacent to the target cell
    // use the neighTgtEdgeTag
    EntityHandle adjTgtEdges[MAXEDGES];
    ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge target tag" );
    // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
    // potatoes some of them will be handles to the initial vertices from source or target meshes
 
    std::vector< EntityHandle > foundIds;
    foundIds.resize( nP );
#ifdef CHECK_CONVEXITY
    int npBefore1 = nP;
    int oldNodes  = 0;
    int otherIntx = 0;
    moab::IntxAreaUtils areaAdaptor;
#endif
    for( int i = 0; i < nP; i++ )
    {
        double* pp = &iP[2 * i];  // iP+2*i
        // project the point back on the sphere
        CartVect pos;
        IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
        int found = 0;
        // first, are they on vertices from target or src?
        // priority is the target mesh (mb2?)
        int j                = 0;
        EntityHandle outNode = (EntityHandle)0;
        for( j = 0; j < nsTgt && !found; j++ )
        {
            // int node = tgtTri.v[j];
            double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
            if( d2 < epsilon_1 / 1000 ) // two orders of magnitude smaller than it should, to avoid concave polygons
            {
 
                foundIds[i] = tgtConn[j];  // no new node
                found       = 1;
#ifdef CHECK_CONVEXITY
                oldNodes++;
#endif
#ifdef ENABLE_DEBUG
                if( dbg_1 )
                    std::cout << "  target node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
                              << " 2d coords:" << tgtCoords2D[2 * j] << "  " << tgtCoords2D[2 * j + 1] << " d2: " << d2
                              << " \n";
#endif
            }
        }
 
        for( j = 0; j < nsBlue && !found; j++ )
        {
            // int node = srcTri.v[j];
            double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
            if( d2 < epsilon_1 / 1000 )
            {
                // suspect is srcConn[j] corresponding in mbOut
 
                foundIds[i] = srcConn[j];  // no new node
                found       = 1;
#ifdef CHECK_CONVEXITY
                oldNodes++;
#endif
#ifdef ENABLE_DEBUG
                if( dbg_1 )
                    std::cout << "  source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
                              << " \n";
#endif
            }
        }
 
        if( !found )
        {
            // find the edge it belongs, first, on the red element
            // look at the minimum area, not at the first below some tolerance
            double minArea = 1.e+38;
            int index_min  = -1;
            for( j = 0; j < nsTgt; j++ )
            {
                int j1      = ( j + 1 ) % nsTgt;
                double area = fabs( IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp ) );
                // how to check if pp is between redCoords2D[j] and redCoords2D[j1] ?
                // they should form a straight line; the sign should be -1
                double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
                                IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
                                IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
                if( area < minArea && checkx < 2 * epsilon_1 )  // round off error or not?
                {
                    index_min = j;
                    minArea   = area;
                }
            }
            // verify that index_min is valid
            assert( index_min >= 0 );
 
            if( minArea < epsilon_1 / 2 )  // we found the smallest area, so we think we found the
                                           // target edge it belongs
            {
                // found the edge; now find if there is a point in the list here
                // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
                int indx = TgtEdges.index( adjTgtEdges[index_min] );
                if( indx < 0 )  // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
                {
                    std::cerr << " error in adjacent target edge: " << mb->id_from_handle( adjTgtEdges[index_min] )
                              << "\n";
                    return MB_FAILURE;
                }
                std::vector< EntityHandle >* expts = extraNodesVec[indx];
                // if the points pp is between extra points, then just give that id
                // if not, create a new point, (check the id)
                // get the coordinates of the extra points so far
                int nbExtraNodesSoFar = expts->size();
                if( nbExtraNodesSoFar > 0 )
                {
                    std::vector< CartVect > coords1;
                    coords1.resize( nbExtraNodesSoFar );
                    mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
                    // std::list<int>::iterator it;
                    for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
                    {
                        // int pnt = *it;
                        double d2 = ( pos - coords1[k] ).length();
                        if( d2 < 2 * epsilon_1 )  // is this below machine precision?
                        {
                            found       = 1;
                            foundIds[i] = ( *expts )[k];
#ifdef CHECK_CONVEXITY
                            otherIntx++;
#endif
                        }
                    }
                }
                if( !found )
                {
                    // create a new point in 2d (at the intersection)
                    // foundIds[i] = m_num2dPoints;
                    // expts.push_back(m_num2dPoints);
                    // need to create a new node in mbOut
                    // this will be on the edge, and it will be added to the local list
                    rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
                    ( *expts ).push_back( outNode );
                    // CID 181168; avoid leak storage error
                    rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
                    foundIds[i] = outNode;
                    found       = 1;
                }
            }
        }
        if( !found )
        {
            std::cout << " target quad: ";
            for( int j1 = 0; j1 < nsTgt; j1++ )
            {
                std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
            }
            std::cout << " a point pp is not on a target quad " << *pp << " " << pp[1] << " target quad "
                      << mb->id_from_handle( tgt ) << " \n";
            return MB_FAILURE;
        }
    }
#ifdef ENABLE_DEBUG
    if( dbg_1 )
    {
        std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
        for( int i1 = 0; i1 < nP; i1++ )
            std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
    }
#endif
    // first, find out if we have nodes collapsed; shrink them
    // we may have to reduce nP
    // it is possible that some nodes are collapsed after intersection only
    // nodes will always be in order (convex intersection)
#ifdef CHECK_CONVEXITY
    int npBefore2 = nP;
#endif
    correct_polygon( &foundIds[0], nP );
    // now we can build the triangles, from P array, with foundIds
    // we will put them in the out set
    if( nP >= 3 )
    {
        EntityHandle polyNew;
        rval = mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew );MB_CHK_ERR( rval );
        rval = mb->add_entities( outSet, &polyNew, 1 );MB_CHK_ERR( rval );
 
        // tag it with the global ids from target and source elements
        int globalID;
        rval = mb->tag_get_data( gid, &src, 1, &globalID );MB_CHK_ERR( rval );
        rval = mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
        // if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
        // " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";
        rval = mb->tag_get_data( gid, &tgt, 1, &globalID );MB_CHK_ERR( rval );
        rval = mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
        // if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
        // " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";
 
        counting++;
        rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
        if( orgSendProcTag )
        {
            int org_proc = -1;
            rval         = mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc );MB_CHK_ERR( rval );
            rval = mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc );MB_CHK_ERR( rval );  // yet another tag
        }
#ifdef CHECK_CONVEXITY
        // each edge should be large enough that we can compute angles between edges
        std::vector< double > coords;
        coords.resize( 3 * nP );
        rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
        std::vector< CartVect > posi( nP );
        rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
 
        for( int k = 0; k < nP; k++ )
        {
            int k1 = ( k + 1 ) % nP;
            int k2 = ( k1 + 1 ) % nP;
            double orientedArea = areaAdaptor. area_spherical_triangle( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
            if( orientedArea < 0 )
            {
                std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
                for( int i = 0; i < nP; i++ )
                {
                    int nexti         = ( i + 1 ) % nP;
                    double lengthEdge = ( posi[i] - posi[nexti] ).length();
                    std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n";
                }
                std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n";
 
                std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
                          << " target, src:" << tgt << " " << src << " \n";
            }
        }
#endif
 
#ifdef ENABLE_DEBUG
        if( dbg_1 )
        {
            std::cout << "Counting: " << counting << "\n";
            std::cout << " polygon " << mb->id_from_handle( polyNew ) << "  nodes: " << nP << " :";
            for( int i1 = 0; i1 < nP; i1++ )
                std::cout << " " << mb->id_from_handle( foundIds[i1] );
            std::cout << " plane: " << plane << "\n";
            std::vector< CartVect > posi( nP );
            mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
            for( int i1 = 0; i1 < nP; i1++ )
                std::cout << foundIds[i1] << " " << posi[i1] << "\n";
 
            std::stringstream fff;
            fff << "file0" << counting << ".vtk";
            rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
        }
#endif
    }
    // else {
    //   std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";
    // }
    // disable_debug();
    return MB_SUCCESS;
}

References moab::Interface::add_entities(), moab::IntxUtils::area2D(), moab::CartVect::array(), moab::Intx2Mesh::correct_polygon(), moab::Intx2Mesh::counting, moab::Intx2Mesh::countTag, moab::Interface::create_element(), moab::Interface::create_vertex(), moab::IntxUtils::dist2(), moab::Intx2Mesh::epsilon_1, ErrorCode, moab::Intx2Mesh::extraNodesVec, moab::Interface::get_coords(), moab::Intx2Mesh::gid, moab::Interface::id_from_handle(), moab::Range::index(), length(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, moab::Intx2Mesh::my_rank, moab::Intx2Mesh::neighTgtEdgeTag, moab::Intx2Mesh::orgSendProcTag, moab::Intx2Mesh::outSet, plane, Rdest, moab::IntxUtils::reverse_gnomonic_projection(), moab::Intx2Mesh::srcConn, moab::Intx2Mesh::srcCoords2D, moab::Intx2Mesh::srcParentTag, moab::Interface::tag_get_data(), moab::Interface::tag_set_data(), moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords2D, moab::Intx2Mesh::TgtEdges, moab::Intx2Mesh::tgtParentTag, and moab::Interface::write_mesh().

set_radius_destination_mesh()

void moab::Intx2MeshOnSphere::set_radius_destination_mesh ( double  radius )

inline

Definition at line 27 of file Intx2MeshOnSphere.hpp.

    {
        Rdest = radius;
    }

References Rdest.

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

set_radius_source_mesh()

void moab::Intx2MeshOnSphere::set_radius_source_mesh ( double  radius )

inline

Definition at line 23 of file Intx2MeshOnSphere.hpp.

    {
        Rsrc = radius;
    }

References Rsrc.

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

setup_tgt_cell()

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

virtual

Implements moab::Intx2Mesh.

Definition at line 38 of file Intx2MeshOnSphere.cpp.

{
 
    // get coordinates of the target quad, to decide the gnomonic plane
    double cellArea = 0;
 
    int num_nodes;
    ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
 
    nsTgt = num_nodes;
    // account for possible padded polygons
    while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
        nsTgt--;
 
    // CartVect coords[4];
    rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
 
    CartVect middle = tgtCoords[0];
    for( int i = 1; i < nsTgt; i++ )
        middle += tgtCoords[i];
    middle = 1. / nsTgt * middle;
 
    IntxUtils::decide_gnomonic_plane( middle, plane );  // output the plane
    for( int j = 0; j < nsTgt; j++ )
    {
        // populate coords in the plane for intersection
        // they should be oriented correctly, positively
        rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
    }
 
    for( int j = 1; j < nsTgt - 1; j++ )
        cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
 
    // take target coords in order and compute area in plane
    return cellArea;
}

References moab::IntxUtils::area2D(), moab::IntxUtils::decide_gnomonic_plane(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::IntxUtils::gnomonic_projection(), moab::Intx2Mesh::mb, MB_CHK_ERR_RET_VAL, plane, Rdest, moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.

Referenced by computeIntersectionBetweenTgtAndSrc().

update_tracer_data()

ErrorCode moab::Intx2MeshOnSphere::update_tracer_data	(	EntityHandle	out_set,
		Tag &	tagElem,
		Tag &	tagArea
	)

TODO: VSM: Its unclear whether we need the source or destination radius here.

Definition at line 520 of file Intx2MeshOnSphere.cpp.

{
    EntityHandle dum = 0;
 
    Tag corrTag;
    ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE,
                                         &dum );  // it should have been created
    MB_CHK_SET_ERR( rval, "can't get correlation tag" );
 
    // get all polygons out of out_set; then see where are they coming from
    Range polys;
    rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
 
    // rs2 is the target range, arrival; rs1 is src, departure;
    // there is a connection between rs1 and rs2, through the corrTag
    // corrTag is __correlation
    // basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly);
    // also,  mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly);
    // we start from rs2 existing, then we have to update something
 
    // tagElem will have multiple tracers
    int numTracers = 0;
    rval           = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
    if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
 
    std::vector< double > currentVals( rs2.size() * numTracers );
    rval = mb->tag_get_data( tagElem, rs2, &currentVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" );
 
    // create new tuple list for tracers to other processors, from remote_cells
#ifdef MOAB_HAVE_MPI
    if( remote_cells )
    {
        int n = remote_cells->get_n();
        if( n > 0 )
        {
            remote_cells_with_tracers = new TupleList();
            remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
                                                   n );  // tracers are in these tuples
            remote_cells_with_tracers->enableWriteAccess();
            for( int i = 0; i < n; i++ )
            {
                remote_cells_with_tracers->vi_wr[2 * i]     = remote_cells->vi_wr[2 * i];
                remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
                //    remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
                remote_cells_with_tracers->vul_wr[i] =
                    remote_cells->vul_wr[i];  // this is the corresponding target cell (arrival)
                for( int k = 0; k < numTracers; k++ )
                    remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0;  // initialize tracers to be transported
                remote_cells_with_tracers->inc_n();
            }
        }
        delete remote_cells;
        remote_cells = NULL;
    }
#endif
    // for each polygon, we have 2 indices: target and source parents
    // we need index source to update index tgt?
    std::vector< double > newValues( rs2.size() * numTracers,
                                     0. );  // initialize with 0 all of them
    // area of the polygon * conc on target (old) current quantity
    // finally, divide by the area of the tgt
    double check_intx_area = 0.;
    moab::IntxAreaUtils intxAreas( this->areaMethod );  // use_lHuiller = true
    for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
    {
        EntityHandle poly = *it;
        int srcIndex, tgtIndex;
        rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
 
        EntityHandle src = rs1[srcIndex - 1];  // big assumption, it should work for meshes where global id is the same
        // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes
        rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" );
        // EntityHandle target = rs2[tgtIndex];
        // big assumption here, target and source are "parallel" ;we should have an index from
        // source to target (so a deformed source corresponds to an arrival tgt)
        /// TODO: VSM: Its unclear whether we need the source or destination radius here.
        double radius = Rsrc;
        double areap  = intxAreas.area_spherical_element( mb, poly, radius );
        check_intx_area += areap;
        // so the departure cell at time t (srcIndex) covers a portion of a tgtCell
        // that quantity will be transported to the tgtCell at time t+dt
        // the source corresponds to a target arrival
        EntityHandle tgtArr;
        rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
        if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
        {
#ifdef MOAB_HAVE_MPI
            if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
            // maybe the element is remote, from another processor
            int global_id_src;
            rval = mb->tag_get_data( gid, &src, 1, &global_id_src );MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding source gid" );
            // find the
            int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
            if( index_in_remote == -1 )
                MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
            for( int k = 0; k < numTracers; k++ )
                remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
                    currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
#endif
        }
        else if( MB_SUCCESS == rval )
        {
            int arrTgtIndex = rs2.index( tgtArr );
            if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
            for( int k = 0; k < numTracers; k++ )
                newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
        }
 
        else
            MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " );
    }
    // now, send back the remote_cells_with_tracers to the processors they came from, with the
    // updated values for the tracer mass in a cell
#ifdef MOAB_HAVE_MPI
    if( remote_cells_with_tracers )
    {
        // so this means that some cells will be sent back with tracer info to the procs they were
        // sent from
        ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
        // now, look at the global id, find the proper "tgt" cell with that index and update its
        // mass
        // remote_cells->print("remote cells after routing");
        int n = remote_cells_with_tracers->get_n();
        for( int j = 0; j < n; j++ )
        {
            EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j];  // entity handle sent back
            int arrTgtIndex      = rs2.index( tgtCell );
            if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
            for( int k = 0; k < numTracers; k++ )
                newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
        }
    }
#endif /* MOAB_HAVE_MPI */
    // now divide by target area (current)
    int j                = 0;
    Range::iterator iter = rs2.begin();
    void* data           = NULL;  // used for stored area
    int count            = 0;
    std::vector< double > total_mass_local( numTracers, 0. );
    while( iter != rs2.end() )
    {
        rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
        double* ptrArea = (double*)data;
        for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
        {
            for( int k = 0; k < numTracers; k++ )
            {
                total_mass_local[k] += newValues[j * numTracers + k];
                newValues[j * numTracers + k] /= ( *ptrArea );
            }
        }
    }
    rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" );
 
#ifdef MOAB_HAVE_MPI
    std::vector< double > total_mass( numTracers, 0. );
    double total_intx_area = 0;
    int mpi_err =
        MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
    if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
    // now reduce total area
    mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
    if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
    if( my_rank == 0 )
    {
        for( int k = 0; k < numTracers; k++ )
            std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n";
        std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
                  << total_intx_area << "\n";
    }
 
    if( remote_cells_with_tracers )
    {
        delete remote_cells_with_tracers;
        remote_cells_with_tracers = NULL;
    }
#else
    for( int k = 0; k < numTracers; k++ )
        std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n";
    std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
              << check_intx_area << "\n";
#endif
    return MB_SUCCESS;
}

References moab::IntxAreaUtils::area_spherical_element(), areaMethod, moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), ErrorCode, moab::Interface::get_entities_by_dimension(), moab::Intx2Mesh::gid, moab::Range::index(), moab::Intx2Mesh::mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TAG_NOT_FOUND, MB_TYPE_HANDLE, moab::Intx2Mesh::my_rank, moab::Intx2Mesh::rs1, moab::Intx2Mesh::rs2, Rsrc, moab::Range::size(), moab::Intx2Mesh::srcParentTag, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_get_length(), moab::Interface::tag_iterate(), moab::Interface::tag_set_data(), and moab::Intx2Mesh::tgtParentTag.

Member Data Documentation

areaMethod

const IntxAreaUtils::AreaMethod moab::Intx2MeshOnSphere::areaMethod

Definition at line 59 of file Intx2MeshOnSphere.hpp.

Referenced by update_tracer_data().

plane

int moab::Intx2MeshOnSphere::plane

private

Definition at line 62 of file Intx2MeshOnSphere.hpp.

Referenced by computeIntersectionBetweenTgtAndSrc(), findNodes(), and setup_tgt_cell().

Rdest

double moab::Intx2MeshOnSphere::Rdest

private

Definition at line 63 of file Intx2MeshOnSphere.hpp.

Referenced by findNodes(), set_radius_destination_mesh(), and setup_tgt_cell().

Rsrc

double moab::Intx2MeshOnSphere::Rsrc

private

Definition at line 63 of file Intx2MeshOnSphere.hpp.

Referenced by computeIntersectionBetweenTgtAndSrc(), set_radius_source_mesh(), and update_tracer_data().

---

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

• Intx2MeshOnSphere.hpp
• Intx2MeshOnSphere.cpp