#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 Attributes
const IntxAreaUtils::AreaMethod	areaMethod
Private Attributes
int	plane
double	Rsrc
double	Rdest

Detailed Description

Definition at line 16 of file Intx2MeshOnSphere.hpp.

Constructor & Destructor Documentation

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

Definition at line 26 of file Intx2MeshOnSphere.cpp.

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

moab::Intx2MeshOnSphere::~Intx2MeshOnSphere ( ) [virtual]

Definition at line 31 of file Intx2MeshOnSphere.cpp.

{}

Member Function Documentation

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 76 of file Intx2MeshOnSphere.cpp.

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.

{
    // 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
}

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

Implements moab::Intx2Mesh.

Definition at line 245 of file Intx2MeshOnSphere.cpp.

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().

{
#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;
#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 )
            {

                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 )
            {
                // 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 =
                area_spherical_triangle_lHuiller( &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;
}

void moab::Intx2MeshOnSphere::set_radius_destination_mesh ( double radius ) [inline]

Definition at line 27 of file Intx2MeshOnSphere.hpp.

References radius, and Rdest.

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

    {
        Rdest = radius;
    }

void moab::Intx2MeshOnSphere::set_radius_source_mesh ( double radius ) [inline]

Definition at line 23 of file Intx2MeshOnSphere.hpp.

References radius, and Rsrc.

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

    {
        Rsrc = radius;
    }

double moab::Intx2MeshOnSphere::setup_tgt_cell	(	EntityHandle	tgt,
		int &	nsTgt
	)		`[virtual]`

Implements moab::Intx2Mesh.

Definition at line 36 of file Intx2MeshOnSphere.cpp.

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().

{

    // 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;
}

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 518 of file Intx2MeshOnSphere.cpp.

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, MPI_COMM_WORLD, moab::Intx2Mesh::my_rank, radius, 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.

Referenced by compute_tracer_case1(), update_tracer(), and update_tracer_test().

{
    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;
}