Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::IntxRllCssphere Class Reference

#include <IntxRllCssphere.hpp>

+ Inheritance diagram for moab::IntxRllCssphere:
+ Collaboration diagram for moab::IntxRllCssphere:

Public Member Functions

 IntxRllCssphere (Interface *mbimpl)
 
virtual ~IntxRllCssphere ()
 
void set_radius (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)
 
- 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)
 

Private Attributes

double R
 
int plane
 
int srcEdgeType [4]
 

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 14 of file IntxRllCssphere.hpp.

Constructor & Destructor Documentation

◆ IntxRllCssphere()

moab::IntxRllCssphere::IntxRllCssphere ( Interface mbimpl)

Definition at line 15 of file IntxRllCssphere.cpp.

15 : Intx2Mesh( mbimpl ), R( 0.0 ), plane( 0 ) {}

◆ ~IntxRllCssphere()

moab::IntxRllCssphere::~IntxRllCssphere ( )
virtual

Definition at line 17 of file IntxRllCssphere.cpp.

17 {}

Member Function Documentation

◆ computeIntersectionBetweenTgtAndSrc()

ErrorCode moab::IntxRllCssphere::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 61 of file IntxRllCssphere.cpp.

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

References moab::IntxUtils::area2D(), moab::IntxUtils::borderPointsOfCSinRLL(), moab::IntxUtils::borderPointsOfXinY2(), moab::GeomUtil::bounding_boxes_overlap(), moab::Intx2Mesh::box_error, moab::IntxUtils::EdgeIntxRllCs(), 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, R, setup_tgt_cell(), moab::IntxUtils::SortAndRemoveDoubles2(), moab::Intx2Mesh::srcConn, moab::Intx2Mesh::srcCoords, moab::Intx2Mesh::srcCoords2D, srcEdgeType, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.

◆ findNodes()

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

Implements moab::Intx2Mesh.

Definition at line 208 of file IntxRllCssphere.cpp.

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

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, plane, R, moab::IntxUtils::reverse_gnomonic_projection(), 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().

◆ set_radius()

void moab::IntxRllCssphere::set_radius ( double  radius)
inline

Definition at line 21 of file IntxRllCssphere.hpp.

22  {
23  R = radius;
24  }

References R.

◆ setup_tgt_cell()

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

Implements moab::Intx2Mesh.

Definition at line 22 of file IntxRllCssphere.cpp.

23 {
24  // get coordinates of the tgt quad, to decide the gnomonic plane
25  double cellArea = 0;
26 
27  int num_nodes;
28  ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );
29 
30  if( MB_SUCCESS != rval ) return 1;
31  nsTgt = num_nodes;
32  // these edges will never be polygons, only quads or triangles
33 
34  // CartVect coords[4];
35  rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );
36  if( MB_SUCCESS != rval ) return 1;
37  CartVect middle = tgtCoords[0];
38  for( int i = 1; i < nsTgt; i++ )
39  middle += tgtCoords[i];
40  middle = 1. / nsTgt * middle;
41 
42  IntxUtils::decide_gnomonic_plane( middle, plane ); // output the plane
43  for( int j = 0; j < nsTgt; j++ )
44  {
45  // populate coords in the plane for intersection
46  // they should be oriented correctly, positively
47  int rc = IntxUtils::gnomonic_projection( tgtCoords[j], R, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );
48  if( rc != 0 ) return 1;
49  }
50 
51  for( int j = 1; j < nsTgt - 1; j++ )
52  cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
53 
54  // take tgt coords in order and compute area in plane
55  return cellArea;
56 }

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_SUCCESS, plane, R, moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.

Referenced by computeIntersectionBetweenTgtAndSrc().

Member Data Documentation

◆ plane

int moab::IntxRllCssphere::plane
private

◆ R

double moab::IntxRllCssphere::R
private

◆ srcEdgeType

int moab::IntxRllCssphere::srcEdgeType[4]
private

Definition at line 50 of file IntxRllCssphere.hpp.

Referenced by computeIntersectionBetweenTgtAndSrc().


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