#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.

33 {}

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;
     MB_CHK_ERR( mb->get_connectivity( src, srcConn, num_nodes ) );
     nsBlue = num_nodes;
     // account for possible padded polygons
     while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
         nsBlue--;
     MB_CHK_ERR( mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) ) );
  
     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++ )
             {
                 MB_CHK_ERR( IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j],
                                                             srcCoords2D[2 * j + 1] ) );
             }
             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++ )
     {
         MB_CHK_ERR(
             IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] ) );
     }
  
 #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
  
     MB_CHK_ERR( IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP ) );
  
     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, 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 248 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];
     MB_CHK_SET_ERR( mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) ), "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
                     MB_CHK_ERR( mb->create_vertex( pos.array(), outNode ) );
                     ( *expts ).push_back( outNode );
                     // CID 181168; avoid leak storage error
                     MB_CHK_ERR( mb->add_entities( outSet, &outNode, 1 ) );
                     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;
         MB_CHK_ERR( mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew ) );
         MB_CHK_ERR( mb->add_entities( outSet, &polyNew, 1 ) );
  
         // tag it with the global ids from target and source elements
         int globalID;
         MB_CHK_ERR( mb->tag_get_data( gid, &src, 1, &globalID ) );
         MB_CHK_ERR( mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID ) );
         // if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
         // " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";
         MB_CHK_ERR( mb->tag_get_data( gid, &tgt, 1, &globalID ) );
         MB_CHK_ERR( mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID ) );
         // if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
         // " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";
  
         counting++;
         MB_CHK_ERR( mb->tag_set_data( countTag, &polyNew, 1, &counting ) );
         if( orgSendProcTag )
         {
             int org_proc = -1;
             MB_CHK_ERR( mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc ) );
             MB_CHK_ERR( mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc ) );  // 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 );
         MB_CHK_ERR( mb->get_coords( &foundIds[0], nP, &coords[0] ) );
         std::vector< CartVect > posi( nP );
         MB_CHK_ERR( mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) ) );
  
         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";
             MB_CHK_ERR( mb->write_mesh( fff.str().c_str(), &outSet, 1 ) );
         }
 #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::IntxAreaUtils::area_spherical_triangle(), 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, 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 radius, and 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 radius, and 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 522 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;
     MB_CHK_SET_ERR( mb->get_entities_by_dimension( out_set, 2, polys ), "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;
     MB_CHK_SET_ERR( mb->tag_get_length( tagElem, numTracers ), "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 );
     MB_CHK_SET_ERR( mb->tag_get_data( tagElem, rs2, &currentVals[0] ), "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 = nullptr;
     }
 #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;
         MB_CHK_SET_ERR( mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex ), "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
         MB_CHK_SET_ERR( mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex ), "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;
             MB_CHK_SET_ERR( mb->tag_get_data( gid, &src, 1, &global_id_src ),
                             "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           = nullptr;  // used for stored area
     int count            = 0;
     std::vector< double > total_mass_local( numTracers, 0. );
     while( iter != rs2.end() )
     {
         MB_CHK_SET_ERR( mb->tag_iterate( tagArea, iter, rs2.end(), count, data ), "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 );
             }
         }
     }
     MB_CHK_SET_ERR( mb->tag_set_data( tagElem, rs2, &newValues[0] ), "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 = nullptr;
     }
 #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, 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.

Member Data Documentation

◆ areaMethod

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

Definition at line 60 of file Intx2MeshOnSphere.hpp.

Referenced by update_tracer_data().

◆ plane

int moab::Intx2MeshOnSphere::plane

private

Definition at line 63 of file Intx2MeshOnSphere.hpp.

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

◆ Rdest

double moab::Intx2MeshOnSphere::Rdest

private

Definition at line 64 of file Intx2MeshOnSphere.hpp.

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

◆ Rsrc

double moab::Intx2MeshOnSphere::Rsrc

private

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

Public Member Functions

Public Attributes

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Intx2MeshOnSphere()

◆ ~Intx2MeshOnSphere()

Member Function Documentation

◆ computeIntersectionBetweenTgtAndSrc()

◆ findNodes()

◆ set_radius_destination_mesh()

◆ set_radius_source_mesh()

◆ setup_tgt_cell()

◆ update_tracer_data()

Member Data Documentation

◆ areaMethod

◆ plane

◆ Rdest

◆ Rsrc