Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
IntxUtils.cpp
Go to the documentation of this file.
1 /*
2  * IntxUtils.cpp
3  *
4  * Created on: Oct 3, 2012
5  */
6 #if defined( _MSC_VER ) || defined( WIN32 ) /* windows */
7 #define _USE_MATH_DEFINES // For M_PI
8 #endif
9 
10 #include <cmath>
11 #include <cassert>
12 #include <iostream>
13 
15 // this is from mbcoupler; maybe it should be moved somewhere in moab src
16 // right now, add a dependency to mbcoupler
17 // #include "ElemUtil.hpp"
18 #include "moab/MergeMesh.hpp"
19 #include "moab/ReadUtilIface.hpp"
20 #include "MBTagConventions.hpp"
21 
22 #define CHECKNEGATIVEAREA
23 #ifdef CHECKNEGATIVEAREA
24 #include <iomanip>
25 #endif
26 #include <queue>
27 #include <map>
28 
29 #ifdef MOAB_HAVE_TEMPESTREMAP
30 #include "GridElements.h"
31 #endif
32 
33 #ifdef MOAB_HAVE_EIGEN3
34 #define EIGEN_NO_DEBUG
35 #define EIGEN_MAX_CPP_VER 11
36 #include "Eigen/Dense"
37 #endif
38 
39 namespace moab
40 {
41 /**
42  * This code defines several utility functions for computing edge intersections and performing geometric operations.
43  *
44  * - `borderPointsOfXinY2`: Computes the border points of a set of points `X` inside another set of points `Y`.
45  * - `SortAndRemoveDoubles2`: Sorts a set of points `P` according to their angles and removes duplicate points.
46  * - `EdgeIntersections2`: Computes the intersections between the edges of two sets of points `blue` and `red`.
47  * - `EdgeIntxRllCs`: Computes the intersections between the edges of a set of points `blue` and a set of points `red` on a specific plane.
48  *
49  * The code also defines some helper structs and functions used by these utility functions.
50  */
51 
52 #define CORRTAGNAME "__correspondent"
53 #define MAXEDGES 10
54 
55 /**
56  * Computes the border points of X in Y2.
57  *
58  * @param X The array of points representing X.
59  * @param nX The number of points in X.
60  * @param Y The array of points representing Y.
61  * @param nY The number of points in Y.
62  * @param P The array to store the border points of X in Y2.
63  * @param side The array to store the side information for each point in X.
64  * @param epsilon_area The epsilon value for area comparison.
65  * @return The number of extra points found.
66  */
67 int IntxUtils::borderPointsOfXinY2( double* X, int nX, double* Y, int nY, double* P, int* side, double epsilon_area )
68 {
69  // 2 triangles, 3 corners, is the corner of X in Y?
70  // Y must have a positive area
71  /*
72  */
73  int extraPoint = 0;
74  for( int i = 0; i < nX; i++ )
75  {
76  // compute double the area of all nY triangles formed by a side of Y and a corner of X; if
77  // one is negative, stop (negative means it is outside; X and Y are all oriented such that
78  // they are positive oriented;
79  // if one area is negative, it means it is outside the convex region, for sure)
80  double* A = X + 2 * i;
81 
82  int inside = 1;
83  for( int j = 0; j < nY; j++ )
84  {
85  double* B = Y + 2 * j;
86  int j1 = ( j + 1 ) % nY;
87  double* C = Y + 2 * j1; // no copy of data
88 
89  double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
90  if( area2 < -epsilon_area )
91  {
92  inside = 0;
93  break;
94  }
95  }
96  if( inside )
97  {
98  side[i] = 1; // so vertex i of X is inside the convex region formed by Y
99  // so side has nX dimension (first array)
100  P[extraPoint * 2] = A[0];
101  P[extraPoint * 2 + 1] = A[1];
102  extraPoint++;
103  }
104  }
105  return extraPoint;
106 }
107 
108 // used to order according to angle, so it can remove double easily
110 {
111  double angle;
112  int index;
113 };
114 
116 {
117  return lhs.angle < rhs.angle;
118 }
119 
120 /**
121  * Sorts and removes duplicate points in the given array.
122  *
123  * Note: nP might be modified too, we will remove duplicates if found
124  *
125  * @param P The array of points to be sorted and checked for duplicates.
126  * @param nP The number of points in P.
127  * @param epsilon_1 The epsilon value for distance comparison.
128  * @return 0 if successful.
129  */
130 int IntxUtils::SortAndRemoveDoubles2( double* P, int& nP, double epsilon_1 )
131 {
132  if( nP < 2 ) return 0; // nothing to do
133 
134  // center of gravity for the points
135  double c[2] = { 0., 0. };
136  int k = 0;
137  for( k = 0; k < nP; k++ )
138  {
139  c[0] += P[2 * k];
140  c[1] += P[2 * k + 1];
141  }
142  c[0] /= nP;
143  c[1] /= nP;
144 
145  // how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
146  // intersection points
147  struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
148 
149  for( k = 0; k < nP; k++ )
150  {
151  double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
152  if( x != 0. || y != 0. )
153  {
154  pairAngleIndex[k].angle = atan2( y, x );
155  }
156  else
157  {
158  pairAngleIndex[k].angle = 0;
159  // it would mean that the cells are touching at a vertex
160  }
161  pairAngleIndex[k].index = k;
162  }
163 
164  // this should be faster than the bubble sort we had before
165  std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
166  // copy now to a new double array
167  double PCopy[10 * MAXEDGES]; // the same dimension as P; very conservative, but faster than
168  // reallocate for a vector
169  for( k = 0; k < nP; k++ ) // index will show where it should go now;
170  {
171  int ck = pairAngleIndex[k].index;
172  PCopy[2 * k] = P[2 * ck];
173  PCopy[2 * k + 1] = P[2 * ck + 1];
174  }
175  // now copy from PCopy over original P location
176  std::copy( PCopy, PCopy + 2 * nP, P );
177 
178  // eliminate duplicates, finally
179 
180  int i = 0, j = 1; // the next one; j may advance faster than i
181  // check the unit
182  // double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
183  // than 1.e-5 cm
184  while( j < nP )
185  {
186  double d2 = dist2( &P[2 * i], &P[2 * j] );
187  if( d2 > epsilon_1 )
188  {
189  i++;
190  P[2 * i] = P[2 * j];
191  P[2 * i + 1] = P[2 * j + 1];
192  }
193  j++;
194  }
195  // test also the last point with the first one (index 0)
196  // the first one could be at -PI; last one could be at +PI, according to atan2 span
197 
198  double d2 = dist2( P, &P[2 * i] ); // check the first and last points (ordered from -pi to +pi)
199  if( d2 > epsilon_1 )
200  {
201  nP = i + 1;
202  }
203  else
204  nP = i; // effectively delete the last point (that would have been the same with first)
205  if( nP == 0 ) nP = 1; // we should be left with at least one point we already tested if nP is 0 originally
206  return 0;
207 }
208 
209 /**
210  * Computes the edge intersections of two elements.
211  *
212  * @param blue The array of points representing the blue element.
213  * @param nsBlue The number of points in the blue element.
214  * @param red The array of points representing the red element.
215  * @param nsRed The number of points in the red element.
216  * @param markb The array to mark the intersecting edges of the blue element.
217  * @param markr The array to mark the intersecting edges of the red element.
218  * @param points The array to store the intersection points.
219  * @param nPoints The number of intersection points found.
220  * @return The error code.
221  */
223  int nsBlue,
224  double* red,
225  int nsRed,
226  int* markb,
227  int* markr,
228  double* points,
229  int& nPoints )
230 {
231  /* EDGEINTERSECTIONS computes edge intersections of two elements
232  [P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
233  and blue ( stored column wise )
234  (point coordinates are stored column-wise, in counter clock
235  order) the points P where their edges intersect. In addition,
236  in n the indices of which neighbors of red are also intersecting
237  with blue are given.
238  */
239 
240  // points is an array with enough slots (24 * 2 doubles)
241  nPoints = 0;
242  for( int i = 0; i < MAXEDGES; i++ )
243  {
244  markb[i] = markr[i] = 0;
245  }
246 
247  for( int i = 0; i < nsBlue; i++ )
248  {
249  for( int j = 0; j < nsRed; j++ )
250  {
251  double b[2];
252  double a[2][2]; // 2*2
253  int iPlus1 = ( i + 1 ) % nsBlue;
254  int jPlus1 = ( j + 1 ) % nsRed;
255  for( int k = 0; k < 2; k++ )
256  {
257  b[k] = red[2 * j + k] - blue[2 * i + k];
258  // row k of a: a(k, 0), a(k, 1)
259  a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
260  a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
261  }
262  double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
263  if( fabs( delta ) > 1.e-14 ) // this is close to machine epsilon
264  {
265  // not parallel
266  double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
267  double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
268  if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
269  {
270  // the intersection is good
271  for( int k = 0; k < 2; k++ )
272  {
273  points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
274  }
275  markb[i] = 1; // so neighbor number i of blue will be considered too.
276  markr[j] = 1; // this will be used in advancing red around blue quad
277  nPoints++;
278  }
279  }
280  // the case delta ~ 0. will be considered by the interior points logic
281  }
282  }
283  return MB_SUCCESS;
284 }
285 
286 /**
287  * Computes the edge intersections between a RLL and CS quad.
288  *
289  * Note: Special function.
290  *
291  * @param blue The array of points representing the blue element.
292  * @param bluec The array of Cartesian coordinates for the blue element.
293  * @param blueEdgeType The array of edge types for the blue element.
294  * @param nsBlue The number of points in the blue element.
295  * @param red The array of points representing the red element.
296  * @param redc The array of Cartesian coordinates for the red element.
297  * @param nsRed The number of points in the red element.
298  * @param markb The array to mark the intersecting edges of the blue element.
299  * @param markr The array to mark the intersecting edges of the red element.
300  * @param plane The plane of intersection.
301  * @param R The radius of the sphere.
302  * @param points The array to store the intersection points.
303  * @param nPoints The number of intersection points found.
304  * @return The error code.
305  */
307  CartVect* bluec,
308  int* blueEdgeType,
309  int nsBlue,
310  double* red,
311  CartVect* redc,
312  int nsRed,
313  int* markb,
314  int* markr,
315  int plane,
316  double R,
317  double* points,
318  int& nPoints )
319 {
320  // if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
321  // everything else the same (except if there are 2 points resulting, which is rare)
322  for( int i = 0; i < 4; i++ )
323  { // always at most 4 , so maybe don't bother
324  markb[i] = markr[i] = 0;
325  }
326 
327  for( int i = 0; i < nsBlue; i++ )
328  {
329  int iPlus1 = ( i + 1 ) % nsBlue;
330  if( blueEdgeType[i] == 0 ) // old style, just 2d
331  {
332  for( int j = 0; j < nsRed; j++ )
333  {
334  double b[2];
335  double a[2][2]; // 2*2
336 
337  int jPlus1 = ( j + 1 ) % nsRed;
338  for( int k = 0; k < 2; k++ )
339  {
340  b[k] = red[2 * j + k] - blue[2 * i + k];
341  // row k of a: a(k, 0), a(k, 1)
342  a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
343  a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
344  }
345  double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
346  if( fabs( delta ) > 1.e-14 )
347  {
348  // not parallel
349  double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
350  double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
351  if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
352  {
353  // the intersection is good
354  for( int k = 0; k < 2; k++ )
355  {
356  points[2 * nPoints + k] =
357  blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
358  }
359  markb[i] = 1; // so neighbor number i of blue will be considered too.
360  markr[j] = 1; // this will be used in advancing red around blue quad
361  nPoints++;
362  }
363  } // if the edges are too "parallel", skip them
364  }
365  }
366  else // edge type is 1, so use 3d intersection
367  {
368  CartVect& C = bluec[i];
369  CartVect& D = bluec[iPlus1];
370  for( int j = 0; j < nsRed; j++ )
371  {
372  int jPlus1 = ( j + 1 ) % nsRed; // nsRed is just 4, forget about it, usually
373  CartVect& A = redc[j];
374  CartVect& B = redc[jPlus1];
375  int np = 0;
376  double E[9];
377  intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
378  if( np == 0 ) continue;
379  if( np >= 2 )
380  {
381  std::cout << "intersection with 2 points :" << A << B << C << D << "\n";
382  }
383  for( int k = 0; k < np; k++ )
384  {
385  gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
386  points[2 * nPoints + 1] );
387  nPoints++;
388  }
389  markb[i] = 1; // so neighbor number i of blue will be considered too.
390  markr[j] = 1; // this will be used in advancing red around blue quad
391  }
392  }
393  }
394  return MB_SUCCESS;
395 }
396 
397 // vec utils related to gnomonic projection on a sphere
398 
399 // vec utils
400 
401 /*
402  *
403  * position on a sphere of radius R
404  * if plane specified, use it; if not, return the plane, and the point in the plane
405  * there are 6 planes, numbered 1 to 6
406  * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
407  *
408  * projection on the plane will preserve the orientation, such that a triangle, quad pointing
409  * outside the sphere will have a positive orientation in the projection plane
410  */
411 void IntxUtils::decide_gnomonic_plane( const CartVect& pos, int& plane )
412 {
413  // decide plane, based on max x, y, z
414  if( fabs( pos[0] ) < fabs( pos[1] ) )
415  {
416  if( fabs( pos[2] ) < fabs( pos[1] ) )
417  {
418  // pos[1] is biggest
419  if( pos[1] > 0 )
420  plane = 2;
421  else
422  plane = 4;
423  }
424  else
425  {
426  // pos[2] is biggest
427  if( pos[2] < 0 )
428  plane = 5;
429  else
430  plane = 6;
431  }
432  }
433  else
434  {
435  if( fabs( pos[2] ) < fabs( pos[0] ) )
436  {
437  // pos[0] is the greatest
438  if( pos[0] > 0 )
439  plane = 1;
440  else
441  plane = 3;
442  }
443  else
444  {
445  // pos[2] is biggest
446  if( pos[2] < 0 )
447  plane = 5;
448  else
449  plane = 6;
450  }
451  }
452  return;
453 }
454 
455 // point on a sphere is projected on one of six planes, decided earlier
456 ErrorCode IntxUtils::gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 )
457 {
458  double alfa = 1.; // the new point will be on line alfa*pos
459 
460  switch( plane )
461  {
462  case 1: {
463  // the plane with x = R; x>y, x>z
464  // c1->y, c2->z
465  alfa = R / pos[0];
466  c1 = alfa * pos[1];
467  c2 = alfa * pos[2];
468  break;
469  }
470  case 2: {
471  // y = R -> zx
472  alfa = R / pos[1];
473  c1 = alfa * pos[2];
474  c2 = alfa * pos[0];
475  break;
476  }
477  case 3: {
478  // x=-R, -> yz
479  alfa = -R / pos[0];
480  c1 = -alfa * pos[1]; // the sign is to preserve orientation
481  c2 = alfa * pos[2];
482  break;
483  }
484  case 4: {
485  // y = -R
486  alfa = -R / pos[1];
487  c1 = -alfa * pos[2]; // the sign is to preserve orientation
488  c2 = alfa * pos[0];
489  break;
490  }
491  case 5: {
492  // z = -R
493  alfa = -R / pos[2];
494  c1 = -alfa * pos[0]; // the sign is to preserve orientation
495  c2 = alfa * pos[1];
496  break;
497  }
498  case 6: {
499  alfa = R / pos[2];
500  c1 = alfa * pos[0];
501  c2 = alfa * pos[1];
502  break;
503  }
504  default:
505  return MB_FAILURE; // error
506  }
507 
508  return MB_SUCCESS; // no error
509 }
510 
511 // given the position on plane (one out of 6), find out the position on sphere
513  const double& c2,
514  double R,
515  int plane,
516  CartVect& pos )
517 {
518 
519  // the new point will be on line beta*pos
520  double len = sqrt( c1 * c1 + c2 * c2 + R * R );
521  double beta = R / len; // it is less than 1, in general
522 
523  switch( plane )
524  {
525  case 1: {
526  // the plane with x = R; x>y, x>z
527  // c1->y, c2->z
528  pos[0] = beta * R;
529  pos[1] = c1 * beta;
530  pos[2] = c2 * beta;
531  break;
532  }
533  case 2: {
534  // y = R -> zx
535  pos[1] = R * beta;
536  pos[2] = c1 * beta;
537  pos[0] = c2 * beta;
538  break;
539  }
540  case 3: {
541  // x=-R, -> yz
542  pos[0] = -R * beta;
543  pos[1] = -c1 * beta; // the sign is to preserve orientation
544  pos[2] = c2 * beta;
545  break;
546  }
547  case 4: {
548  // y = -R
549  pos[1] = -R * beta;
550  pos[2] = -c1 * beta; // the sign is to preserve orientation
551  pos[0] = c2 * beta;
552  break;
553  }
554  case 5: {
555  // z = -R
556  pos[2] = -R * beta;
557  pos[0] = -c1 * beta; // the sign is to preserve orientation
558  pos[1] = c2 * beta;
559  break;
560  }
561  case 6: {
562  pos[2] = R * beta;
563  pos[0] = c1 * beta;
564  pos[1] = c2 * beta;
565  break;
566  }
567  default:
568  return MB_FAILURE; // error
569  }
570 
571  return MB_SUCCESS; // no error
572 }
573 
574 void IntxUtils::gnomonic_unroll( double& c1, double& c2, double R, int plane )
575 {
576  double tmp;
577  switch( plane )
578  {
579  case 1:
580  break; // nothing
581  case 2: // rotate + 90
582  tmp = c1;
583  c1 = -c2;
584  c2 = tmp;
585  c1 = c1 + 2 * R;
586  break;
587  case 3:
588  c1 = c1 + 4 * R;
589  break;
590  case 4: // rotate with -90 x-> -y; y -> x
591 
592  tmp = c1;
593  c1 = c2;
594  c2 = -tmp;
595  c1 = c1 - 2 * R;
596  break;
597  case 5: // South Pole
598  // rotate 180 then move to (-2, -2)
599  c1 = -c1 - 2. * R;
600  c2 = -c2 - 2. * R;
601  break;
602  ;
603  case 6: // North Pole
604  c1 = c1 - 2. * R;
605  c2 = c2 + 2. * R;
606  break;
607  }
608  return;
609 }
610 // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
612  EntityHandle inSet,
613  double R,
614  bool centers_only,
615  EntityHandle& outSet )
616 {
617  std::string parTagName( "PARALLEL_PARTITION" );
618  Tag part_tag;
619  Range partSets;
620  ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
621  if( MB_SUCCESS == rval && part_tag != 0 )
622  {
623  rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
624  }
625  rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
626  // Get all entities of dimension 2
627  Range inputRange; // get
628  rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
629  rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
630 
631  std::map< EntityHandle, int > partsAssign;
632  std::map< int, EntityHandle > newPartSets;
633  if( !partSets.empty() )
634  {
635  // get all cells, and assign parts
636  for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
637  {
638  EntityHandle pSet = *setIt;
639  Range ents;
640  rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
641  int val;
642  rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
643  // create a new set with the same part id tag, in the outSet
644  EntityHandle newPartSet;
645  rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
646  rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
647  newPartSets[val] = newPartSet;
648  rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
649  for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
650  {
651  partsAssign[*it] = val;
652  }
653  }
654  }
655 
656  if( centers_only )
657  {
658  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
659  {
661  EntityHandle cell = *it;
662  rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
663  int plane;
664  decide_gnomonic_plane( center, plane );
665  double c[3];
666  c[2] = 0.;
667  gnomonic_projection( center, R, plane, c[0], c[1] );
668 
669  gnomonic_unroll( c[0], c[1], R, plane );
670 
672  rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
673  rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
674  }
675  }
676  else
677  {
678  // distribute the cells to 6 planes, based on the center
679  Range subranges[6];
680  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
681  {
683  EntityHandle cell = *it;
684  rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
685  int plane;
686  decide_gnomonic_plane( center, plane );
687  subranges[plane - 1].insert( cell );
688  }
689  for( int i = 1; i <= 6; i++ )
690  {
691  Range verts;
692  rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
693  std::map< EntityHandle, EntityHandle > corr;
694  for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
695  {
696  CartVect vect;
697  EntityHandle v = *vt;
698  rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
699  double c[3];
700  c[2] = 0.;
701  gnomonic_projection( vect, R, i, c[0], c[1] );
702  gnomonic_unroll( c[0], c[1], R, i );
704  rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
705  corr[v] = vertex; // for new connectivity
706  }
707  EntityHandle new_conn[20]; // max edges in 2d ?
708  for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
709  {
710  EntityHandle eh = *eit;
711  const EntityHandle* conn = NULL;
712  int num_nodes;
713  rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
714  // build a new vertex array
715  for( int j = 0; j < num_nodes; j++ )
716  new_conn[j] = corr[conn[j]];
717  EntityType type = mb->type_from_handle( eh );
718  EntityHandle newCell;
719  rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
720  rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
721  std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
722  if( mit != partsAssign.end() )
723  {
724  int val = mit->second;
725  rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
726  }
727  }
728  }
729  }
730 
731  return MB_SUCCESS;
732 }
733 void IntxUtils::transform_coordinates( double* avg_position, int projection_type )
734 {
735  if( projection_type == 1 )
736  {
737  double R =
738  avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
739  R = sqrt( R );
740  double lat = asin( avg_position[2] / R );
741  double lon = atan2( avg_position[1], avg_position[0] );
742  avg_position[0] = lon;
743  avg_position[1] = lat;
744  avg_position[2] = R;
745  }
746  else if( projection_type == 2 ) // gnomonic projection
747  {
748  CartVect pos( avg_position );
749  int gplane;
750  IntxUtils::decide_gnomonic_plane( pos, gplane );
751 
752  IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
753  avg_position[2] = 0;
754  IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
755  }
756 }
757 /*
758  *
759  use physical_constants, only : dd_pi
760  type(cartesian3D_t), intent(in) :: cart
761  type(spherical_polar_t) :: sphere
762 
763  sphere%r=distance(cart)
764  sphere%lat=ASIN(cart%z/sphere%r)
765  sphere%lon=0
766 
767  ! ==========================================================
768  ! enforce three facts:
769  !
770  ! 1) lon at poles is defined to be zero
771  !
772  ! 2) Grid points must be separated by about .01 Meter (on earth)
773  ! from pole to be considered "not the pole".
774  !
775  ! 3) range of lon is { 0<= lon < 2*pi }
776  !
777  ! ==========================================================
778 
779  if (distance(cart) >= DIST_THRESHOLD) then
780  sphere%lon=ATAN2(cart%y,cart%x)
781  if (sphere%lon<0) then
782  sphere%lon=sphere%lon+2*DD_PI
783  end if
784  end if
785 
786  end function cart_to_spherical
787  */
789 {
790  SphereCoords res;
791  res.R = cart3d.length();
792  if( res.R < 0 )
793  {
794  res.lon = res.lat = 0.;
795  return res;
796  }
797  res.lat = asin( cart3d[2] / res.R );
798  res.lon = atan2( cart3d[1], cart3d[0] );
799  if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
800  // there are some defines it depends on :(
801  // #if defined __USE_BSD || defined __USE_XOPEN ???
802 
803  return res;
804 }
805 
806 /*
807  * ! ===================================================================
808  ! spherical_to_cart:
809  ! converts spherical polar {lon,lat} to 3D cartesian {x,y,z}
810  ! on unit sphere
811  ! ===================================================================
812 
813  function spherical_to_cart(sphere) result (cart)
814 
815  type(spherical_polar_t), intent(in) :: sphere
816  type(cartesian3D_t) :: cart
817 
818  cart%x=sphere%r*COS(sphere%lat)*COS(sphere%lon)
819  cart%y=sphere%r*COS(sphere%lat)*SIN(sphere%lon)
820  cart%z=sphere%r*SIN(sphere%lat)
821 
822  end function spherical_to_cart
823  */
825 {
826  CartVect res;
827  res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
828  res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
829  res[2] = sc.R * sin( sc.lat ); // z
830  return res;
831 }
832 
834 {
835  Range nodes;
836  ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
837  if( rval != moab::MB_SUCCESS ) return rval;
838 
839  // one by one, get the node and project it on the sphere, with a radius given
840  // the center of the sphere is at 0,0,0
841  for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
842  {
843  EntityHandle nd = *nit;
844  CartVect pos;
845  rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
846  if( rval != moab::MB_SUCCESS ) return rval;
847  double len = pos.length();
848  if( len == 0. ) return MB_FAILURE;
849  pos = R / len * pos;
850  rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
851  if( rval != moab::MB_SUCCESS ) return rval;
852  }
853  return MB_SUCCESS;
854 }
855 
856 // assume they are one the same sphere
857 double IntxAreaUtils::spherical_angle( const double* A, const double* B, const double* C, double Radius )
858 {
859  // the angle by definition is between the planes OAB and OBC
860  CartVect a( A );
861  CartVect b( B );
862  CartVect c( C );
863  double err1 = a.length_squared() - Radius * Radius;
864  if( fabs( err1 ) > 0.0001 )
865  {
866  std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n";
867  }
868  CartVect normalOAB = a * b;
869  CartVect normalOCB = c * b;
870  return angle( normalOAB, normalOCB );
871 }
872 
873 // could be bigger than M_PI;
874 // angle at B could be bigger than M_PI, if the orientation is such that ABC points toward the
875 // interior
876 double IntxUtils::oriented_spherical_angle( const double* A, const double* B, const double* C )
877 {
878  // assume the same radius, sphere at origin
879  CartVect a( A ), b( B ), c( C );
880  CartVect normalOAB = a * b;
881  CartVect normalOCB = c * b;
882  CartVect orient = ( c - b ) * ( a - b );
883  double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
884  if( ang != ang )
885  {
886  // signal of a nan
887  std::cout << a << " " << b << " " << c << "\n";
888  std::cout << ang << "\n";
889  }
890  if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
891 
892  return ang;
893 }
894 
895 double IntxAreaUtils::area_spherical_triangle( const double* A, const double* B, const double* C, double Radius )
896 {
897  switch( m_eAreaMethod )
898  {
899  case Girard:
900  return area_spherical_triangle_girard( A, B, C, Radius );
901 #ifdef MOAB_HAVE_TEMPESTREMAP
902  case GaussQuadrature:
903  return area_spherical_triangle_GQ( A, B, C );
904 #endif
905  case lHuiller:
906  default:
907  return area_spherical_triangle_lHuiller( A, B, C, Radius );
908  }
909 }
910 
911 double IntxAreaUtils::area_spherical_polygon( const double* A, int N, double Radius, int* sign )
912 {
913  switch( m_eAreaMethod )
914  {
915  case Girard:
916  return area_spherical_polygon_girard( A, N, Radius );
917 #ifdef MOAB_HAVE_TEMPESTREMAP
918  case GaussQuadrature:
919  return area_spherical_polygon_GQ( A, N ) * Radius * Radius; //area_spherical_polygon_GQ normalizes
920 #endif
921  case lHuiller:
922  default:
923  return area_spherical_polygon_lHuiller( A, N, Radius, sign );
924  }
925 }
926 
927 double IntxAreaUtils::area_spherical_triangle_girard( const double* A, const double* B, const double* C, double Radius )
928 {
929  double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
930  spherical_angle( C, A, B, Radius ) - M_PI;
931  double area = Radius * Radius * correction;
932  // now, is it negative or positive? is it pointing toward the center or outward?
933  CartVect a( A ), b( B ), c( C );
934  CartVect abc = ( b - a ) * ( c - a );
935  if( abc % a > 0 ) // dot product positive, means ABC points out
936  return area;
937  else
938  return -area;
939 }
940 
941 double IntxAreaUtils::area_spherical_polygon_girard( const double* A, int N, double Radius )
942 {
943  // this should work for non-convex polygons too
944  // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
945  //
946  if( N <= 2 ) return 0.;
947  double sum_angles = 0.;
948  for( int i = 0; i < N; i++ )
949  {
950  int i1 = ( i + 1 ) % N;
951  int i2 = ( i + 2 ) % N;
952  sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
953  }
954  double correction = sum_angles - ( N - 2 ) * M_PI;
955  return Radius * Radius * correction;
956 }
957 
958 double IntxAreaUtils::area_spherical_polygon_lHuiller( const double* A, int N, double Radius, int* sign )
959 {
960  // This should work for non-convex polygons too
961  // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
962  // We also assume that the orientation is positive;
963  // If negative orientation, the area will be negative
964  if( N <= 2 ) return 0.;
965 
966  int lsign = 1; // assume positive orientain
967  double area = 0.;
968  for( int i = 1; i < N - 1; i++ )
969  {
970  int i1 = i + 1;
971  double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
972  if( areaTriangle < 0 )
973  lsign = -1; // signal that we have at least one triangle with negative orientation ;
974  // possible nonconvex polygon
975  area += areaTriangle;
976  }
977  if( sign ) *sign = lsign;
978 
979  return area;
980 }
981 
982 #ifdef MOAB_HAVE_TEMPESTREMAP
983 
984 double IntxAreaUtils::area_spherical_polygon_GQ( const double* A, int N )
985 {
986  // this should work for non-convex polygons too
987  // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
988  // We also assume that the orientation is positive;
989  // If negative orientation, the area can be negative
990  if( N <= 2 ) return 0.;
991 
992  // assume positive orientation
993  double area = 0.;
994  for( int i = 1; i < N - 1; i++ )
995  {
996  area += area_spherical_triangle_GQ( A, A + 3 * i, A + 3 * ( i + 1 ) );
997  }
998  return area;
999 }
1000 
1001 template < typename Derived >
1002 Eigen::Array< typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime > shift(
1003  const Eigen::ArrayBase< Derived >& array,
1004  int positions )
1005 {
1006  Eigen::Array< typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime > result = array;
1007  if( positions > 0 )
1008  {
1009  result.segment( positions, array.size() - positions ) = array.head( array.size() - positions );
1010  result.head( positions ).setZero();
1011  }
1012  else if( positions < 0 )
1013  {
1014  result.head( array.size() + positions ) = array.tail( array.size() + positions );
1015  result.tail( -positions ).setZero();
1016  }
1017  return result;
1018 }
1019 
1020 double IntxAreaUtils::area_spherical_triangle_GQ( const double* inode1, const double* inode2, const double* inode3 )
1021 {
1022 #if defined( MOAB_HAVE_EIGEN3 )
1023  typedef Eigen::Map< const Eigen::Vector3d > V3d;
1024  const V3d node1( inode1 );
1025  const V3d node2( inode2 );
1026  const V3d node3( inode3 );
1027  const int nOrder = 6;
1028 
1029  // If we change the quadrature order, use the call: GaussQuadrature::GetPoints(nOrder, 0.0, 1.0, dG, dW);
1030  const double dG[6] = { 0.03376524289842397, 0.1693953067668678, 0.3806904069584016,
1031  0.6193095930415985, 0.8306046932331322, 0.966234757101576 };
1032  const double dW[6] = { 0.08566224618958521, 0.1803807865240693, 0.2339569672863455,
1033  0.2339569672863455, 0.1803807865240693, 0.08566224618958521 };
1034 
1035  double dFaceArea = 0.0;
1036  Eigen::Vector3d dF, dF2, dDaF, dDbF, dDaG, dDbG;
1037  double nodeCross[3];
1038 
1039  // Calculate area at quadrature node and sum it up
1040  for( int p = 0; p < nOrder; p++ )
1041  {
1042  for( int q = 0; q < nOrder; q++ )
1043  {
1044 
1045  const double dA = dG[p];
1046  const double dB = dG[q];
1047 
1048  // V3d dF = (1.0 - dB) * (((1.0 - dA) * node1) + (dA * node2)) + (dB * node3);
1049  dF = ( ( ( 1.0 - dB ) * ( 1.0 - dA ) ) * node1 ) + ( ( ( 1.0 - dB ) * dA ) * node2 ) + ( dB * node3 );
1050  dF2 = dF.array().square();
1051 
1052  dDaF = ( node1 - node2 );
1053 
1054  // dDbF = -( 1.0 - dA ) * node1 - dA * node2 + node3;
1055  dDbF = ( node3 - node1 ) + dA * dDaF;
1056  dDaF *= ( dB - 1.0 );
1057 
1058  const double dDenomTerm = std::pow( dF.norm(), -3.0 );
1059 
1060  // Eigen::Vector3d temp1 = dF2;
1061  // temp1( 2 ) += dF2( 0 ); // temp1 = [dF2(0), dF2(1), dF2(0)+dF2(2)]
1062  // Eigen::Vector3d temp2 = dDaF.cwiseProduct( dF ); // temp2 = [dDaF(0)*dF(0), dDaF(1)*dF(1), dDaF(2)*dF(2)]
1063  // dDaG = dDaF.cwiseProduct( temp1.segment( 1, 2 ) ) - temp2.segment( 1, 2 );
1064  // Eigen::Vector3d temp3 = dDbF.cwiseProduct( dF ); // temp3 = [dDbF(0)*dF(0), dDbF(1)*dF(1), dDbF(2)*dF(2)]
1065  // dDbG = dDbF.cwiseProduct( temp1.segment( 1, 2 ) ) - temp3.segment( 1, 2 );
1066 
1067  // Eigen::Vector3d dF2_shifted = dF2 + Eigen::Vector3d( dF2( 1 ), dF2( 2 ), dF2( 0 ) );
1068  // Eigen::Vector3d dDaF_shifted = Eigen::Vector3d( dDaF( 1 ), dDaF( 2 ), dDaF( 0 ) ).cwiseProduct( dF );
1069 
1070  // dDaG = dDaF.cwiseProduct( dF2_shifted ) - dF.cwiseProduct( dDaF_shifted );
1071 
1072  // Eigen::Vector3d dDbF_shifted = Eigen::Vector3d( dDbF( 1 ), dDbF( 2 ), dDbF( 0 ) ).cwiseProduct( dF );
1073 
1074  // dDbG = dDbF.cwiseProduct( dF2_shifted ) - dF.cwiseProduct( dDbF_shifted );
1075 
1076  dDaG( 0 ) = dDaF( 0 ) * ( dF2( 1 ) + dF2( 2 ) ) - dF( 0 ) * ( dDaF( 1 ) * dF( 1 ) + dDaF( 2 ) * dF( 2 ) );
1077  dDaG( 1 ) = dDaF( 1 ) * ( dF2( 0 ) + dF2( 2 ) ) - dF( 1 ) * ( dDaF( 0 ) * dF( 0 ) + dDaF( 2 ) * dF( 2 ) );
1078  dDaG( 2 ) = dDaF( 2 ) * ( dF2( 0 ) + dF2( 1 ) ) - dF( 2 ) * ( dDaF( 0 ) * dF( 0 ) + dDaF( 1 ) * dF( 1 ) );
1079 
1080  // dDbG = dDbF.cwiseProduct( ( dF2.array().shift( -1 ) + dF2.array().shift( -2 ) ) ) + dF.cwiseProduct( dDbF );
1081  dDbG( 0 ) = dDbF( 0 ) * ( dF2( 1 ) + dF2( 2 ) ) - dF( 0 ) * ( dDbF( 1 ) * dF( 1 ) + dDbF( 2 ) * dF( 2 ) );
1082  dDbG( 1 ) = dDbF( 1 ) * ( dF2( 0 ) + dF2( 2 ) ) - dF( 1 ) * ( dDbF( 0 ) * dF( 0 ) + dDbF( 2 ) * dF( 2 ) );
1083  dDbG( 2 ) = dDbF( 2 ) * ( dF2( 0 ) + dF2( 1 ) ) - dF( 2 ) * ( dDbF( 0 ) * dF( 0 ) + dDbF( 1 ) * dF( 1 ) );
1084 
1085  // Scale the vectors by the denominator
1086  dDaG *= dDenomTerm;
1087  dDbG *= dDenomTerm;
1088 
1089  // Cross product gives local Jacobian: dGaG x dDbG
1090  nodeCross[0] = dDaG( 1 ) * dDbG( 2 ) - dDaG( 2 ) * dDbG( 1 );
1091  nodeCross[1] = dDaG( 2 ) * dDbG( 0 ) - dDaG( 0 ) * dDbG( 2 );
1092  nodeCross[2] = dDaG( 0 ) * dDbG( 1 ) - dDaG( 1 ) * dDbG( 0 );
1093 
1094  const double dJacobian =
1095  std::sqrt( nodeCross[0] * nodeCross[0] + nodeCross[1] * nodeCross[1] + nodeCross[2] * nodeCross[2] );
1096 
1097  // dFaceArea += 2.0 * dW[p] * dW[q] * (1.0 - dG[q]) * dJacobian;
1098  dFaceArea += dW[p] * dW[q] * dJacobian;
1099  }
1100  }
1101 
1102  return dFaceArea;
1103 #else
1104  /* compute the area by using Gauss-Quadratures; use TR interfaces directly */
1105  Face face( 3 );
1106  NodeVector nodes( 3 );
1107  nodes[0] = Node( inode1[0], inode1[1], inode1[2] );
1108  nodes[1] = Node( inode2[0], inode2[1], inode2[2] );
1109  nodes[2] = Node( inode3[0], inode3[1], inode3[2] );
1110  face.SetNode( 0, 0 );
1111  face.SetNode( 1, 1 );
1112  face.SetNode( 2, 2 );
1113  return CalculateFaceArea( face, nodes );
1114 #endif
1115 }
1116 
1117 #endif
1118 
1119 /*
1120  * l'Huiller's formula for spherical triangle
1121  * http://williams.best.vwh.net/avform.htm
1122  * a, b, c are arc measures in radians, too
1123  * A, B, C are angles on the sphere, for which we already have formula
1124  * c
1125  * A -------B
1126  * \ |
1127  * \ |
1128  * \b |a
1129  * \ |
1130  * \ |
1131  * \ |
1132  * \C|
1133  * \|
1134  *
1135  * (The angle at B is not necessarily a right angle)
1136  *
1137  * sin(a) sin(b) sin(c)
1138  * ----- = ------ = ------
1139  * sin(A) sin(B) sin(C)
1140  *
1141  * In terms of the sides (this is excess, as before, but numerically stable)
1142  *
1143  * E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
1144  */
1146  const double* ptB,
1147  const double* ptC,
1148  double Radius )
1149 {
1150 
1151  // now, a is angle BOC, O is origin
1152  CartVect vA( ptA ), vB( ptB ), vC( ptC );
1153  double a = angle_robust( vB, vC );
1154  double b = angle_robust( vC, vA );
1155  double c = angle_robust( vA, vB );
1156  int sign = 1;
1157  // if( fabs( ( vA * vB ) % vC ) < 1e-17 ) sign = -1;
1158  if( ( vA * vB ) % vC < 0 ) sign = -1;
1159  double s = ( a + b + c ) / 2;
1160  double a1 = ( s - a ) / 2;
1161  double b1 = ( s - b ) / 2;
1162  double c1 = ( s - c ) / 2;
1163 #ifdef MOAB_HAVE_TEMPESTREMAP
1164  if( fabs( a1 ) < 1.e-14 || fabs( b1 ) < 1.e-14 || fabs( c1 ) < 1.e-14 )
1165  {
1166  double area = area_spherical_triangle_GQ( ptA, ptB, ptC ) * sign;
1167 #ifdef VERBOSE
1168  std::cout << " very obtuse angle, use TR to compute area " << " a1:" << a1 << " b1:" << b1 << " c1:" << c1
1169  << "\n";
1170  std::cout << " area with TR: " << area << "\n";
1171 #endif
1172  return area;
1173  }
1174 #endif
1175  double tmp = tan( s / 2 ) * tan( a1 ) * tan( b1 ) * tan( c1 );
1176  if( tmp < 0. ) tmp = 0.;
1177 
1178  double E = 4 * atan( sqrt( tmp ) );
1179  if( E != E ) std::cout << " NaN at spherical triangle area \n";
1180 
1181  double area = sign * E * Radius * Radius;
1182 
1183 #ifdef CHECKNEGATIVEAREA
1184  if( area < 0 )
1185  {
1186  std::cout << "negative area: " << area << "\n";
1187  std::cout << std::setprecision( 15 );
1188  std::cout << "vA: " << vA << "\n";
1189  std::cout << "vB: " << vB << "\n";
1190  std::cout << "vC: " << vC << "\n";
1191  std::cout << "sign: " << sign << "\n";
1192  std::cout << " a: " << a << "\n";
1193  std::cout << " b: " << b << "\n";
1194  std::cout << " c: " << c << "\n";
1195  }
1196 #endif
1197 
1198  return area;
1199 }
1200 
1202 {
1203  // Get all entities of dimension 2
1204  Range inputRange;
1205  ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1206 
1207  // Filter by elements that are owned by current process
1208  std::vector< int > ownerinfo( inputRange.size(), -1 );
1209  Tag intxOwnerTag;
1210  rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
1211  if( MB_SUCCESS == rval )
1212  {
1213  rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1214  }
1215 
1216  // compare total area with 4*M_PI * R^2
1217  int ie = 0;
1218  double total_area = 0.;
1219  for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
1220  {
1221 
1222  // All zero/positive owner data represents ghosted elems
1223  if( ownerinfo[ie++] >= 0 ) continue;
1224 
1225  EntityHandle eh = *eit;
1226  const double elem_area = this->area_spherical_element( mb, eh, R );
1227 
1228  // check whether the area of the spherical element is positive.
1229  if( elem_area <= 0 )
1230  {
1231  std::cout << "Area of element " << mb->id_from_handle( eh ) << " is = " << elem_area << "\n";
1232  mb->list_entity( eh );
1233  }
1234  assert( elem_area > 0 );
1235 
1236  // sum up the contribution
1237  total_area += elem_area;
1238  }
1239 
1240  // return total mesh area
1241  return total_area;
1242 }
1243 
1245 {
1246  // get the nodes, then the coordinates
1247  const EntityHandle* verts;
1248  int nsides;
1249  ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1250 
1251  // account for possible padded polygons
1252  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1253  nsides--;
1254 
1255  // get coordinates
1256  std::vector< double > coords( 3 * nsides );
1257  rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1258 
1259  // compute and return the area of the polygonal element
1260  return area_spherical_polygon( &coords[0], nsides, R );
1261 }
1262 
1264 {
1265  SphereCoords sph1 = cart_to_spherical( p1 );
1266  SphereCoords sph2 = cart_to_spherical( p2 );
1267  // radius should be the same
1268  return sph1.R *
1269  acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
1270 }
1271 
1272 //
1273 /**
1274  * @brief Enforces convexity for a given set of polygons.
1275  *
1276  * This function checks each polygon in the input set and computes the angles of each vertex.
1277  * If a reflex angle is found, the polygon is broken into triangles and added back to the set.
1278  * This process continues until all polygons in the set are convex.
1279  *
1280  * @param mb The interface to the MOAB instance.
1281  * @param lset The handle of the input set containing the polygons.
1282  * @param my_rank The rank of the local process.
1283  * @return The error code indicating the success or failure of the operation.
1284  */
1286 {
1287  Range inputRange;
1288  ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
1289 
1290  Tag corrTag = 0;
1291  EntityHandle dumH = 0;
1292  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
1293  if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
1294 
1295  Tag gidTag;
1296  rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
1297 
1298  std::vector< double > coords;
1299  coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
1300  // we should create a queue with new polygons that need processing for reflex angles
1301  // (obtuse)
1302  std::queue< EntityHandle > newPolys;
1303  int brokenPolys = 0;
1304  Range::iterator eit = inputRange.begin();
1305  while( eit != inputRange.end() || !newPolys.empty() )
1306  {
1307  EntityHandle eh;
1308  if( eit != inputRange.end() )
1309  {
1310  eh = *eit;
1311  ++eit;
1312  }
1313  else
1314  {
1315  eh = newPolys.front();
1316  newPolys.pop();
1317  }
1318  // get the nodes, then the coordinates
1319  const EntityHandle* verts;
1320  int num_nodes;
1321  rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
1322  int nsides = num_nodes;
1323  // account for possible padded polygons
1324  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1325  nsides--;
1326  EntityHandle corrHandle = 0;
1327  if( corrTag )
1328  {
1329  rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
1330  }
1331  int gid = 0;
1332  rval = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
1333  coords.resize( 3 * nsides );
1334  if( nsides < 4 ) continue; // if already triangles, don't bother
1335  // get coordinates
1336  rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
1337  // compute each angle
1338  bool alreadyBroken = false;
1339 
1340  for( int i = 0; i < nsides; i++ )
1341  {
1342  double* A = &coords[3 * i];
1343  double* B = &coords[3 * ( ( i + 1 ) % nsides )];
1344  double* C = &coords[3 * ( ( i + 2 ) % nsides )];
1345  double angle = IntxUtils::oriented_spherical_angle( A, B, C );
1346  if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
1347  {
1348  if( alreadyBroken )
1349  {
1350  mb->list_entities( &eh, 1 );
1351  mb->list_entities( verts, nsides );
1352  double* D = &coords[3 * ( ( i + 3 ) % nsides )];
1353  std::cout << "ABC: " << angle << " \n";
1354  std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
1355  std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
1356  std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
1357  std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
1358 
1359  return MB_FAILURE;
1360  }
1361  // the bad angle is at i+1;
1362  // create 1 triangle and one polygon; add the polygon to the input range, so
1363  // it will be processed too
1364  // also, add both to the set :) and remove the original polygon from the set
1365  // break the next triangle, even though not optimal
1366  // so create the triangle i+1, i+2, i+3; remove i+2 from original list
1367  // even though not optimal in general, it is good enough.
1368  EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
1369  verts[( i + 3 ) % nsides] };
1370  // create a polygon with num_nodes-1 vertices, and connectivity
1371  // verts[i+1], verts[i+3], (all except i+2)
1372  std::vector< EntityHandle > conn( nsides - 1 );
1373  for( int j = 1; j < nsides; j++ )
1374  {
1375  conn[j - 1] = verts[( i + j + 2 ) % nsides];
1376  }
1377  EntityHandle newElement;
1378  rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
1379 
1380  rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
1381  if( corrTag )
1382  {
1383  rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
1384  }
1385  rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
1386  if( nsides == 4 )
1387  {
1388  // create another triangle
1389  rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
1390  }
1391  else
1392  {
1393  // create another polygon, and add it to the inputRange
1394  rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
1395  newPolys.push( newElement ); // because it has less number of edges, the
1396  // reverse should work to find it.
1397  }
1398  rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
1399  if( corrTag )
1400  {
1401  rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
1402  }
1403  rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
1404  rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
1405  brokenPolys++;
1406  alreadyBroken = true; // get out of the loop, element is broken
1407  }
1408  }
1409  }
1410  if( brokenPolys > 0 )
1411  {
1412  std::cout << "on local process " << my_rank << ", " << brokenPolys
1413  << " concave polygons were decomposed in convex ones \n";
1414 #ifdef VERBOSE
1415  std::stringstream fff;
1416  fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
1417  rval = mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 );MB_CHK_ERR( rval );
1418  std::cout << "wrote new file set: " << fff.str() << "\n";
1419 #endif
1420  }
1421  return MB_SUCCESS;
1422 }
1423 
1424 // looking at quad connectivity, collapse to triangle if 2 nodes equal
1425 // then delete the old quad
1427 {
1428  Range quads;
1429  ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
1430  Tag gid;
1431  gid = mb->globalId_tag();
1432  for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
1433  {
1434  EntityHandle quad = *qit;
1435  const EntityHandle* conn4 = NULL;
1436  int num_nodes = 0;
1437  rval = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
1438  for( int i = 0; i < num_nodes; i++ )
1439  {
1440  int next_node_index = ( i + 1 ) % num_nodes;
1441  if( conn4[i] == conn4[next_node_index] )
1442  {
1443  // form a triangle and delete the quad
1444  // first get the global id, to set it on triangle later
1445  int global_id = 0;
1446  rval = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
1447  int i2 = ( i + 2 ) % num_nodes;
1448  int i3 = ( i + 3 ) % num_nodes;
1449  EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
1450  EntityHandle tri;
1451  rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
1452  mb->add_entities( set, &tri, 1 );
1453  mb->remove_entities( set, &quad, 1 );
1454  mb->delete_entities( &quad, 1 );
1455  rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
1456  }
1457  }
1458  }
1459  return MB_SUCCESS;
1460 }
1461 
1463 {
1464  Range cells2d;
1465  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );MB_CHK_ERR( rval );
1466  for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
1467  {
1468  EntityHandle cell = *qit;
1469  const EntityHandle* conn = NULL;
1470  int num_nodes = 0;
1471  rval = mb->get_connectivity( cell, conn, num_nodes );MB_CHK_ERR( rval );
1472  if( num_nodes < 3 ) return MB_FAILURE;
1473 
1474  double coords[9];
1475  rval = mb->get_coords( conn, 3, coords );MB_CHK_ERR( rval );
1476 
1477  double area;
1478  if( R > 0 )
1479  area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
1480  else
1481  area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
1482  if( area < 0 )
1483  {
1484  // compute all area, do not revert if total area is positive
1485  std::vector< double > coords2( 3 * num_nodes );
1486  // get coordinates
1487  rval = mb->get_coords( conn, num_nodes, &coords2[0] );MB_CHK_ERR( rval );
1488  double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
1489  if( totArea < 0 )
1490  {
1491  std::vector< EntityHandle > newconn( num_nodes );
1492  for( int i = 0; i < num_nodes; i++ )
1493  {
1494  newconn[num_nodes - 1 - i] = conn[i];
1495  }
1496  rval = mb->set_connectivity( cell, &newconn[0], num_nodes );MB_CHK_ERR( rval );
1497  }
1498  else
1499  {
1500  std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
1501  }
1502  }
1503  }
1504  return MB_SUCCESS;
1505 }
1506 
1507 // distance along a great circle on a sphere of radius 1
1508 // page 4
1509 double IntxUtils::distance_on_sphere( double la1, double te1, double la2, double te2 )
1510 {
1511  return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
1512 }
1513 
1514 /*
1515  * given 2 great circle arcs, AB and CD, compute the unique intersection point, if it exists
1516  * in between
1517  */
1518 ErrorCode IntxUtils::intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E )
1519 {
1520  // first verify A, B, C, D are on the same sphere
1521  double R2 = R * R;
1522  const double Tolerance = 1.e-12 * R2;
1523 
1524  CartVect a( A ), b( B ), c( C ), d( D );
1525 
1526  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1527  fabs( d.length_squared() - R2 ) >
1528  10 * Tolerance )
1529  return MB_FAILURE;
1530 
1531  CartVect n1 = a * b;
1532  if( n1.length_squared() < Tolerance ) return MB_FAILURE;
1533 
1534  CartVect n2 = c * d;
1535  if( n2.length_squared() < Tolerance ) return MB_FAILURE;
1536  CartVect n3 = n1 * n2;
1537  n3.normalize();
1538 
1539  n3 = R * n3;
1540  // the intersection is either n3 or -n3
1541  CartVect n4 = a * n3, n5 = n3 * b;
1542  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1543  {
1544  // n3 is good for ab, see if it is good for cd
1545  n4 = c * n3;
1546  n5 = n3 * d;
1547  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1548  {
1549  E[0] = n3[0];
1550  E[1] = n3[1];
1551  E[2] = n3[2];
1552  }
1553  else
1554  return MB_FAILURE;
1555  }
1556  else
1557  {
1558  // try -n3
1559  n3 = -n3;
1560  n4 = a * n3, n5 = n3 * b;
1561  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1562  {
1563  // n3 is good for ab, see if it is good for cd
1564  n4 = c * n3;
1565  n5 = n3 * d;
1566  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1567  {
1568  E[0] = n3[0];
1569  E[1] = n3[1];
1570  E[2] = n3[2];
1571  }
1572  else
1573  return MB_FAILURE;
1574  }
1575  else
1576  return MB_FAILURE;
1577  }
1578 
1579  return MB_SUCCESS;
1580 }
1581 
1582 // verify that result is in between a and b on a great circle arc, and between c and d on a constant
1583 // latitude arc
1584 static bool verify( CartVect a, CartVect b, CartVect c, CartVect d, double x, double y, double z )
1585 {
1586  // to check, the point has to be between a and b on a great arc, and between c and d on a const
1587  // lat circle
1588  CartVect s( x, y, z );
1589  CartVect n1 = a * b;
1590  CartVect n2 = a * s;
1591  CartVect n3 = s * b;
1592  if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
1593 
1594  // do the same for c, d, s, in plane z=0
1595  c[2] = d[2] = s[2] = 0.; // bring everything in the same plane, z=0;
1596 
1597  n1 = c * d;
1598  n2 = c * s;
1599  n3 = s * d;
1600  if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
1601 
1602  return true;
1603 }
1604 
1606  double* B,
1607  double* C,
1608  double* D,
1609  double R,
1610  double* E,
1611  int& np )
1612 {
1613  const double distTol = R * 1.e-6;
1614  const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
1615  np = 0; // number of points in intersection
1616  CartVect a( A ), b( B ), c( C ), d( D );
1617  // check input first
1618  double R2 = R * R;
1619  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1620  fabs( d.length_squared() - R2 ) >
1621  10 * Tolerance )
1622  return MB_FAILURE;
1623 
1624  if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
1625  if( ( c - d ).length_squared() < Tolerance ) // edges are too short
1626  return MB_FAILURE;
1627 
1628  // CD is the const latitude arc
1629  if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
1630  return MB_FAILURE;
1631 
1632  if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
1633 
1634  // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
1635  // great circle arc AB normal to the AB circle:
1636  CartVect n1 = a * b; // the normal to the great circle arc (circle)
1637  // solve the system of equations:
1638  /*
1639  * n1%(x, y, z) = 0 // on the great circle
1640  * z = C[2];
1641  * x^2+y^2+z^2 = R^2
1642  */
1643  double z = C[2];
1644  if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
1645  {
1646  // it is the Equator; check if the const lat edge is Equator too
1647  if( fabs( C[2] ) > distTol )
1648  {
1649  return MB_FAILURE; // no intx, too far from Eq
1650  }
1651  else
1652  {
1653  // all points are on the equator
1654  //
1655  CartVect cd = c * d;
1656  // by convention, c<d, positive is from c to d
1657  // is a or b between c , d?
1658  CartVect ca = c * a;
1659  CartVect ad = a * d;
1660  CartVect cb = c * b;
1661  CartVect bd = b * d;
1662  bool agtc = ( ca % cd >= -Tolerance ); // a>c?
1663  bool dgta = ( ad % cd >= -Tolerance ); // d>a?
1664  bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
1665  bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
1666  if( agtc )
1667  {
1668  if( dgta )
1669  {
1670  // a is for sure a point
1671  E[0] = a[0];
1672  E[1] = a[1];
1673  E[2] = a[2];
1674  np++;
1675  if( bgtc )
1676  {
1677  if( dgtb )
1678  {
1679  // b is also in between c and d
1680  E[3] = b[0];
1681  E[4] = b[1];
1682  E[5] = b[2];
1683  np++;
1684  }
1685  else
1686  {
1687  // then order is c a d b, intx is ad
1688  E[3] = d[0];
1689  E[4] = d[1];
1690  E[5] = d[2];
1691  np++;
1692  }
1693  }
1694  else
1695  {
1696  // b is less than c, so b c a d, intx is ac
1697  E[3] = c[0];
1698  E[4] = c[1];
1699  E[5] = c[2];
1700  np++; // what if E[0] is E[3]?
1701  }
1702  }
1703  else // c < d < a
1704  {
1705  if( dgtb ) // d is for sure in
1706  {
1707  E[0] = d[0];
1708  E[1] = d[1];
1709  E[2] = d[2];
1710  np++;
1711  if( bgtc ) // c<b<d<a
1712  {
1713  // another point is b
1714  E[3] = b[0];
1715  E[4] = b[1];
1716  E[5] = b[2];
1717  np++;
1718  }
1719  else // b<c<d<a
1720  {
1721  // another point is c
1722  E[3] = c[0];
1723  E[4] = c[1];
1724  E[5] = c[2];
1725  np++;
1726  }
1727  }
1728  else
1729  {
1730  // nothing, order is c, d < a, b
1731  }
1732  }
1733  }
1734  else // a < c < d
1735  {
1736  if( bgtc )
1737  {
1738  // c is for sure in
1739  E[0] = c[0];
1740  E[1] = c[1];
1741  E[2] = c[2];
1742  np++;
1743  if( dgtb )
1744  {
1745  // a < c < b < d; second point is b
1746  E[3] = b[0];
1747  E[4] = b[1];
1748  E[5] = b[2];
1749  np++;
1750  }
1751  else
1752  {
1753  // a < c < d < b; second point is d
1754  E[3] = d[0];
1755  E[4] = d[1];
1756  E[5] = d[2];
1757  np++;
1758  }
1759  }
1760  else // a, b < c < d
1761  {
1762  // nothing
1763  }
1764  }
1765  }
1766  // for the 2 points selected, see if it is only one?
1767  // no problem, maybe it will be collapsed later anyway
1768  if( np > 0 ) return MB_SUCCESS;
1769  return MB_FAILURE; // no intersection
1770  }
1771  {
1772  if( fabs( n1[0] ) <= fabs( n1[1] ) )
1773  {
1774  // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
1775  // (u+v*x)^2+x^2=R2-z^2
1776  // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
1777  // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
1778  // x1,2 =
1779  double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
1780  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1781  double delta = b1 * b1 - 4 * a1 * c1;
1782  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1783  if( delta > Tolerance ) // 2 solutions possible
1784  {
1785  double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1786  double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1787  double y1 = u + v * x1;
1788  double y2 = u + v * x2;
1789  if( verify( a, b, c, d, x1, y1, z ) )
1790  {
1791  E[0] = x1;
1792  E[1] = y1;
1793  E[2] = z;
1794  np++;
1795  }
1796  if( verify( a, b, c, d, x2, y2, z ) )
1797  {
1798  E[3 * np + 0] = x2;
1799  E[3 * np + 1] = y2;
1800  E[3 * np + 2] = z;
1801  np++;
1802  }
1803  }
1804  else
1805  {
1806  // one solution
1807  double x1 = -b1 / 2 / a1;
1808  double y1 = u + v * x1;
1809  if( verify( a, b, c, d, x1, y1, z ) )
1810  {
1811  E[0] = x1;
1812  E[1] = y1;
1813  E[2] = z;
1814  np++;
1815  }
1816  }
1817  }
1818  else
1819  {
1820  // resolve eq in y, reverse
1821  // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
1822  // (u+v*y)^2+y^2 -R2+z^2 =0
1823  // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
1824  //
1825  // x1,2 =
1826  double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
1827  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1828  double delta = b1 * b1 - 4 * a1 * c1;
1829  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1830  if( delta > Tolerance ) // 2 solutions possible
1831  {
1832  double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1833  double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1834  double x1 = u + v * y1;
1835  double x2 = u + v * y2;
1836  if( verify( a, b, c, d, x1, y1, z ) )
1837  {
1838  E[0] = x1;
1839  E[1] = y1;
1840  E[2] = z;
1841  np++;
1842  }
1843  if( verify( a, b, c, d, x2, y2, z ) )
1844  {
1845  E[3 * np + 0] = x2;
1846  E[3 * np + 1] = y2;
1847  E[3 * np + 2] = z;
1848  np++;
1849  }
1850  }
1851  else
1852  {
1853  // one solution
1854  double y1 = -b1 / 2 / a1;
1855  double x1 = u + v * y1;
1856  if( verify( a, b, c, d, x1, y1, z ) )
1857  {
1858  E[0] = x1;
1859  E[1] = y1;
1860  E[2] = z;
1861  np++;
1862  }
1863  }
1864  }
1865  }
1866 
1867  if( np <= 0 ) return MB_FAILURE;
1868  return MB_SUCCESS;
1869 }
1870 
1871 #if 0
1872 ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1)
1873 {
1874  Range cells;
1875  ErrorCode rval = mb->get_entities_by_dimension(sf1, 2, cells);
1876  if (MB_SUCCESS!= rval)
1877  return rval;
1878  Range edges;
1879  rval = mb->get_adjacencies(cells, 1, true, edges, Interface::UNION);
1880  if (MB_SUCCESS!= rval)
1881  return rval;
1882 
1883  Tag edgeTypeTag;
1884  int default_int=0;
1885  rval = mb->tag_get_handle("edge_type", 1, MB_TYPE_INTEGER, edgeTypeTag,
1886  MB_TAG_DENSE | MB_TAG_CREAT, &default_int);
1887  if (MB_SUCCESS!= rval)
1888  return rval;
1889  // add edges to the set? not yet, maybe later
1890  // if edge horizontal, set value to 1
1891  int type_constant_lat=1;
1892  for (Range::iterator eit=edges.begin(); eit!=edges.end(); ++eit)
1893  {
1894  EntityHandle edge = *eit;
1895  const EntityHandle *conn=0;
1896  int num_n=0;
1897  rval = mb->get_connectivity(edge, conn, num_n );
1898  if (MB_SUCCESS!= rval)
1899  return rval;
1900  double coords[6];
1901  rval = mb->get_coords(conn, 2, coords);
1902  if (MB_SUCCESS!= rval)
1903  return rval;
1904  if (fabs( coords[2]-coords[5] )< 1.e-6 )
1905  {
1906  rval = mb->tag_set_data(edgeTypeTag, &edge, 1, &type_constant_lat);
1907  if (MB_SUCCESS!= rval)
1908  return rval;
1909  }
1910  }
1911 
1912  return MB_SUCCESS;
1913 }
1914 #endif
1915 
1916 // decide in a different metric if the corners of CS quad are
1917 // in the interior of an RLL quad
1919  double* red2dc,
1920  int nsRed,
1921  CartVect* bluec,
1922  int nsBlue,
1923  int* blueEdgeType,
1924  double* P,
1925  int* side,
1926  double epsil )
1927 {
1928  int extraPoints = 0;
1929  // first decide the blue z coordinates
1930  CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
1931  for( int i = 0; i < nsBlue; i++ )
1932  {
1933  if( blueEdgeType[i] == 0 )
1934  {
1935  int iP1 = ( i + 1 ) % nsBlue;
1936  if( bluec[i][2] > bluec[iP1][2] )
1937  {
1938  A = bluec[i];
1939  B = bluec[iP1];
1940  C = bluec[( i + 2 ) % nsBlue];
1941  D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
1942  break;
1943  }
1944  }
1945  }
1946  if( nsBlue == 3 && B[2] < 0 )
1947  {
1948  // select D to be C
1949  D = C;
1950  C = B; // B is the south pole then
1951  }
1952  // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
1953  // AB is const longitude, BC and DA constant latitude
1954  // check now each of the red points if they are inside this rectangle
1955  for( int i = 0; i < nsRed; i++ )
1956  {
1957  CartVect& X = redc[i];
1958  if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
1959  // now decide if it is between the planes OAB and OCD
1960  if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
1961  {
1962  side[i] = 1; //
1963  // it means point X is in the rectangle that we want , on the sphere
1964  // pass the coords 2d
1965  P[extraPoints * 2] = red2dc[2 * i];
1966  P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
1967  extraPoints++;
1968  }
1969  }
1970  return extraPoints;
1971 }
1972 
1974 {
1975  ReadUtilIface* read_iface;
1976  ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
1977  // create the handle tag for the corresponding element / vertex
1978 
1979  EntityHandle dum = 0;
1980  Tag corrTag = 0; // it will be created here
1981  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
1982 
1983  // give the same global id to new verts and cells created in the lagr(departure) mesh
1984  Tag gid = mb->globalId_tag();
1985 
1986  Range quads;
1987  rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
1988 
1989  Range connecVerts;
1990  rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
1991 
1992  std::map< EntityHandle, EntityHandle > newNodes;
1993 
1994  std::vector< double* > coords;
1995  EntityHandle start_vert, start_elem, *connect;
1996  int num_verts = connecVerts.size();
1997  rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
1998  if( MB_SUCCESS != rval ) return rval;
1999  // fill it up
2000  int i = 0;
2001  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
2002  {
2003  EntityHandle oldV = *vit;
2004  CartVect posi;
2005  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
2006 
2007  int global_id;
2008  rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
2009  EntityHandle new_vert = start_vert + i;
2010  // Cppcheck warning (false positive): variable coords is assigned a value that is never used
2011  coords[0][i] = posi[0];
2012  coords[1][i] = posi[1];
2013  coords[2][i] = posi[2];
2014 
2015  newNodes[oldV] = new_vert;
2016  // set also the correspondent tag :)
2017  rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
2018 
2019  // also the other side
2020  // need to check if we really need this; the new vertex will never need the old vertex
2021  // we have the global id which is the same
2022  rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
2023  // set the global id on the corresponding vertex the same as the initial vertex
2024  rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
2025  }
2026  // now create new quads in order (in a sequence)
2027 
2028  rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
2029  if( MB_SUCCESS != rval ) return rval;
2030  int ie = 0;
2031  for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
2032  {
2033  EntityHandle q = *it;
2034  int nnodes;
2035  const EntityHandle* conn;
2036  rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
2037  int global_id;
2038  rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
2039 
2040  for( int ii = 0; ii < nnodes; ii++ )
2041  {
2042  EntityHandle v1 = conn[ii];
2043  connect[4 * ie + ii] = newNodes[v1];
2044  }
2045  EntityHandle newElement = start_elem + ie;
2046 
2047  // set the corresponding tag; not sure we need this one, from old to new
2048  rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
2049  rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
2050 
2051  // set the global id
2052  rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
2053 
2054  rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
2055  }
2056 
2057  rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
2058 
2059  return MB_SUCCESS;
2060 }
2061 
2063  EntityHandle file_set,
2064  double merge_tol,
2065  std::vector< Tag >& tagList )
2066 {
2067  Range verts;
2068  ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
2069  rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
2070 
2071  MergeMesh mm( mb );
2072 
2073  // remove the vertices from the set, before merging
2074 
2075  rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
2076 
2077  // now correct vertices that are repeated in polygons
2078  rval = remove_padded_vertices( mb, file_set, tagList );
2079  return MB_SUCCESS;
2080 }
2081 
2082 ErrorCode IntxUtils::remove_padded_vertices( Interface* mb, EntityHandle file_set, std::vector< Tag >& tagList )
2083 {
2084 
2085  // now correct vertices that are repeated in polygons
2086  Range cells;
2087  ErrorCode rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
2088 
2089  Range verts;
2090  rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
2091 
2092  Range modifiedCells; // will be deleted at the end; keep the gid
2093  Range newCells;
2094 
2095  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
2096  {
2097  EntityHandle cell = *cit;
2098  const EntityHandle* connec = NULL;
2099  int num_verts = 0;
2100  rval = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
2101 
2102  std::vector< EntityHandle > newConnec;
2103  newConnec.push_back( connec[0] ); // at least one vertex
2104  int index = 0;
2105  int new_size = 1;
2106  while( index < num_verts - 2 )
2107  {
2108  int next_index = ( index + 1 );
2109  if( connec[next_index] != newConnec[new_size - 1] )
2110  {
2111  newConnec.push_back( connec[next_index] );
2112  new_size++;
2113  }
2114  index++;
2115  }
2116  // add the last one only if different from previous and first node
2117  if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
2118  {
2119  newConnec.push_back( connec[num_verts - 1] );
2120  new_size++;
2121  }
2122  if( new_size < num_verts )
2123  {
2124  // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
2125  modifiedCells.insert( cell );
2126  // create a new cell with type triangle, quad or polygon
2127  EntityType type = MBTRI;
2128  if( new_size == 3 )
2129  type = MBTRI;
2130  else if( new_size == 4 )
2131  type = MBQUAD;
2132  else if( new_size > 4 )
2133  type = MBPOLYGON;
2134 
2135  // create new cell
2136  EntityHandle newCell;
2137  rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
2138  // set the old id to the new element
2139  newCells.insert( newCell );
2140  double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
2141  for( size_t i = 0; i < tagList.size(); i++ )
2142  {
2143  rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
2144  rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
2145  }
2146  }
2147  }
2148 
2149  rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
2150  rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
2151  rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
2152  rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
2153 
2154  return MB_SUCCESS;
2155 }
2156 
2157 } // namespace moab