moab::Intx2MeshInPlane Class Reference

#include <Intx2MeshInPlane.hpp>

Inheritance diagram for moab::Intx2MeshInPlane:

Collaboration diagram for moab::Intx2MeshInPlane:

Public Member Functions
	Intx2MeshInPlane (Interface *mbimpl)
virtual	~Intx2MeshInPlane ()
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)
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 ()
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)

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
const EntityHandle *	tgtConn
	for coverage mesh, will store the original sender More...
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 15 of file Intx2MeshInPlane.hpp.

Constructor & Destructor Documentation

Intx2MeshInPlane()

moab::Intx2MeshInPlane::Intx2MeshInPlane

(

Interface *

mbimpl

)

Definition at line 15 of file Intx2MeshInPlane.cpp.

: Intx2Mesh( mbimpl ) {}

~Intx2MeshInPlane()

moab::Intx2MeshInPlane::~Intx2MeshInPlane ( )

virtual

Definition at line 17 of file Intx2MeshInPlane.cpp.

{}

Member Function Documentation

computeIntersectionBetweenTgtAndSrc()

ErrorCode moab::Intx2MeshInPlane::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 47 of file Intx2MeshInPlane.cpp.

{
 
 int num_nodes = 0;
 ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
 
 nsSrc = num_nodes;
 rval = mb->get_coords( srcConn, num_nodes, &( 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 )
 {
 setup_tgt_cell( tgt, nsTgt ); // we do not need area here
 // look at the boxes formed with vertices; if they are far away, return false early
 if( !GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsSrc, box_error ) )
 return MB_SUCCESS; // no error, but no intersection, decide early to get out
 }
#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 < nsSrc; j++ )
 {
 std::cout << srcCoords[j] << "\n";
 }
 mb->list_entities( &tgt, 1 );
 mb->list_entities( &src, 1 );
 }
#endif
 
 for( int j = 0; j < nsSrc; j++ )
 {
 srcCoords2D[2 * j] = srcCoords[j][0]; // x coordinate,
 srcCoords2D[2 * j + 1] = srcCoords[j][1]; // y coordinate
 }
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 // std::cout << "gnomonic plane: " << plane << "\n";
 std::cout << " tgt \n";
 for( int j = 0; j < nsTgt; j++ )
 {
 std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n ";
 }
 std::cout << " src\n";
 for( int j = 0; j < nsSrc; j++ )
 {
 std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
 }
 }
#endif
 
 rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 for( int k = 0; k < 3; k++ )
 {
 std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
 }
 }
#endif
 
 int side[MAXEDGES] = { 0 }; // this refers to what side? src or tgt?
 int extraPoints =
 IntxUtils::borderPointsOfXinY2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
 if( extraPoints >= 1 )
 {
 for( int k = 0; k < nsSrc; k++ )
 {
 if( side[k] )
 {
 // this means that vertex k of src 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 + nsSrc - 1 ) % nsSrc] =
 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;
 }
 }
 }
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 for( int k = 0; k < 3; k++ )
 {
 std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
 }
 }
#endif
 nP += extraPoints;
 
 extraPoints =
 IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsSrc, &( P[2 * nP] ), side, epsilon_area );
 if( extraPoints >= 1 )
 {
 for( int k = 0; k < nsTgt; k++ )
 {
 if( side[k] )
 {
 // this is to mark that tgt 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
 }
 }
 }
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 for( int k = 0; k < 3; k++ )
 {
 std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
 }
 }
#endif
 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] );
 }
 
 return MB_SUCCESS; // no error
}

