Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 ()
 
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)
 

Additional Inherited Members

- Protected Attributes inherited from moab::Intx2Mesh
Interfacemb
 
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 EntityHandletgtConn
 
const EntityHandlesrcConn
 
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.

15 : Intx2Mesh( mbimpl ) {}

◆ ~Intx2MeshInPlane()

moab::Intx2MeshInPlane::~Intx2MeshInPlane ( )
virtual

Definition at line 17 of file Intx2MeshInPlane.cpp.

17 {}

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.

57 {
58 
59  int num_nodes = 0;
60  ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
61 
62  nsSrc = num_nodes;
63  rval = mb->get_coords( srcConn, num_nodes, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
64 
65  area = 0.;
66  nP = 0; // number of intersection points we are marking the boundary of src!
67  if( check_boxes_first )
68  {
69  setup_tgt_cell( tgt, nsTgt ); // we do not need area here
70  // look at the boxes formed with vertices; if they are far away, return false early
72  return MB_SUCCESS; // no error, but no intersection, decide early to get out
73  }
74 #ifdef ENABLE_DEBUG
75  if( dbg_1 )
76  {
77  std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
78  for( int j = 0; j < nsTgt; j++ )
79  {
80  std::cout << tgtCoords[j] << "\n";
81  }
82  std::cout << "src " << mb->id_from_handle( src ) << "\n";
83  for( int j = 0; j < nsSrc; j++ )
84  {
85  std::cout << srcCoords[j] << "\n";
86  }
87  mb->list_entities( &tgt, 1 );
88  mb->list_entities( &src, 1 );
89  }
90 #endif
91 
92  for( int j = 0; j < nsSrc; j++ )
93  {
94  srcCoords2D[2 * j] = srcCoords[j][0]; // x coordinate,
95  srcCoords2D[2 * j + 1] = srcCoords[j][1]; // y coordinate
96  }
97 #ifdef ENABLE_DEBUG
98  if( dbg_1 )
99  {
100  // std::cout << "gnomonic plane: " << plane << "\n";
101  std::cout << " tgt \n";
102  for( int j = 0; j < nsTgt; j++ )
103  {
104  std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n ";
105  }
106  std::cout << " src\n";
107  for( int j = 0; j < nsSrc; j++ )
108  {
109  std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
110  }
111  }
112 #endif
113 
114  rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
115 #ifdef ENABLE_DEBUG
116  if( dbg_1 )
117  {
118  for( int k = 0; k < 3; k++ )
119  {
120  std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
121  }
122  }
123 #endif
124 
125  int side[MAXEDGES] = { 0 }; // this refers to what side? src or tgt?
126  int extraPoints =
127  IntxUtils::borderPointsOfXinY2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
128  if( extraPoints >= 1 )
129  {
130  for( int k = 0; k < nsSrc; k++ )
131  {
132  if( side[k] )
133  {
134  // this means that vertex k of src is inside convex tgt; mark edges k-1 and k in
135  // src,
136  // as being "intersected" by tgt; (even though they might not be intersected by
137  // other edges, the fact that their apex is inside, is good enough)
138  markb[k] = 1;
139  markb[( k + nsSrc - 1 ) % nsSrc] =
140  1; // it is the previous edge, actually, but instead of doing -1, it is
141  // better to do modulo +3 (modulo 4)
142  // null side b for next call
143  side[k] = 0;
144  }
145  }
146  }
147 #ifdef ENABLE_DEBUG
148  if( dbg_1 )
149  {
150  for( int k = 0; k < 3; k++ )
151  {
152  std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
153  }
154  }
155 #endif
156  nP += extraPoints;
157 
158  extraPoints =
159  IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsSrc, &( P[2 * nP] ), side, epsilon_area );
160  if( extraPoints >= 1 )
161  {
162  for( int k = 0; k < nsTgt; k++ )
163  {
164  if( side[k] )
165  {
166  // this is to mark that tgt edges k-1 and k are intersecting src
167  markr[k] = 1;
168  markr[( k + nsTgt - 1 ) % nsTgt] =
169  1; // it is the previous edge, actually, but instead of doing -1, it is
170  // better to do modulo +3 (modulo 4)
171  // null side b for next call
172  }
173  }
174  }
175 #ifdef ENABLE_DEBUG
176  if( dbg_1 )
177  {
178  for( int k = 0; k < 3; k++ )
179  {
180  std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
181  }
182  }
183 #endif
184  nP += extraPoints;
185 
186  // now sort and orient the points in P, such that they are forming a convex polygon
187  // this will be the foundation of our new mesh
188  // this works if the polygons are convex
189  IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 ); // nP should be at most 8 in the end ?
190  // if there are more than 3 points, some area will be positive
191 
192  if( nP >= 3 )
193  {
194  for( int k = 1; k < nP - 1; k++ )
195  area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
196  }
197 
198  return MB_SUCCESS; // no error
199 }

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.

