Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
IntxRllCssphere.cpp
Go to the documentation of this file.
1 /*
2  * IntxRllCssphere.cpp
3  *
4  * Created on: Aug 8, 2014
5  * Author: iulian
6  */
7 
9 #include "moab/GeomUtil.hpp"
11 
12 namespace moab
13 {
14 
15 IntxRllCssphere::IntxRllCssphere( Interface* mbimpl ) : Intx2Mesh( mbimpl ), R( 0.0 ), plane( 0 ) {}
16 
18 
19 /*
20  * return also the area for robustness verification
21  */
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 }
57 
58 /* the elements are convex for sure, then do a gnomonic projection of both,
59  * compute intersection in the plane, then go back to the sphere for the points
60  * */
62  EntityHandle src,
63  double* P,
64  int& nP,
65  double& area,
66  int markb[MAXEDGES],
67  int markr[MAXEDGES],
68  int& nsSrc,
69  int& nsTgt,
70  bool check_boxes_first )
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 }
202 
203 // this method will also construct the triangles/polygons in the new mesh
204 // if we accept planar polygons, we just save them
205 // also, we could just create new vertices every time, and merge only in the end;
206 // could be too expensive, and the tolerance for merging could be an
207 // interesting topic
208 ErrorCode IntxRllCssphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP )
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 }
417 
418 } /* namespace moab */