References moab::IntxUtils::area2D(), moab::IntxUtils::borderPointsOfXinY2(), moab::GeomUtil::bounding_boxes_overlap(), moab::Intx2Mesh::box_error, moab::IntxUtils::EdgeIntersections2(), moab::Intx2Mesh::epsilon_1, moab::Intx2Mesh::epsilon_area, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::id_from_handle(), moab::Interface::list_entities(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_ERR, MB_SUCCESS, 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::Intx2MeshInPlane::findNodes ( EntityHandle  tgt, int  nsTgt, EntityHandle  src, int  nsSrc, double *  iP, int  nP  )

virtual

Implements moab::Intx2Mesh.

Definition at line 206 of file Intx2MeshInPlane.cpp.

{
 // except for gnomonic projection, everything is the same as spherical intx
 // start copy
 // first of all, check against tgt and src vertices
 //
#ifdef ENABLE_DEBUG
 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 tgt triangle; the extra points will be on those edges, saved as
 // lists (unordetgt)
 
 // first get the list of edges adjacent to the tgt 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 tgt tag" );
 // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
 // potatoes
 
 // these will be in the new mesh, mbOut
 // some of them will be handles to the initial vertices from src or tgt meshes (lagr or euler)
 
 EntityHandle* foundIds = new EntityHandle[nP];
 for( int i = 0; i < nP; i++ )
 {
 double* pp = &iP[2 * i]; // iP+2*i
 //
 CartVect pos( pp[0], pp[1], 0. );
 int found = 0;
 // first, are they on vertices from tgt or src?
 // priority is the tgt 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 ENABLE_DEBUG
 if( dbg_1 )
 std::cout << " tgt 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 < nsSrc && !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 ENABLE_DEBUG
 if( dbg_1 )
 std::cout << " src node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2 << " \n";
#endif
 }
 }
 if( !found )
 {
 // find the edge it belongs, first, on the tgt element
 //
 for( j = 0; j < nsTgt; j++ )
 {
 int j1 = ( j + 1 ) % nsTgt;
 double area = IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp );
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 std::cout << " edge " << j << ": " << mb->id_from_handle( adjTgtEdges[j] ) << " " << tgtConn[j]
 << " " << tgtConn[j1] << " area : " << area << "\n";
#endif
 if( fabs( area ) < epsilon_1 / 2 )
 {
 // 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[j] );
 if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
 {
 std::cerr << " error in adjacent tgt edge: " << mb->id_from_handle( adjTgtEdges[j] ) << "\n";
 delete[] foundIds;
 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 )
 {
 CartVect* coords1 = new CartVect[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 d3 = ( pos - coords1[k] ).length_squared();
 if( d3 < epsilon_1 )
 {
 found = 1;
 foundIds[i] = ( *expts )[k];
#ifdef ENABLE_DEBUG
 if( dbg_1 ) std::cout << " found node:" << foundIds[i] << std::endl;
#endif
 }
 }
 delete[] coords1;
 }
 
 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->create_vertex( pos.array(), outNode );
 ( *expts ).push_back( outNode );
 foundIds[i] = outNode;
 found = 1;
#ifdef ENABLE_DEBUG
 if( dbg_1 ) std::cout << " new node: " << outNode << std::endl;
#endif
 }
 }
 }
 }
 if( !found )
 {
 std::cout << " tgt polygon: ";
 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 tgt polygon " << *pp << " " << pp[1] << " tgt polygon "
 << mb->id_from_handle( tgt ) << " \n";
 delete[] foundIds;
 return MB_FAILURE;
 }
 }
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 std::cout << " candidate polygon: nP " << nP << "\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)
 correct_polygon( foundIds, 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->create_element( MBPOLYGON, foundIds, nP, polyNew );
 mb->add_entities( outSet, &polyNew, 1 );
 
 // tag it with the index ids from tgt and src sets
 int id = rs1.index( src ); // index starts from 0
 mb->tag_set_data( srcParentTag, &polyNew, 1, &id );
 id = rs2.index( tgt );
 mb->tag_set_data( tgtParentTag, &polyNew, 1, &id );
 
 counting++;
 mb->tag_set_data( countTag, &polyNew, 1, &counting );
 
#ifdef ENABLE_DEBUG
 if( dbg_1 )
 {
 
 std::cout << "Count: " << counting + 1 << "\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 << "\n";
 std::vector< CartVect > posi( nP );
 mb->get_coords( foundIds, nP, &( posi[0][0] ) );
 for( int i1 = 0; i1 < nP; i1++ )
 std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << posi[i1] << "\n";
 
 std::stringstream fff;
 fff << "file0" << counting << ".vtk";
 mb->write_mesh( fff.str().c_str(), &outSet, 1 );
 }
#endif
 }
 delete[] foundIds;
 foundIds = NULL;
 return MB_SUCCESS;
 // end copy
}

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::Interface::id_from_handle(), moab::Range::index(), length_squared(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, moab::Intx2Mesh::neighTgtEdgeTag, moab::Intx2Mesh::outSet, moab::Intx2Mesh::rs1, moab::Intx2Mesh::rs2, 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().

setup_tgt_cell()

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

virtual

Implements moab::Intx2Mesh.

Definition at line 19 of file Intx2MeshInPlane.cpp.

{
 // the points will be at most ?; they will describe a convex patch, after the points will be
 // ordered and collapsed (eliminate doubles) the area is not really required get coordinates of
 // the tgt quad
 double cellArea = 0;
 int num_nodes;
 ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );
 if( MB_SUCCESS != rval ) return 1.; // it should be an error
 
 nsTgt = num_nodes;
 
 rval = mb->get_coords( tgtConn, num_nodes, &( tgtCoords[0][0] ) );
 if( MB_SUCCESS != rval ) return 1.; // it should be an error
 
 for( int j = 0; j < nsTgt; j++ )
 {
 // populate coords in the plane for intersection
 // they should be oriented correctly, positively
 tgtCoords2D[2 * j] = tgtCoords[j][0]; // x coordinate,
 tgtCoords2D[2 * j + 1] = tgtCoords[j][1]; // y coordinate
 }
 for( int j = 1; j < nsTgt - 1; j++ )
 cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
 
 return cellArea;
}

References moab::IntxUtils::area2D(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Intx2Mesh::mb, MB_SUCCESS, moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.

Referenced by computeIntersectionBetweenTgtAndSrc().

---

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

• Intx2MeshInPlane.hpp
• Intx2MeshInPlane.cpp