207 {
208  // except for gnomonic projection, everything is the same as spherical intx
209  // start copy
210  // first of all, check against tgt and src vertices
211  //
212 #ifdef ENABLE_DEBUG
213  if( dbg_1 )
214  {
215  std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
216  << "\n";
217  for( int n = 0; n < nP; n++ )
218  std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
219  }
220 #endif
221 
222  // get the edges for the tgt triangle; the extra points will be on those edges, saved as
223  // lists (unordetgt)
224 
225  // first get the list of edges adjacent to the tgt cell
226  // use the neighTgtEdgeTag
227  EntityHandle adjTgtEdges[MAXEDGES];
228  ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge tgt tag" );
229  // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
230  // potatoes
231 
232  // these will be in the new mesh, mbOut
233  // some of them will be handles to the initial vertices from src or tgt meshes (lagr or euler)
234 
235  EntityHandle* foundIds = new EntityHandle[nP];
236  for( int i = 0; i < nP; i++ )
237  {
238  double* pp = &iP[2 * i]; // iP+2*i
239  //
240  CartVect pos( pp[0], pp[1], 0. );
241  int found = 0;
242  // first, are they on vertices from tgt or src?
243  // priority is the tgt mesh (mb2?)
244  int j = 0;
245  EntityHandle outNode = (EntityHandle)0;
246  for( j = 0; j < nsTgt && !found; j++ )
247  {
248  // int node = tgtTri.v[j];
249  double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
250  if( d2 < epsilon_1 )
251  {
252 
253  foundIds[i] = tgtConn[j]; // no new node
254  found = 1;
255 #ifdef ENABLE_DEBUG
256  if( dbg_1 )
257  std::cout << " tgt node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
258  << " 2d coords:" << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << " d2: " << d2
259  << " \n";
260 #endif
261  }
262  }
263 
264  for( j = 0; j < nsSrc && !found; j++ )
265  {
266  // int node = srcTri.v[j];
267  double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
268  if( d2 < epsilon_1 )
269  {
270  // suspect is srcConn[j] corresponding in mbOut
271 
272  foundIds[i] = srcConn[j]; // no new node
273  found = 1;
274 #ifdef ENABLE_DEBUG
275  if( dbg_1 )
276  std::cout << " src node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2 << " \n";
277 #endif
278  }
279  }
280  if( !found )
281  {
282  // find the edge it belongs, first, on the tgt element
283  //
284  for( j = 0; j < nsTgt; j++ )
285  {
286  int j1 = ( j + 1 ) % nsTgt;
287  double area = IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp );
288 #ifdef ENABLE_DEBUG
289  if( dbg_1 )
290  std::cout << " edge " << j << ": " << mb->id_from_handle( adjTgtEdges[j] ) << " " << tgtConn[j]
291  << " " << tgtConn[j1] << " area : " << area << "\n";
292 #endif
293  if( fabs( area ) < epsilon_1 / 2 )
294  {
295  // found the edge; now find if there is a point in the list here
296  // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
297  int indx = TgtEdges.index( adjTgtEdges[j] );
298  if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
299  {
300  std::cerr << " error in adjacent tgt edge: " << mb->id_from_handle( adjTgtEdges[j] ) << "\n";
301  delete[] foundIds;
302  return MB_FAILURE;
303  }
304  std::vector< EntityHandle >* expts = extraNodesVec[indx];
305  // if the points pp is between extra points, then just give that id
306  // if not, create a new point, (check the id)
307  // get the coordinates of the extra points so far
308  int nbExtraNodesSoFar = expts->size();
309  if( nbExtraNodesSoFar > 0 )
310  {
311  CartVect* coords1 = new CartVect[nbExtraNodesSoFar];
312  mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
313  // std::list<int>::iterator it;
314  for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
315  {
316  // int pnt = *it;
317  double d3 = ( pos - coords1[k] ).length_squared();
318  if( d3 < epsilon_1 )
319  {
320  found = 1;
321  foundIds[i] = ( *expts )[k];
322 #ifdef ENABLE_DEBUG
323  if( dbg_1 ) std::cout << " found node:" << foundIds[i] << std::endl;
324 #endif
325  }
326  }
327  delete[] coords1;
328  }
329 
330  if( !found )
331  {
332  // create a new point in 2d (at the intersection)
333  // foundIds[i] = m_num2dPoints;
334  // expts.push_back(m_num2dPoints);
335  // need to create a new node in mbOut
336  // this will be on the edge, and it will be added to the local list
337  mb->create_vertex( pos.array(), outNode );
338  ( *expts ).push_back( outNode );
339  foundIds[i] = outNode;
340  found = 1;
341 #ifdef ENABLE_DEBUG
342  if( dbg_1 ) std::cout << " new node: " << outNode << std::endl;
343 #endif
344  }
345  }
346  }
347  }
348  if( !found )
349  {
350  std::cout << " tgt polygon: ";
351  for( int j1 = 0; j1 < nsTgt; j1++ )
352  {
353  std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
354  }
355  std::cout << " a point pp is not on a tgt polygon " << *pp << " " << pp[1] << " tgt polygon "
356  << mb->id_from_handle( tgt ) << " \n";
357  delete[] foundIds;
358  return MB_FAILURE;
359  }
360  }
361 #ifdef ENABLE_DEBUG
362  if( dbg_1 )
363  {
364  std::cout << " candidate polygon: nP " << nP << "\n";
365  for( int i1 = 0; i1 < nP; i1++ )
366  std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
367  }
368 #endif
369  // first, find out if we have nodes collapsed; shrink them
370  // we may have to reduce nP
371  // it is possible that some nodes are collapsed after intersection only
372  // nodes will always be in order (convex intersection)
373  correct_polygon( foundIds, nP );
374  // now we can build the triangles, from P array, with foundIds
375  // we will put them in the out set
376  if( nP >= 3 )
377  {
378  EntityHandle polyNew;
379  mb->create_element( MBPOLYGON, foundIds, nP, polyNew );
380  mb->add_entities( outSet, &polyNew, 1 );
381 
382  // tag it with the index ids from tgt and src sets
383  int id = rs1.index( src ); // index starts from 0
384  mb->tag_set_data( srcParentTag, &polyNew, 1, &id );
385  id = rs2.index( tgt );
386  mb->tag_set_data( tgtParentTag, &polyNew, 1, &id );
387 
388  counting++;
389  mb->tag_set_data( countTag, &polyNew, 1, &counting );
390 
391 #ifdef ENABLE_DEBUG
392  if( dbg_1 )
393  {
394 
395  std::cout << "Count: " << counting + 1 << "\n";
396  std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
397  for( int i1 = 0; i1 < nP; i1++ )
398  std::cout << " " << mb->id_from_handle( foundIds[i1] );
399  std::cout << "\n";
400  std::vector< CartVect > posi( nP );
401  mb->get_coords( foundIds, nP, &( posi[0][0] ) );
402  for( int i1 = 0; i1 < nP; i1++ )
403  std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << posi[i1] << "\n";
404 
405  std::stringstream fff;
406  fff << "file0" << counting << ".vtk";
407  mb->write_mesh( fff.str().c_str(), &outSet, 1 );
408  }
409 #endif
410  }
411  delete[] foundIds;
412  foundIds = NULL;
413  return MB_SUCCESS;
414  // end copy
415 }

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.

20 {
21  // the points will be at most ?; they will describe a convex patch, after the points will be
22  // ordered and collapsed (eliminate doubles) the area is not really required get coordinates of
23  // the tgt quad
24  double cellArea = 0;
25  int num_nodes;
26  ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );
27  if( MB_SUCCESS != rval ) return 1.; // it should be an error
28 
29  nsTgt = num_nodes;
30 
31  rval = mb->get_coords( tgtConn, num_nodes, &( tgtCoords[0][0] ) );
32  if( MB_SUCCESS != rval ) return 1.; // it should be an error
33 
34  for( int j = 0; j < nsTgt; j++ )
35  {
36  // populate coords in the plane for intersection
37  // they should be oriented correctly, positively
38  tgtCoords2D[2 * j] = tgtCoords[j][0]; // x coordinate,
39  tgtCoords2D[2 * j + 1] = tgtCoords[j][1]; // y coordinate
40  }
41  for( int j = 1; j < nsTgt - 1; j++ )
42  cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
43 
44  return cellArea;
45 }

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: