Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
Classes | Static Public Member Functions | List of all members
moab::IntxUtils Class Reference

#include <IntxUtils.hpp>

Classes

struct  SphereCoords
 

Static Public Member Functions

static double dist2 (double *a, double *b)
 
static double area2D (double *a, double *b, double *c)
 
static int borderPointsOfXinY2 (double *X, int nX, double *Y, int nY, double *P, int *side, double epsilon_area)
 
static int SortAndRemoveDoubles2 (double *P, int &nP, double epsilon)
 
static ErrorCode EdgeIntersections2 (double *blue, int nsBlue, double *red, int nsRed, int *markb, int *markr, double *points, int &nPoints)
 
static ErrorCode EdgeIntxRllCs (double *blue, CartVect *bluec, int *blueEdgeType, int nsBlue, double *red, CartVect *redc, int nsRed, int *markb, int *markr, int plane, double Radius, double *points, int &nPoints)
 
static void decide_gnomonic_plane (const CartVect &pos, int &oPlane)
 
static ErrorCode gnomonic_projection (const CartVect &pos, double R, int plane, double &c1, double &c2)
 
static ErrorCode gnomonic_projection_generalized (const CartVect &pos, const CartVect axis[3], double &c1, double &c2)
 
static ErrorCode global_gnomonic_projection_general (Interface *mb, EntityHandle inSet, CartVect P, EntityHandle &outSet)
 
static ErrorCode reverse_gnomonic_projection (const double &c1, const double &c2, double R, int plane, CartVect &pos)
 
static void gnomonic_unroll (double &c1, double &c2, double R, int plane)
 
static ErrorCode global_gnomonic_projection (Interface *mb, EntityHandle inSet, double R, bool centers_only, EntityHandle &outSet)
 
static ErrorCode gnomonic_projection_plane_at_point (CartVect P, CartVect &u, CartVect &v)
 
static void transform_coordinates (double *avg_position, int projection_type)
 
static SphereCoords cart_to_spherical (CartVect &)
 
static CartVect spherical_to_cart (SphereCoords &)
 
static ErrorCode ScaleToRadius (Interface *mb, EntityHandle set, double R)
 
static double distance_on_great_circle (CartVect &p1, CartVect &p2)
 
static ErrorCode enforce_convexity (Interface *mb, EntityHandle set, int rank=0)
 Enforces convexity for a given set of polygons. More...
 
static double oriented_spherical_angle (const double *A, const double *B, const double *C)
 
static ErrorCode fix_degenerate_quads (Interface *mb, EntityHandle set)
 
static double distance_on_sphere (double la1, double te1, double la2, double te2)
 
static ErrorCode intersect_great_circle_arcs (double *A, double *B, double *C, double *D, double R, double *E)
 
static ErrorCode intersect_great_circle_arc_with_clat_arc (double *A, double *B, double *C, double *D, double R, double *E, int &np)
 
static int borderPointsOfCSinRLL (CartVect *redc, double *red2dc, int nsRed, CartVect *bluec, int nsBlue, int *blueEdgeType, double *P, int *side, double epsil)
 
static ErrorCode deep_copy_set_with_quads (Interface *mb, EntityHandle source_set, EntityHandle dest_set)
 
static ErrorCode remove_duplicate_vertices (Interface *mb, EntityHandle file_set, double merge_tol, std::vector< Tag > &tagList)
 
static ErrorCode remove_padded_vertices (Interface *mb, EntityHandle file_set, std::vector< Tag > &tagList)
 
static ErrorCode max_diagonal (Interface *mb, Range cells, int max_edges, double &diagonal)
 

Detailed Description

Definition at line 15 of file IntxUtils.hpp.

Member Function Documentation

◆ area2D()

static double moab::IntxUtils::area2D ( double *  a,
double *  b,
double *  c 
)
inlinestatic

Definition at line 25 of file IntxUtils.hpp.

26  {
27  // (b-a)x(c-a) / 2
28  return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
29  }

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshEdges::EdgeSplits(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), moab::IntxAreaUtils::positive_orientation(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), and moab::IntxRllCssphere::setup_tgt_cell().

◆ borderPointsOfCSinRLL()

int moab::IntxUtils::borderPointsOfCSinRLL ( CartVect *  redc,
double *  red2dc,
int  nsRed,
CartVect *  bluec,
int  nsBlue,
int *  blueEdgeType,
double *  P,
int *  side,
double  epsil 
)
static

Definition at line 2200 of file IntxUtils.cpp.

2209 {
2210  int extraPoints = 0;
2211  // first decide the blue z coordinates
2212  CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
2213  for( int i = 0; i < nsBlue; i++ )
2214  {
2215  if( blueEdgeType[i] == 0 )
2216  {
2217  int iP1 = ( i + 1 ) % nsBlue;
2218  if( bluec[i][2] > bluec[iP1][2] )
2219  {
2220  A = bluec[i];
2221  B = bluec[iP1];
2222  C = bluec[( i + 2 ) % nsBlue];
2223  D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
2224  break;
2225  }
2226  }
2227  }
2228  if( nsBlue == 3 && B[2] < 0 )
2229  {
2230  // select D to be C
2231  D = C;
2232  C = B; // B is the south pole then
2233  }
2234  // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
2235  // AB is const longitude, BC and DA constant latitude
2236  // check now each of the red points if they are inside this rectangle
2237  for( int i = 0; i < nsRed; i++ )
2238  {
2239  CartVect& X = redc[i];
2240  if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
2241  // now decide if it is between the planes OAB and OCD
2242  if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
2243  {
2244  side[i] = 1; //
2245  // it means point X is in the rectangle that we want , on the sphere
2246  // pass the coords 2d
2247  P[extraPoints * 2] = red2dc[2 * i];
2248  P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
2249  extraPoints++;
2250  }
2251  }
2252  return extraPoints;
2253 }

Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

◆ borderPointsOfXinY2()

int moab::IntxUtils::borderPointsOfXinY2 ( double *  X,
int  nX,
double *  Y,
int  nY,
double *  P,
int *  side,
double  epsilon_area 
)
static

Computes the border points of X in Y2.

Parameters
XThe array of points representing X.
nXThe number of points in X.
YThe array of points representing Y.
nYThe number of points in Y.
PThe array to store the border points of X in Y2.
sideThe array to store the side information for each point in X.
epsilon_areaThe epsilon value for area comparison.
Returns
The number of extra points found.

Definition at line 66 of file IntxUtils.cpp.

67 {
68  // 2 triangles, 3 corners, is the corner of X in Y?
69  // Y must have a positive area
70  /*
71  */
72  int extraPoint = 0;
73  for( int i = 0; i < nX; i++ )
74  {
75  // compute double the area of all nY triangles formed by a side of Y and a corner of X; if
76  // one is negative, stop (negative means it is outside; X and Y are all oriented such that
77  // they are positive oriented;
78  // if one area is negative, it means it is outside the convex region, for sure)
79  double* A = X + 2 * i;
80 
81  int inside = 1;
82  for( int j = 0; j < nY; j++ )
83  {
84  double* B = Y + 2 * j;
85  int j1 = ( j + 1 ) % nY;
86  double* C = Y + 2 * j1; // no copy of data
87 
88  double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
89  if( area2 < -epsilon_area )
90  {
91  inside = 0;
92  break;
93  }
94  }
95  if( inside )
96  {
97  side[i] = 1; // so vertex i of X is inside the convex region formed by Y
98  // so side has nX dimension (first array)
99  P[extraPoint * 2] = A[0];
100  P[extraPoint * 2 + 1] = A[1];
101  extraPoint++;
102  }
103  }
104  return extraPoint;
105 }

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), and moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

◆ cart_to_spherical()

IntxUtils::SphereCoords moab::IntxUtils::cart_to_spherical ( CartVect &  cart3d)
static

Definition at line 1072 of file IntxUtils.cpp.

1073 {
1074  SphereCoords res;
1075  res.R = cart3d.length();
1076  if( res.R < 0 )
1077  {
1078  res.lon = res.lat = 0.;
1079  return res;
1080  }
1081  res.lat = asin( cart3d[2] / res.R );
1082  res.lon = atan2( cart3d[1], cart3d[0] );
1083  if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
1084  // there are some defines it depends on :(
1085  // #if defined __USE_BSD || defined __USE_XOPEN ???
1086 
1087  return res;
1088 }

References moab::IntxUtils::SphereCoords::lat, moab::CartVect::length(), moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

Referenced by distance_on_great_circle(), and main().

◆ decide_gnomonic_plane()

void moab::IntxUtils::decide_gnomonic_plane ( const CartVect &  pos,
int &  oPlane 
)
static

Definition at line 410 of file IntxUtils.cpp.

411 {
412  // decide plane, based on max x, y, z
413  if( fabs( pos[0] ) < fabs( pos[1] ) )
414  {
415  if( fabs( pos[2] ) < fabs( pos[1] ) )
416  {
417  // pos[1] is biggest
418  if( pos[1] > 0 )
419  plane = 2;
420  else
421  plane = 4;
422  }
423  else
424  {
425  // pos[2] is biggest
426  if( pos[2] < 0 )
427  plane = 5;
428  else
429  plane = 6;
430  }
431  }
432  else
433  {
434  if( fabs( pos[2] ) < fabs( pos[0] ) )
435  {
436  // pos[0] is the greatest
437  if( pos[0] > 0 )
438  plane = 1;
439  else
440  plane = 3;
441  }
442  else
443  {
444  // pos[2] is biggest
445  if( pos[2] < 0 )
446  plane = 5;
447  else
448  plane = 6;
449  }
450  }
451  return;
452 }

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshEdges::EdgeSplits(), global_gnomonic_projection(), gnomonic_projection_plane_at_point(), moab::Intx2MeshOnSphere::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().

◆ deep_copy_set_with_quads()

ErrorCode moab::IntxUtils::deep_copy_set_with_quads ( Interface *  mb,
EntityHandle  source_set,
EntityHandle  dest_set 
)
static

Definition at line 2255 of file IntxUtils.cpp.

2256 {
2257  ReadUtilIface* read_iface;
2258  MB_CHK_ERR( mb->query_interface( read_iface ) );
2259  // create the handle tag for the corresponding element / vertex
2260 
2261  EntityHandle dum = 0;
2262  Tag corrTag = 0; // it will be created here
2263  MB_CHK_ERR( mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum ) );
2264 
2265  // give the same global id to new verts and cells created in the lagr(departure) mesh
2266  Tag gid = mb->globalId_tag();
2267 
2268  Range quads;
2269  MB_CHK_ERR( mb->get_entities_by_type( source_set, MBQUAD, quads ) );
2270 
2271  Range connecVerts;
2272  MB_CHK_ERR( mb->get_connectivity( quads, connecVerts ) );
2273 
2274  std::map< EntityHandle, EntityHandle > newNodes;
2275 
2276  std::vector< double* > coords;
2277  EntityHandle start_vert, start_elem, *connect;
2278  int num_verts = connecVerts.size();
2279  MB_CHK_ERR( read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords ) );
2280 
2281  // fill it up
2282  int i = 0;
2283  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
2284  {
2285  EntityHandle oldV = *vit;
2286  CartVect posi;
2287  MB_CHK_ERR( mb->get_coords( &oldV, 1, &( posi[0] ) ) );
2288 
2289  int global_id;
2290  MB_CHK_ERR( mb->tag_get_data( gid, &oldV, 1, &global_id ) );
2291  EntityHandle new_vert = start_vert + i;
2292  // Cppcheck warning (false positive): variable coords is assigned a value that is never used
2293  coords[0][i] = posi[0];
2294  coords[1][i] = posi[1];
2295  coords[2][i] = posi[2];
2296 
2297  newNodes[oldV] = new_vert;
2298  // set also the correspondent tag :)
2299  MB_CHK_ERR( mb->tag_set_data( corrTag, &oldV, 1, &new_vert ) );
2300 
2301  // also the other side
2302  // need to check if we really need this; the new vertex will never need the old vertex
2303  // we have the global id which is the same
2304  MB_CHK_ERR( mb->tag_set_data( corrTag, &new_vert, 1, &oldV ) );
2305  // set the global id on the corresponding vertex the same as the initial vertex
2306  MB_CHK_ERR( mb->tag_set_data( gid, &new_vert, 1, &global_id ) );
2307  }
2308  // now create new quads in order (in a sequence)
2309 
2310  MB_CHK_ERR( read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect ) );
2311 
2312  int ie = 0;
2313  for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
2314  {
2315  EntityHandle q = *it;
2316  int nnodes;
2317  const EntityHandle* conn;
2318  MB_CHK_ERR( mb->get_connectivity( q, conn, nnodes ) );
2319  int global_id;
2320  MB_CHK_ERR( mb->tag_get_data( gid, &q, 1, &global_id ) );
2321 
2322  for( int ii = 0; ii < nnodes; ii++ )
2323  {
2324  EntityHandle v1 = conn[ii];
2325  connect[4 * ie + ii] = newNodes[v1];
2326  }
2327  EntityHandle newElement = start_elem + ie;
2328 
2329  // set the corresponding tag; not sure we need this one, from old to new
2330  MB_CHK_ERR( mb->tag_set_data( corrTag, &q, 1, &newElement ) );
2331  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &q ) );
2332 
2333  // set the global id
2334  MB_CHK_ERR( mb->tag_set_data( gid, &newElement, 1, &global_id ) );
2335 
2336  MB_CHK_ERR( mb->add_entities( dest_set, &newElement, 1 ) );
2337  }
2338 
2339  MB_CHK_ERR( read_iface->update_adjacencies( start_elem, quads.size(), 4, connect ) );
2340 
2341  return MB_SUCCESS;
2342 }

References moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), moab::ReadUtilIface::get_element_connect(), moab::ReadUtilIface::get_node_coords(), mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, MBQUAD, moab::Range::size(), and moab::ReadUtilIface::update_adjacencies().

◆ dist2()

static double moab::IntxUtils::dist2 ( double *  a,
double *  b 
)
inlinestatic

Definition at line 19 of file IntxUtils.hpp.

20  {
21  double abx = b[0] - a[0], aby = b[1] - a[1];
22  return sqrt( abx * abx + aby * aby );
23  }

Referenced by moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and SortAndRemoveDoubles2().

◆ distance_on_great_circle()

double moab::IntxUtils::distance_on_great_circle ( CartVect &  p1,
CartVect &  p2 
)
static

Definition at line 1546 of file IntxUtils.cpp.

1547 {
1548  SphereCoords sph1 = cart_to_spherical( p1 );
1549  SphereCoords sph2 = cart_to_spherical( p2 );
1550  // radius should be the same
1551  return sph1.R *
1552  acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
1553 }

References cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

◆ distance_on_sphere()

double moab::IntxUtils::distance_on_sphere ( double  la1,
double  te1,
double  la2,
double  te2 
)
static

Definition at line 1791 of file IntxUtils.cpp.

1792 {
1793  return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
1794 }

◆ EdgeIntersections2()

ErrorCode moab::IntxUtils::EdgeIntersections2 ( double *  blue,
int  nsBlue,
double *  red,
int  nsRed,
int *  markb,
int *  markr,
double *  points,
int &  nPoints 
)
static

Computes the edge intersections of two elements.

Parameters
blueThe array of points representing the blue element.
nsBlueThe number of points in the blue element.
redThe array of points representing the red element.
nsRedThe number of points in the red element.
markbThe array to mark the intersecting edges of the blue element.
markrThe array to mark the intersecting edges of the red element.
pointsThe array to store the intersection points.
nPointsThe number of intersection points found.
Returns
The error code.

Definition at line 221 of file IntxUtils.cpp.

229 {
230  /* EDGEINTERSECTIONS computes edge intersections of two elements
231  [P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
232  and blue ( stored column wise )
233  (point coordinates are stored column-wise, in counter clock
234  order) the points P where their edges intersect. In addition,
235  in n the indices of which neighbors of red are also intersecting
236  with blue are given.
237  */
238 
239  // points is an array with enough slots (24 * 2 doubles)
240  nPoints = 0;
241  for( int i = 0; i < MAXEDGES; i++ )
242  {
243  markb[i] = markr[i] = 0;
244  }
245 
246  for( int i = 0; i < nsBlue; i++ )
247  {
248  for( int j = 0; j < nsRed; j++ )
249  {
250  double b[2];
251  double a[2][2]; // 2*2
252  int iPlus1 = ( i + 1 ) % nsBlue;
253  int jPlus1 = ( j + 1 ) % nsRed;
254  for( int k = 0; k < 2; k++ )
255  {
256  b[k] = red[2 * j + k] - blue[2 * i + k];
257  // row k of a: a(k, 0), a(k, 1)
258  a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
259  a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
260  }
261  double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
262  if( fabs( delta ) > 1.e-14 ) // this is close to machine epsilon
263  {
264  // not parallel
265  double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
266  double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
267  if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
268  {
269  // the intersection is good
270  for( int k = 0; k < 2; k++ )
271  {
272  points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
273  }
274  markb[i] = 1; // so neighbor number i of blue will be considered too.
275  markr[j] = 1; // this will be used in advancing red around blue quad
276  nPoints++;
277  }
278  }
279  // the case delta ~ 0. will be considered by the interior points logic
280  }
281  }
282  return MB_SUCCESS;
283 }

References MAXEDGES, and MB_SUCCESS.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().

◆ EdgeIntxRllCs()

ErrorCode moab::IntxUtils::EdgeIntxRllCs ( double *  blue,
CartVect *  bluec,
int *  blueEdgeType,
int  nsBlue,
double *  red,
CartVect *  redc,
int  nsRed,
int *  markb,
int *  markr,
int  plane,
double  R,
double *  points,
int &  nPoints 
)
static

Computes the edge intersections between a RLL and CS quad.

Note: Special function.

Parameters
blueThe array of points representing the blue element.
bluecThe array of Cartesian coordinates for the blue element.
blueEdgeTypeThe array of edge types for the blue element.
nsBlueThe number of points in the blue element.
redThe array of points representing the red element.
redcThe array of Cartesian coordinates for the red element.
nsRedThe number of points in the red element.
markbThe array to mark the intersecting edges of the blue element.
markrThe array to mark the intersecting edges of the red element.
planeThe plane of intersection.
RThe radius of the sphere.
pointsThe array to store the intersection points.
nPointsThe number of intersection points found.
Returns
The error code.

Definition at line 305 of file IntxUtils.cpp.

318 {
319  // if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
320  // everything else the same (except if there are 2 points resulting, which is rare)
321  for( int i = 0; i < 4; i++ )
322  { // always at most 4 , so maybe don't bother
323  markb[i] = markr[i] = 0;
324  }
325 
326  for( int i = 0; i < nsBlue; i++ )
327  {
328  int iPlus1 = ( i + 1 ) % nsBlue;
329  if( blueEdgeType[i] == 0 ) // old style, just 2d
330  {
331  for( int j = 0; j < nsRed; j++ )
332  {
333  double b[2];
334  double a[2][2]; // 2*2
335 
336  int jPlus1 = ( j + 1 ) % nsRed;
337  for( int k = 0; k < 2; k++ )
338  {
339  b[k] = red[2 * j + k] - blue[2 * i + k];
340  // row k of a: a(k, 0), a(k, 1)
341  a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
342  a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
343  }
344  double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
345  if( fabs( delta ) > 1.e-14 )
346  {
347  // not parallel
348  double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
349  double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
350  if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
351  {
352  // the intersection is good
353  for( int k = 0; k < 2; k++ )
354  {
355  points[2 * nPoints + k] =
356  blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
357  }
358  markb[i] = 1; // so neighbor number i of blue will be considered too.
359  markr[j] = 1; // this will be used in advancing red around blue quad
360  nPoints++;
361  }
362  } // if the edges are too "parallel", skip them
363  }
364  }
365  else // edge type is 1, so use 3d intersection
366  {
367  CartVect& C = bluec[i];
368  CartVect& D = bluec[iPlus1];
369  for( int j = 0; j < nsRed; j++ )
370  {
371  int jPlus1 = ( j + 1 ) % nsRed; // nsRed is just 4, forget about it, usually
372  CartVect& A = redc[j];
373  CartVect& B = redc[jPlus1];
374  int np = 0;
375  double E[9];
376  intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
377  if( np == 0 ) continue;
378  if( np >= 2 )
379  {
380  std::cout << "intersection with 2 points :" << A << B << C << D << "\n";
381  }
382  for( int k = 0; k < np; k++ )
383  {
384  gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
385  points[2 * nPoints + 1] );
386  nPoints++;
387  }
388  markb[i] = 1; // so neighbor number i of blue will be considered too.
389  markr[j] = 1; // this will be used in advancing red around blue quad
390  }
391  }
392  }
393  return MB_SUCCESS;
394 }

References moab::CartVect::array(), moab::E, gnomonic_projection(), intersect_great_circle_arc_with_clat_arc(), MB_SUCCESS, and moab::R.

Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

◆ enforce_convexity()

ErrorCode moab::IntxUtils::enforce_convexity ( Interface *  mb,
EntityHandle  lset,
int  my_rank = 0 
)
static

Enforces convexity for a given set of polygons.

This function checks each polygon in the input set and computes the angles of each vertex. If a reflex angle is found, the polygon is broken into triangles and added back to the set. This process continues until all polygons in the set are convex.

Parameters
mbThe interface to the MOAB instance.
lsetThe handle of the input set containing the polygons.
my_rankThe rank of the local process.
Returns
The error code indicating the success or failure of the operation.

Definition at line 1568 of file IntxUtils.cpp.

1569 {
1570  Range inputRange;
1571  MB_CHK_ERR( mb->get_entities_by_dimension( lset, 2, inputRange ) );
1572 
1573  Tag corrTag = nullptr;
1574  EntityHandle dumH = 0;
1575  // no need to check return error
1576  mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
1577 
1578  Tag gidTag = mb->globalId_tag();
1579 
1580  std::vector< double > coords;
1581  coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
1582  // we should create a queue with new polygons that need processing for reflex angles
1583  // (obtuse)
1584  std::queue< EntityHandle > newPolys;
1585  int brokenPolys = 0;
1586  Range::iterator eit = inputRange.begin();
1587  while( eit != inputRange.end() || !newPolys.empty() )
1588  {
1589  EntityHandle eh;
1590  if( eit != inputRange.end() )
1591  {
1592  eh = *eit;
1593  ++eit;
1594  }
1595  else
1596  {
1597  eh = newPolys.front();
1598  newPolys.pop();
1599  }
1600  // get the nodes, then the coordinates
1601  const EntityHandle* verts;
1602  int num_nodes;
1603  MB_CHK_ERR( mb->get_connectivity( eh, verts, num_nodes ) );
1604  int nsides = num_nodes;
1605  // account for possible padded polygons
1606  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1607  nsides--;
1608  EntityHandle corrHandle = 0;
1609  if( corrTag )
1610  {
1611  MB_CHK_ERR( mb->tag_get_data( corrTag, &eh, 1, &corrHandle ) );
1612  }
1613  int gid = 0;
1614  MB_CHK_ERR( mb->tag_get_data( gidTag, &eh, 1, &gid ) );
1615  coords.resize( 3 * nsides );
1616  if( nsides < 4 ) continue; // if already triangles, don't bother
1617  // get coordinates
1618  MB_CHK_ERR( mb->get_coords( verts, nsides, &coords[0] ) );
1619  // compute each angle
1620  bool alreadyBroken = false;
1621 
1622  for( int i = 0; i < nsides; i++ )
1623  {
1624  double* A = &coords[3 * i];
1625  double* B = &coords[3 * ( ( i + 1 ) % nsides )];
1626  double* C = &coords[3 * ( ( i + 2 ) % nsides )];
1627  double angle = IntxUtils::oriented_spherical_angle( A, B, C );
1628  if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
1629  {
1630  if( alreadyBroken )
1631  {
1632  mb->list_entities( &eh, 1 );
1633  mb->list_entities( verts, nsides );
1634  double* D = &coords[3 * ( ( i + 3 ) % nsides )];
1635  std::cout << "ABC: " << angle << " \n";
1636  std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
1637  std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
1638  std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
1639  std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
1640 
1641  return MB_FAILURE;
1642  }
1643  // the bad angle is at i+1;
1644  // create 1 triangle and one polygon; add the polygon to the input range, so
1645  // it will be processed too
1646  // also, add both to the set :) and remove the original polygon from the set
1647  // break the next triangle, even though not optimal
1648  // so create the triangle i+1, i+2, i+3; remove i+2 from original list
1649  // even though not optimal in general, it is good enough.
1650  EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
1651  verts[( i + 3 ) % nsides] };
1652  // create a polygon with num_nodes-1 vertices, and connectivity
1653  // verts[i+1], verts[i+3], (all except i+2)
1654  std::vector< EntityHandle > conn( nsides - 1 );
1655  for( int j = 1; j < nsides; j++ )
1656  {
1657  conn[j - 1] = verts[( i + j + 2 ) % nsides];
1658  }
1659  EntityHandle newElement;
1660  MB_CHK_ERR( mb->create_element( MBTRI, conn3, 3, newElement ) );
1661 
1662  MB_CHK_ERR( mb->add_entities( lset, &newElement, 1 ) );
1663  if( corrTag )
1664  {
1665  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &corrHandle ) );
1666  }
1667  MB_CHK_ERR( mb->tag_set_data( gidTag, &newElement, 1, &gid ) );
1668  if( nsides == 4 )
1669  {
1670  // create another triangle
1671  MB_CHK_ERR( mb->create_element( MBTRI, &conn[0], 3, newElement ) );
1672  }
1673  else
1674  {
1675  // create another polygon, and add it to the inputRange
1676  MB_CHK_ERR( mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement ) );
1677  newPolys.push( newElement ); // because it has less number of edges, the
1678  // reverse should work to find it.
1679  }
1680  MB_CHK_ERR( mb->add_entities( lset, &newElement, 1 ) );
1681  if( corrTag )
1682  {
1683  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &corrHandle ) );
1684  }
1685  MB_CHK_ERR( mb->tag_set_data( gidTag, &newElement, 1, &gid ) );
1686  MB_CHK_ERR( mb->remove_entities( lset, &eh, 1 ) );
1687  brokenPolys++;
1688  alreadyBroken = true; // get out of the loop, element is broken
1689  }
1690  }
1691  }
1692  if( brokenPolys > 0 )
1693  {
1694  std::cout << "on local process " << my_rank << ", " << brokenPolys
1695  << " concave polygons were decomposed in convex ones \n";
1696 #ifdef VERBOSE
1697  std::stringstream fff;
1698  fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
1699  MB_CHK_ERR( mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 ) );
1700  std::cout << "wrote new file set: " << fff.str() << "\n";
1701 #endif
1702  }
1703  return MB_SUCCESS;
1704 }

References moab::angle(), moab::Range::begin(), CORRTAGNAME, moab::Range::end(), MAXEDGES, mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TYPE_HANDLE, MBPOLYGON, MBTRI, and oriented_spherical_angle().

Referenced by main().

◆ fix_degenerate_quads()

ErrorCode moab::IntxUtils::fix_degenerate_quads ( Interface *  mb,
EntityHandle  set 
)
static

Definition at line 1708 of file IntxUtils.cpp.

1709 {
1710  Range quads;
1711  MB_CHK_ERR( mb->get_entities_by_type( set, MBQUAD, quads ) );
1712 
1713  Tag gid = mb->globalId_tag();
1714  for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
1715  {
1716  EntityHandle quad = *qit;
1717  const EntityHandle* conn4 = nullptr;
1718  int num_nodes = 0;
1719  MB_CHK_ERR( mb->get_connectivity( quad, conn4, num_nodes ) );
1720  for( int i = 0; i < num_nodes; i++ )
1721  {
1722  int next_node_index = ( i + 1 ) % num_nodes;
1723  if( conn4[i] == conn4[next_node_index] )
1724  {
1725  // form a triangle and delete the quad
1726  // first get the global id, to set it on triangle later
1727  int global_id = 0;
1728  MB_CHK_ERR( mb->tag_get_data( gid, &quad, 1, &global_id ) );
1729  int i2 = ( i + 2 ) % num_nodes;
1730  int i3 = ( i + 3 ) % num_nodes;
1731  EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
1732  EntityHandle tri;
1733  MB_CHK_ERR( mb->create_element( MBTRI, conn3, 3, tri ) );
1734  MB_CHK_ERR( mb->add_entities( set, &tri, 1 ) );
1735  MB_CHK_ERR( mb->remove_entities( set, &quad, 1 ) );
1736  MB_CHK_ERR( mb->delete_entities( &quad, 1 ) );
1737  MB_CHK_ERR( mb->tag_set_data( gid, &tri, 1, &global_id ) );
1738  }
1739  }
1740  }
1741  return MB_SUCCESS;
1742 }

References moab::Range::begin(), moab::Range::end(), mb, MB_CHK_ERR, MB_SUCCESS, MBQUAD, and MBTRI.

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().

◆ global_gnomonic_projection()

ErrorCode moab::IntxUtils::global_gnomonic_projection ( Interface *  mb,
EntityHandle  inSet,
double  R,
bool  centers_only,
EntityHandle &  outSet 
)
static

Definition at line 840 of file IntxUtils.cpp.

845 {
846  std::string parTagName( "PARALLEL_PARTITION" );
847  Tag part_tag;
848  Tag gidTag = mb->globalId_tag();
849  Tag targetParentTag, sourceParentTag;
850  mb->tag_get_handle( "TargetParent", targetParentTag );
851  mb->tag_get_handle( "SourceParent", sourceParentTag );
852  bool intxMesh = false;
853  if( targetParentTag != nullptr && sourceParentTag != nullptr )
854  intxMesh = true; // interested in source and target parent tags then
855  Range partSets;
856  ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
857  if( MB_SUCCESS == rval && part_tag != 0 )
858  {
859  MB_CHK_ERR(
860  mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, nullptr, 1, partSets, Interface::UNION ) );
861  }
862  MB_CHK_ERR( ScaleToRadius( mb, inSet, 1.0 ) );
863  // Get all entities of dimension 2
864  Range inputRange; // get
865  MB_CHK_ERR( mb->get_entities_by_dimension( inSet, 1, inputRange ) );
866  MB_CHK_ERR( mb->get_entities_by_dimension( inSet, 2, inputRange ) );
867 
868  std::map< EntityHandle, int > partsAssign;
869  std::map< int, EntityHandle > newPartSets;
870  if( !partSets.empty() )
871  {
872  // get all cells, and assign parts
873  for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
874  {
875  EntityHandle pSet = *setIt;
876  Range ents;
877  MB_CHK_ERR( mb->get_entities_by_handle( pSet, ents ) );
878  int val;
879  MB_CHK_ERR( mb->tag_get_data( part_tag, &pSet, 1, &val ) );
880  // create a new set with the same part id tag, in the outSet
881  EntityHandle newPartSet;
882  MB_CHK_ERR( mb->create_meshset( MESHSET_SET, newPartSet ) );
883  MB_CHK_ERR( mb->tag_set_data( part_tag, &newPartSet, 1, &val ) );
884  newPartSets[val] = newPartSet;
885  MB_CHK_ERR( mb->add_entities( outSet, &newPartSet, 1 ) );
886  for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
887  {
888  partsAssign[*it] = val;
889  }
890  }
891  }
892 
893  if( centers_only )
894  {
895  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
896  {
897  CartVect center;
898  EntityHandle cell = *it;
899  MB_CHK_ERR( mb->get_coords( &cell, 1, center.array() ) );
900  int globalID = 0;
901  if( !intxMesh )
902  {
903  MB_CHK_SET_ERR( mb->tag_get_data( gidTag, &cell, 1, &globalID ), "can't get id tag on cell" );
904  }
905 
906  int plane;
907  decide_gnomonic_plane( center, plane );
908  double c[3];
909  c[2] = 0.;
910  gnomonic_projection( center, R, plane, c[0], c[1] );
911 
912  gnomonic_unroll( c[0], c[1], R, plane );
913 
914  EntityHandle vertex;
915  MB_CHK_ERR( mb->create_vertex( c, vertex ) );
916 
917  if( !intxMesh )
918  {
919  MB_CHK_SET_ERR( mb->tag_set_data( gidTag, &vertex, 1, &globalID ), "can't set id tag on center" );
920  }
921  MB_CHK_ERR( mb->add_entities( outSet, &vertex, 1 ) );
922  }
923  }
924  else
925  {
926  // distribute the cells to 6 planes, based on the center
927  Range subranges[6];
928  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
929  {
930  CartVect center;
931  EntityHandle cell = *it;
932  MB_CHK_ERR( mb->get_coords( &cell, 1, center.array() ) );
933  int plane;
934  decide_gnomonic_plane( center, plane );
935  subranges[plane - 1].insert( cell ); // includes edges if they exist
936  }
937  for( int i = 1; i <= 6; i++ )
938  {
939  Range verts;
940  MB_CHK_ERR( mb->get_connectivity( subranges[i - 1], verts ) );
941  std::map< EntityHandle, EntityHandle > corr;
942  for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
943  {
944  CartVect vect;
945  EntityHandle v = *vt;
946  MB_CHK_ERR( mb->get_coords( &v, 1, vect.array() ) );
947  double c[3];
948  c[2] = 0.;
949  gnomonic_projection( vect, R, i, c[0], c[1] );
950  gnomonic_unroll( c[0], c[1], R, i );
951  EntityHandle vertex;
952  MB_CHK_ERR( mb->create_vertex( c, vertex ) );
953 
954  int vID;
955  if( !intxMesh )
956  {
957  MB_CHK_SET_ERR( mb->tag_get_data( gidTag, &v, 1, &vID ), "can't get id tag on vertex" );
958  // new vertex will get old ID
959  MB_CHK_SET_ERR( mb->tag_set_data( gidTag, &vertex, 1, &vID ), "can't get id tag on vertex" );
960  }
961  corr[v] = vertex; // for new connectivity
962  }
963  EntityHandle new_conn[20]; // max edges in 2d ?
964  for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
965  {
966  EntityHandle eh = *eit;
967  const EntityHandle* conn = nullptr;
968  int num_nodes;
969  MB_CHK_ERR( mb->get_connectivity( eh, conn, num_nodes ) );
970 
971  // build a new vertex array
972  for( int j = 0; j < num_nodes; j++ )
973  new_conn[j] = corr[conn[j]];
974 
975  EntityType type = mb->type_from_handle( eh );
976  EntityHandle newCell;
977  MB_CHK_ERR( mb->create_element( type, new_conn, num_nodes, newCell ) );
978  MB_CHK_ERR( mb->add_entities( outSet, &newCell, 1 ) );
979 
980  int eID;
981  if( !intxMesh )
982  {
983  MB_CHK_SET_ERR( mb->tag_get_data( gidTag, &eh, 1, &eID ), "can't get id tag on entity handle" );
984  // new vertex will get old ID
985  MB_CHK_SET_ERR( mb->tag_set_data( gidTag, &newCell, 1, &eID ), "can't set id tag on new cell" );
986  }
987  else
988  {
989  // look for parent tags if intx mesh targetParentTag , sourceParentTag
990  if( type >= moab::MBPOLYGON )
991  {
992  MB_CHK_SET_ERR( mb->tag_get_data( targetParentTag, &eh, 1, &eID ),
993  "can't get parent tag on entity handle" );
994  MB_CHK_SET_ERR( mb->tag_set_data( targetParentTag, &newCell, 1, &eID ),
995  "can't set parent tag on entity handle" );
996  MB_CHK_SET_ERR( mb->tag_get_data( sourceParentTag, &eh, 1, &eID ),
997  "can't get parent tag on entity handle" );
998  MB_CHK_SET_ERR( mb->tag_set_data( sourceParentTag, &newCell, 1, &eID ),
999  "can't set parent tag on entity handle" );
1000  }
1001  }
1002 
1003  std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
1004  if( mit != partsAssign.end() )
1005  {
1006  int val = mit->second;
1007  MB_CHK_ERR( mb->add_entities( newPartSets[val], &newCell, 1 ) );
1008  }
1009  }
1010  }
1011  }
1012 
1013  return MB_SUCCESS;
1014 }

References moab::CartVect::array(), moab::Range::begin(), center(), decide_gnomonic_plane(), moab::Range::empty(), moab::Range::end(), ErrorCode, gnomonic_projection(), gnomonic_unroll(), moab::Range::insert(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBENTITYSET, MBPOLYGON, MESHSET_SET, moab::R, ScaleToRadius(), and moab::Interface::UNION.

◆ global_gnomonic_projection_general()

ErrorCode moab::IntxUtils::global_gnomonic_projection_general ( Interface *  mb,
EntityHandle  inSet,
CartVect  P,
EntityHandle &  outSet 
)
static

Definition at line 716 of file IntxUtils.cpp.

720 {
721  std::string parTagName( "PARALLEL_PARTITION" );
722  Tag part_tag;
723  Tag gidTag = mb->globalId_tag();
724  Tag targetParentTag, sourceParentTag;
725  mb->tag_get_handle( "TargetParent", targetParentTag );
726  mb->tag_get_handle( "SourceParent", sourceParentTag );
727  bool intxMesh = false;
728  if( targetParentTag != nullptr && sourceParentTag != nullptr )
729  intxMesh = true; // interested in source and target parent tags then
730  Range partSets;
731  ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
732  if( MB_SUCCESS == rval && part_tag != 0 )
733  {
734  rval =
735  mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, nullptr, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
736  }
737  rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
738  // Get all entities of dimension 2
739  Range inputRange; // get
740  rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
741  rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
742 
743  std::map< EntityHandle, int > partsAssign;
744  std::map< int, EntityHandle > newPartSets;
745  if( !partSets.empty() )
746  {
747  // get all cells, and assign parts
748  for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
749  {
750  EntityHandle pSet = *setIt;
751  Range ents;
752  rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
753  int val;
754  rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
755  // create a new set with the same part id tag, in the outSet
756  EntityHandle newPartSet;
757  rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
758  rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
759  newPartSets[val] = newPartSet;
760  rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
761  for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
762  {
763  partsAssign[*it] = val;
764  }
765  }
766  }
767 
768  // decide gnomonic plane
769  CartVect axis[3];
770  axis[0] = P;
771  IntxUtils::gnomonic_projection_plane_at_point( axis[0], axis[1], axis[2] );
772  // project all vertices, and then create new cells
773 
774  Range verts;
775  rval = mb->get_connectivity( inputRange, verts );MB_CHK_ERR( rval );
776  std::map< EntityHandle, EntityHandle > corr;
777  for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
778  {
779  CartVect vect;
780  EntityHandle v = *vt;
781  rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
782  double c[3];
783  c[2] = 0.;
784  IntxUtils::gnomonic_projection_generalized( vect, axis, c[0], c[1] );
785 
786  EntityHandle vertex;
787  rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
788  int vID;
789  if( !intxMesh )
790  {
791  rval = mb->tag_get_data( gidTag, &v, 1, &vID );MB_CHK_SET_ERR( rval, "can't get id tag on vertex" );
792  // new vertex will get old ID
793  rval = mb->tag_set_data( gidTag, &vertex, 1, &vID );MB_CHK_SET_ERR( rval, "can't get id tag on vertex" );
794  }
795  corr[v] = vertex; // for new connectivity
796  }
797  EntityHandle new_conn[20]; // max edges in 2d ?
798  for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
799  {
800  EntityHandle eh = *eit;
801  const EntityHandle* conn = nullptr;
802  int num_nodes;
803  rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
804  // build a new vertex array
805  for( int j = 0; j < num_nodes; j++ )
806  new_conn[j] = corr[conn[j]];
807  EntityType type = mb->type_from_handle( eh );
808  EntityHandle newCell;
809  rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
810  rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
811  int eID;
812  if( !intxMesh )
813  {
814  rval = mb->tag_get_data( gidTag, &eh, 1, &eID );MB_CHK_SET_ERR( rval, "can't get id tag on entity handle" );
815  // new vertex will get old ID
816  rval = mb->tag_set_data( gidTag, &newCell, 1, &eID );MB_CHK_SET_ERR( rval, "can't set id tag on new cell" );
817  }
818  else
819  {
820  // look for parent tags if intx mesh targetParentTag , sourceParentTag
821  if( type >= moab::MBPOLYGON )
822  {
823  rval = mb->tag_get_data( targetParentTag, &eh, 1, &eID );MB_CHK_SET_ERR( rval, "can't get parent tag on entity handle" );
824  rval = mb->tag_set_data( targetParentTag, &newCell, 1, &eID );MB_CHK_SET_ERR( rval, "can't set parent tag on entity handle" );
825  rval = mb->tag_get_data( sourceParentTag, &eh, 1, &eID );MB_CHK_SET_ERR( rval, "can't get parent tag on entity handle" );
826  rval = mb->tag_set_data( sourceParentTag, &newCell, 1, &eID );MB_CHK_SET_ERR( rval, "can't set parent tag on entity handle" );
827  }
828  }
829  std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
830  if( mit != partsAssign.end() )
831  {
832  int val = mit->second;
833  rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
834  }
835  }
836  return MB_SUCCESS;
837 }

References moab::CartVect::array(), moab::Range::begin(), moab::Range::empty(), moab::Range::end(), ErrorCode, gnomonic_projection_generalized(), gnomonic_projection_plane_at_point(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBENTITYSET, MBPOLYGON, MESHSET_SET, ScaleToRadius(), and moab::Interface::UNION.

Referenced by main().

◆ gnomonic_projection()

ErrorCode moab::IntxUtils::gnomonic_projection ( const CartVect &  pos,
double  R,
int  plane,
double &  c1,
double &  c2 
)
static

Definition at line 559 of file IntxUtils.cpp.

560 {
561  double alfa = 1.; // the new point will be on line alfa*pos
562 
563  switch( plane )
564  {
565  case 1: {
566  // the plane with x = R; x>y, x>z
567  // c1->y, c2->z
568  alfa = R / pos[0];
569  c1 = alfa * pos[1];
570  c2 = alfa * pos[2];
571  break;
572  }
573  case 2: {
574  // y = R -> zx
575  alfa = R / pos[1];
576  c1 = alfa * pos[2];
577  c2 = alfa * pos[0];
578  break;
579  }
580  case 3: {
581  // x=-R, -> yz
582  alfa = -R / pos[0];
583  c1 = -alfa * pos[1]; // the sign is to preserve orientation
584  c2 = alfa * pos[2];
585  break;
586  }
587  case 4: {
588  // y = -R
589  alfa = -R / pos[1];
590  c1 = -alfa * pos[2]; // the sign is to preserve orientation
591  c2 = alfa * pos[0];
592  break;
593  }
594  case 5: {
595  // z = -R
596  alfa = -R / pos[2];
597  c1 = -alfa * pos[0]; // the sign is to preserve orientation
598  c2 = alfa * pos[1];
599  break;
600  }
601  case 6: {
602  alfa = R / pos[2];
603  c1 = alfa * pos[0];
604  c2 = alfa * pos[1];
605  break;
606  }
607  default:
608  return MB_FAILURE; // error
609  }
610 
611  return MB_SUCCESS; // no error
612 }

References MB_SUCCESS, and moab::R.

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), EdgeIntxRllCs(), moab::Intx2MeshEdges::EdgeSplits(), global_gnomonic_projection(), moab::Intx2MeshOnSphere::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().

◆ gnomonic_projection_generalized()

ErrorCode moab::IntxUtils::gnomonic_projection_generalized ( const CartVect &  pos,
const CartVect  axis[3],
double &  c1,
double &  c2 
)
static

Definition at line 542 of file IntxUtils.cpp.

546 {
547  double ang = angle( pos, axis[0] );
548  if( ang > 1.57 ) // pi/2 do not project if very close to hemisphere
549  return MB_FAILURE;
550  // solve the equation in plane, (alfa * pos - axis[0]) % axis[0] = 0.0
551  double alpha = axis[0] % axis[0] / ( pos % axis[0] ); // we know this denominator is greater than 0
552  CartVect planeVect = alpha * pos - axis[0]; // axis[0] is P
553  c1 = planeVect % axis[1];
554  c2 = planeVect % axis[2];
555  return MB_SUCCESS;
556 }

References moab::angle(), and MB_SUCCESS.

Referenced by global_gnomonic_projection_general().

◆ gnomonic_projection_plane_at_point()

ErrorCode moab::IntxUtils::gnomonic_projection_plane_at_point ( CartVect  P,
CartVect &  u,
CartVect &  v 
)
static

Definition at line 454 of file IntxUtils.cpp.

455 {
456 
457  double d = P.length();
458  if( d == 0.0 )
459  {
460  MB_CHK_SET_ERR( MB_FAILURE, "point P is at the origin" );
461  }
462  double x = P[0];
463  double y = P[1];
464  double z = P[2];
465  // easy cases
466  if( x == 0.0 && y == 0.0 )
467  {
468  if( z > 0. )
469  {
470  u = CartVect( 1., 0., 0. );
471  v = CartVect( 0., 1., 0. ); // gnomonic plane 6
472  }
473  else
474  {
475  u = CartVect( 0., 1., 0. );
476  v = CartVect( 1., 0., 0. ); // gnomonic plane 5
477  }
478  return MB_SUCCESS;
479  }
480  if( x == 0.0 && z == 0.0 )
481  {
482  if( y > 0. )
483  {
484  u = CartVect( -1., 0., 0. );
485  v = CartVect( 0., 0., 1. ); // gnomonic plane 2
486  }
487  else
488  {
489  u = CartVect( 0., 0., 1. );
490  v = CartVect( -1., 0., 0. ); // gnomonic plane 4
491  }
492  return MB_SUCCESS;
493  }
494  if( z == 0.0 && y == 0.0 )
495  {
496  if( x > 0. )
497  {
498  u = CartVect( 0., 1., 0. );
499  v = CartVect( 0., 0., 1. ); // gnomonic plane 1
500  }
501  else
502  {
503  u = CartVect( 0., 0., 1. );
504  v = CartVect( 0., 1., 0. ); // gnomonic plane 3
505  }
506  return MB_SUCCESS;
507  }
508  int plane;
509  IntxUtils::decide_gnomonic_plane( P, plane );
510  if( 1 == plane ) // towards x > 0
511  {
512  u = CartVect( 1., 0., 0. ) * P;
513  }
514 
515  if( 2 == plane ) // towards y > 0
516  {
517  u = CartVect( 0., 1., 0. ) * P;
518  }
519  if( 3 == plane ) // towards x < 0
520  {
521  u = CartVect( -1., 0., 0. ) * P;
522  }
523  if( 4 == plane ) // towards y < 0
524  {
525  u = CartVect( 0., -1., 0. ) * P;
526  }
527  if( 5 == plane ) // towards z < 0
528  {
529  u = CartVect( 0., 0., -1. ) * P;
530  }
531  if( 6 == plane ) // towards z > 0
532  {
533  u = CartVect( 0., 0., 1. ) * P;
534  }
535  v = P * u;
536  u.normalize();
537  v.normalize();
538 
539  return MB_SUCCESS;
540 }

References decide_gnomonic_plane(), moab::CartVect::length(), MB_CHK_SET_ERR, MB_SUCCESS, and moab::CartVect::normalize().

Referenced by global_gnomonic_projection_general().

◆ gnomonic_unroll()

void moab::IntxUtils::gnomonic_unroll ( double &  c1,
double &  c2,
double  R,
int  plane 
)
static

Definition at line 677 of file IntxUtils.cpp.

678 {
679  double tmp;
680  switch( plane )
681  {
682  case 1:
683  break; // nothing
684  case 2: // rotate + 90
685  tmp = c1;
686  c1 = -c2;
687  c2 = tmp;
688  c1 = c1 + 2 * R;
689  break;
690  case 3:
691  c1 = c1 + 4 * R;
692  break;
693  case 4: // rotate with -90 x-> -y; y -> x
694 
695  tmp = c1;
696  c1 = c2;
697  c2 = -tmp;
698  c1 = c1 - 2 * R;
699  break;
700  case 5: // South Pole
701  // rotate 180 then move to (-2, -2)
702  c1 = -c1 - 2. * R;
703  c2 = -c2 - 2. * R;
704  break;
705  ;
706  case 6: // North Pole
707  c1 = c1 - 2. * R;
708  c2 = c2 + 2. * R;
709  break;
710  }
711  return;
712 }

References moab::R.

Referenced by global_gnomonic_projection(), and transform_coordinates().

◆ intersect_great_circle_arc_with_clat_arc()

ErrorCode moab::IntxUtils::intersect_great_circle_arc_with_clat_arc ( double *  A,
double *  B,
double *  C,
double *  D,
double  R,
double *  E,
int &  np 
)
static

Definition at line 1887 of file IntxUtils.cpp.

1894 {
1895  const double distTol = R * 1.e-6;
1896  const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
1897  np = 0; // number of points in intersection
1898  CartVect a( A ), b( B ), c( C ), d( D );
1899  // check input first
1900  double R2 = R * R;
1901  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1902  fabs( d.length_squared() - R2 ) >
1903  10 * Tolerance )
1904  return MB_FAILURE;
1905 
1906  if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
1907  if( ( c - d ).length_squared() < Tolerance ) // edges are too short
1908  return MB_FAILURE;
1909 
1910  // CD is the const latitude arc
1911  if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
1912  return MB_FAILURE;
1913 
1914  if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
1915 
1916  // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
1917  // great circle arc AB normal to the AB circle:
1918  CartVect n1 = a * b; // the normal to the great circle arc (circle)
1919  // solve the system of equations:
1920  /*
1921  * n1%(x, y, z) = 0 // on the great circle
1922  * z = C[2];
1923  * x^2+y^2+z^2 = R^2
1924  */
1925  double z = C[2];
1926  if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
1927  {
1928  // it is the Equator; check if the const lat edge is Equator too
1929  if( fabs( C[2] ) > distTol )
1930  {
1931  return MB_FAILURE; // no intx, too far from Eq
1932  }
1933  else
1934  {
1935  // all points are on the equator
1936  //
1937  CartVect cd = c * d;
1938  // by convention, c<d, positive is from c to d
1939  // is a or b between c , d?
1940  CartVect ca = c * a;
1941  CartVect ad = a * d;
1942  CartVect cb = c * b;
1943  CartVect bd = b * d;
1944  bool agtc = ( ca % cd >= -Tolerance ); // a>c?
1945  bool dgta = ( ad % cd >= -Tolerance ); // d>a?
1946  bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
1947  bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
1948  if( agtc )
1949  {
1950  if( dgta )
1951  {
1952  // a is for sure a point
1953  E[0] = a[0];
1954  E[1] = a[1];
1955  E[2] = a[2];
1956  np++;
1957  if( bgtc )
1958  {
1959  if( dgtb )
1960  {
1961  // b is also in between c and d
1962  E[3] = b[0];
1963  E[4] = b[1];
1964  E[5] = b[2];
1965  np++;
1966  }
1967  else
1968  {
1969  // then order is c a d b, intx is ad
1970  E[3] = d[0];
1971  E[4] = d[1];
1972  E[5] = d[2];
1973  np++;
1974  }
1975  }
1976  else
1977  {
1978  // b is less than c, so b c a d, intx is ac
1979  E[3] = c[0];
1980  E[4] = c[1];
1981  E[5] = c[2];
1982  np++; // what if E[0] is E[3]?
1983  }
1984  }
1985  else // c < d < a
1986  {
1987  if( dgtb ) // d is for sure in
1988  {
1989  E[0] = d[0];
1990  E[1] = d[1];
1991  E[2] = d[2];
1992  np++;
1993  if( bgtc ) // c<b<d<a
1994  {
1995  // another point is b
1996  E[3] = b[0];
1997  E[4] = b[1];
1998  E[5] = b[2];
1999  np++;
2000  }
2001  else // b<c<d<a
2002  {
2003  // another point is c
2004  E[3] = c[0];
2005  E[4] = c[1];
2006  E[5] = c[2];
2007  np++;
2008  }
2009  }
2010  else
2011  {
2012  // nothing, order is c, d < a, b
2013  }
2014  }
2015  }
2016  else // a < c < d
2017  {
2018  if( bgtc )
2019  {
2020  // c is for sure in
2021  E[0] = c[0];
2022  E[1] = c[1];
2023  E[2] = c[2];
2024  np++;
2025  if( dgtb )
2026  {
2027  // a < c < b < d; second point is b
2028  E[3] = b[0];
2029  E[4] = b[1];
2030  E[5] = b[2];
2031  np++;
2032  }
2033  else
2034  {
2035  // a < c < d < b; second point is d
2036  E[3] = d[0];
2037  E[4] = d[1];
2038  E[5] = d[2];
2039  np++;
2040  }
2041  }
2042  else // a, b < c < d
2043  {
2044  // nothing
2045  }
2046  }
2047  }
2048  // for the 2 points selected, see if it is only one?
2049  // no problem, maybe it will be collapsed later anyway
2050  if( np > 0 ) return MB_SUCCESS;
2051  return MB_FAILURE; // no intersection
2052  }
2053  {
2054  if( fabs( n1[0] ) <= fabs( n1[1] ) )
2055  {
2056  // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
2057  // (u+v*x)^2+x^2=R2-z^2
2058  // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
2059  // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
2060  // x1,2 =
2061  double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
2062  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
2063  double delta = b1 * b1 - 4 * a1 * c1;
2064  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
2065  if( delta > Tolerance ) // 2 solutions possible
2066  {
2067  double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
2068  double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
2069  double y1 = u + v * x1;
2070  double y2 = u + v * x2;
2071  if( verify( a, b, c, d, x1, y1, z ) )
2072  {
2073  E[0] = x1;
2074  E[1] = y1;
2075  E[2] = z;
2076  np++;
2077  }
2078  if( verify( a, b, c, d, x2, y2, z ) )
2079  {
2080  E[3 * np + 0] = x2;
2081  E[3 * np + 1] = y2;
2082  E[3 * np + 2] = z;
2083  np++;
2084  }
2085  }
2086  else
2087  {
2088  // one solution
2089  double x1 = -b1 / 2 / a1;
2090  double y1 = u + v * x1;
2091  if( verify( a, b, c, d, x1, y1, z ) )
2092  {
2093  E[0] = x1;
2094  E[1] = y1;
2095  E[2] = z;
2096  np++;
2097  }
2098  }
2099  }
2100  else
2101  {
2102  // resolve eq in y, reverse
2103  // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
2104  // (u+v*y)^2+y^2 -R2+z^2 =0
2105  // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
2106  //
2107  // x1,2 =
2108  double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
2109  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
2110  double delta = b1 * b1 - 4 * a1 * c1;
2111  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
2112  if( delta > Tolerance ) // 2 solutions possible
2113  {
2114  double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
2115  double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
2116  double x1 = u + v * y1;
2117  double x2 = u + v * y2;
2118  if( verify( a, b, c, d, x1, y1, z ) )
2119  {
2120  E[0] = x1;
2121  E[1] = y1;
2122  E[2] = z;
2123  np++;
2124  }
2125  if( verify( a, b, c, d, x2, y2, z ) )
2126  {
2127  E[3 * np + 0] = x2;
2128  E[3 * np + 1] = y2;
2129  E[3 * np + 2] = z;
2130  np++;
2131  }
2132  }
2133  else
2134  {
2135  // one solution
2136  double y1 = -b1 / 2 / a1;
2137  double x1 = u + v * y1;
2138  if( verify( a, b, c, d, x1, y1, z ) )
2139  {
2140  E[0] = x1;
2141  E[1] = y1;
2142  E[2] = z;
2143  np++;
2144  }
2145  }
2146  }
2147  }
2148 
2149  if( np <= 0 ) return MB_FAILURE;
2150  return MB_SUCCESS;
2151 }

References moab::E, moab::CartVect::length_squared(), length_squared(), MB_SUCCESS, moab::R, and moab::verify().

Referenced by EdgeIntxRllCs().

◆ intersect_great_circle_arcs()

ErrorCode moab::IntxUtils::intersect_great_circle_arcs ( double *  A,
double *  B,
double *  C,
double *  D,
double  R,
double *  E 
)
static

Definition at line 1800 of file IntxUtils.cpp.

1801 {
1802  // first verify A, B, C, D are on the same sphere
1803  double R2 = R * R;
1804  const double Tolerance = 1.e-12 * R2;
1805 
1806  CartVect a( A ), b( B ), c( C ), d( D );
1807 
1808  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1809  fabs( d.length_squared() - R2 ) >
1810  10 * Tolerance )
1811  return MB_FAILURE;
1812 
1813  CartVect n1 = a * b;
1814  if( n1.length_squared() < Tolerance ) return MB_FAILURE;
1815 
1816  CartVect n2 = c * d;
1817  if( n2.length_squared() < Tolerance ) return MB_FAILURE;
1818  CartVect n3 = n1 * n2;
1819  n3.normalize();
1820 
1821  n3 = R * n3;
1822  // the intersection is either n3 or -n3
1823  CartVect n4 = a * n3, n5 = n3 * b;
1824  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1825  {
1826  // n3 is good for ab, see if it is good for cd
1827  n4 = c * n3;
1828  n5 = n3 * d;
1829  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1830  {
1831  E[0] = n3[0];
1832  E[1] = n3[1];
1833  E[2] = n3[2];
1834  }
1835  else
1836  return MB_FAILURE;
1837  }
1838  else
1839  {
1840  // try -n3
1841  n3 = -n3;
1842  n4 = a * n3, n5 = n3 * b;
1843  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1844  {
1845  // n3 is good for ab, see if it is good for cd
1846  n4 = c * n3;
1847  n5 = n3 * d;
1848  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1849  {
1850  E[0] = n3[0];
1851  E[1] = n3[1];
1852  E[2] = n3[2];
1853  }
1854  else
1855  return MB_FAILURE;
1856  }
1857  else
1858  return MB_FAILURE;
1859  }
1860 
1861  return MB_SUCCESS;
1862 }

References moab::E, moab::CartVect::length_squared(), MB_SUCCESS, moab::CartVect::normalize(), and moab::R.

◆ max_diagonal()

ErrorCode moab::IntxUtils::max_diagonal ( Interface *  mb,
Range  cells,
int  max_edges,
double &  diagonal 
)
static

Definition at line 2441 of file IntxUtils.cpp.

2442 {
2443  diagonal = 0.0;
2444  std::vector< CartVect > coords( max_edges ); // maximum number of nodes in a cell? hard coded ?
2445  for( auto it = cells.begin(); it != cells.end(); ++it )
2446  {
2447  // get the connectivity, then the coordinates
2448  EntityHandle cell = *it;
2449  const EntityHandle* connec = nullptr;
2450  int num_verts = 0;
2451  MB_CHK_SET_ERR( mb->get_connectivity( cell, connec, num_verts ), "Failed to get connectivity" );
2452  MB_CHK_SET_ERR( mb->get_coords( connec, num_verts, &( coords[0][0] ) ), "Failed to get coordinates" );
2453  // compute the max diagonal, in a double loop
2454  for( int i = 0; i < num_verts - 1; i++ )
2455  {
2456  for( int j = i + 1; j < num_verts; j++ )
2457  {
2458  double len_sq = ( coords[i] - coords[j] ).length_squared();
2459  if( len_sq > diagonal ) diagonal = len_sq;
2460  }
2461  }
2462  }
2463 
2464  // return the root of the diagonal
2465  diagonal = std::sqrt( diagonal );
2466  return MB_SUCCESS;
2467 }

References moab::Range::begin(), moab::Range::end(), length_squared(), mb, MB_CHK_SET_ERR, and MB_SUCCESS.

◆ oriented_spherical_angle()

double moab::IntxUtils::oriented_spherical_angle ( const double *  A,
const double *  B,
const double *  C 
)
static

Definition at line 1160 of file IntxUtils.cpp.

1161 {
1162  // assume the same radius, sphere at origin
1163  CartVect a( A ), b( B ), c( C );
1164  CartVect normalOAB = a * b;
1165  CartVect normalOCB = c * b;
1166  CartVect orient = ( c - b ) * ( a - b );
1167  double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
1168  if( ang != ang )
1169  {
1170  // signal of a nan
1171  std::cout << a << " " << b << " " << c << "\n";
1172  std::cout << ang << "\n";
1173  }
1174  if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
1175 
1176  return ang;
1177 }

References moab::angle().

Referenced by moab::IntxAreaUtils::area_spherical_polygon_girard(), and enforce_convexity().

◆ remove_duplicate_vertices()

ErrorCode moab::IntxUtils::remove_duplicate_vertices ( Interface *  mb,
EntityHandle  file_set,
double  merge_tol,
std::vector< Tag > &  tagList 
)
static

Definition at line 2344 of file IntxUtils.cpp.

2348 {
2349  Range verts;
2350  MB_CHK_ERR( mb->get_entities_by_dimension( file_set, 0, verts ) );
2351  MB_CHK_ERR( mb->remove_entities( file_set, verts ) );
2352 
2353  MergeMesh mm( mb );
2354 
2355  // remove the vertices from the set, before merging
2356 
2357  MB_CHK_ERR( mm.merge_all( file_set, merge_tol ) );
2358 
2359  // now correct vertices that are repeated in polygons
2360  MB_CHK_ERR( remove_padded_vertices( mb, file_set, tagList ) );
2361  return MB_SUCCESS;
2362 }

References mb, MB_CHK_ERR, MB_SUCCESS, moab::MergeMesh::merge_all(), and remove_padded_vertices().

Referenced by iMOAB_MergeVertices().

◆ remove_padded_vertices()

ErrorCode moab::IntxUtils::remove_padded_vertices ( Interface *  mb,
EntityHandle  file_set,
std::vector< Tag > &  tagList 
)
static

Definition at line 2364 of file IntxUtils.cpp.

2365 {
2366 
2367  // now correct vertices that are repeated in polygons
2368  Range cells;
2369  MB_CHK_ERR( mb->get_entities_by_dimension( file_set, 2, cells ) );
2370 
2371  Range verts;
2372  MB_CHK_ERR( mb->get_connectivity( cells, verts ) );
2373 
2374  Range modifiedCells; // will be deleted at the end; keep the gid
2375  Range newCells;
2376 
2377  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
2378  {
2379  EntityHandle cell = *cit;
2380  const EntityHandle* connec = nullptr;
2381  int num_verts = 0;
2382  MB_CHK_SET_ERR( mb->get_connectivity( cell, connec, num_verts ), "Failed to get connectivity" );
2383 
2384  std::vector< EntityHandle > newConnec;
2385  newConnec.push_back( connec[0] ); // at least one vertex
2386  int index = 0;
2387  int new_size = 1;
2388  while( index < num_verts - 2 )
2389  {
2390  int next_index = ( index + 1 );
2391  if( connec[next_index] != newConnec[new_size - 1] )
2392  {
2393  newConnec.push_back( connec[next_index] );
2394  new_size++;
2395  }
2396  index++;
2397  }
2398  // add the last one only if different from previous and first node
2399  if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
2400  {
2401  newConnec.push_back( connec[num_verts - 1] );
2402  new_size++;
2403  }
2404  if( new_size < num_verts )
2405  {
2406  // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
2407  modifiedCells.insert( cell );
2408  // create a new cell with type triangle, quad or polygon
2409  EntityType type = MBTRI;
2410  if( new_size == 3 )
2411  type = MBTRI;
2412  else if( new_size == 4 )
2413  type = MBQUAD;
2414  else if( new_size > 4 )
2415  type = MBPOLYGON;
2416 
2417  // create new cell
2418  EntityHandle newCell;
2419  MB_CHK_SET_ERR( mb->create_element( type, &newConnec[0], new_size, newCell ), "Failed to create new cell" );
2420  // set the old id to the new element
2421  newCells.insert( newCell );
2422  double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
2423  for( size_t i = 0; i < tagList.size(); i++ )
2424  {
2425  MB_CHK_SET_ERR( mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) ),
2426  "Failed to get tag value" );
2427  MB_CHK_SET_ERR( mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) ),
2428  "Failed to set tag value on new cell" );
2429  }
2430  }
2431  }
2432 
2433  MB_CHK_SET_ERR( mb->remove_entities( file_set, modifiedCells ), "Failed to remove old cells from file set" );
2434  MB_CHK_SET_ERR( mb->delete_entities( modifiedCells ), "Failed to delete old cells" );
2435  MB_CHK_SET_ERR( mb->add_entities( file_set, newCells ), "Failed to add new cells to file set" );
2436  MB_CHK_SET_ERR( mb->add_entities( file_set, verts ), "Failed to add verts to the file set" );
2437 
2438  return MB_SUCCESS;
2439 }

References moab::Range::begin(), moab::Range::end(), moab::index, moab::Range::insert(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, MBQUAD, and MBTRI.

Referenced by moab::NCHelperDomain::create_mesh(), moab::NCHelperScrip::create_mesh(), and remove_duplicate_vertices().

◆ reverse_gnomonic_projection()

ErrorCode moab::IntxUtils::reverse_gnomonic_projection ( const double &  c1,
const double &  c2,
double  R,
int  plane,
CartVect &  pos 
)
static

Definition at line 615 of file IntxUtils.cpp.

620 {
621 
622  // the new point will be on line beta*pos
623  double len = sqrt( c1 * c1 + c2 * c2 + R * R );
624  double beta = R / len; // it is less than 1, in general
625 
626  switch( plane )
627  {
628  case 1: {
629  // the plane with x = R; x>y, x>z
630  // c1->y, c2->z
631  pos[0] = beta * R;
632  pos[1] = c1 * beta;
633  pos[2] = c2 * beta;
634  break;
635  }
636  case 2: {
637  // y = R -> zx
638  pos[1] = R * beta;
639  pos[2] = c1 * beta;
640  pos[0] = c2 * beta;
641  break;
642  }
643  case 3: {
644  // x=-R, -> yz
645  pos[0] = -R * beta;
646  pos[1] = -c1 * beta; // the sign is to preserve orientation
647  pos[2] = c2 * beta;
648  break;
649  }
650  case 4: {
651  // y = -R
652  pos[1] = -R * beta;
653  pos[2] = -c1 * beta; // the sign is to preserve orientation
654  pos[0] = c2 * beta;
655  break;
656  }
657  case 5: {
658  // z = -R
659  pos[2] = -R * beta;
660  pos[0] = -c1 * beta; // the sign is to preserve orientation
661  pos[1] = c2 * beta;
662  break;
663  }
664  case 6: {
665  pos[2] = R * beta;
666  pos[0] = c1 * beta;
667  pos[1] = c2 * beta;
668  break;
669  }
670  default:
671  return MB_FAILURE; // error
672  }
673 
674  return MB_SUCCESS; // no error
675 }

References MB_SUCCESS, and moab::R.

Referenced by moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and moab::BoundBox::update_box_spherical_elem().

◆ ScaleToRadius()

ErrorCode moab::IntxUtils::ScaleToRadius ( Interface *  mb,
EntityHandle  set,
double  R 
)
static

Definition at line 1117 of file IntxUtils.cpp.

1118 {
1119  Range nodes;
1120  ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
1121  if( rval != moab::MB_SUCCESS ) return rval;
1122 
1123  // one by one, get the node and project it on the sphere, with a radius given
1124  // the center of the sphere is at 0,0,0
1125  for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
1126  {
1127  EntityHandle nd = *nit;
1128  CartVect pos;
1129  rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
1130  if( rval != moab::MB_SUCCESS ) return rval;
1131  double len = pos.length();
1132  if( len == 0. ) return MB_FAILURE;
1133  pos = R / len * pos;
1134  rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
1135  if( rval != moab::MB_SUCCESS ) return rval;
1136  }
1137  return MB_SUCCESS;
1138 }

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::CartVect::length(), mb, MB_SUCCESS, MBVERTEX, and moab::R.

Referenced by CreateTempestMesh(), global_gnomonic_projection(), global_gnomonic_projection_general(), and main().

◆ SortAndRemoveDoubles2()

int moab::IntxUtils::SortAndRemoveDoubles2 ( double *  P,
int &  nP,
double  epsilon_1 
)
static

Sorts and removes duplicate points in the given array.

Note: nP might be modified too, we will remove duplicates if found

Parameters
PThe array of points to be sorted and checked for duplicates.
nPThe number of points in P.
epsilon_1The epsilon value for distance comparison.
Returns
0 if successful.

Definition at line 129 of file IntxUtils.cpp.

130 {
131  if( nP < 2 ) return 0; // nothing to do
132 
133  // center of gravity for the points
134  double c[2] = { 0., 0. };
135  int k = 0;
136  for( k = 0; k < nP; k++ )
137  {
138  c[0] += P[2 * k];
139  c[1] += P[2 * k + 1];
140  }
141  c[0] /= nP;
142  c[1] /= nP;
143 
144  // how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
145  // intersection points
146  struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
147 
148  for( k = 0; k < nP; k++ )
149  {
150  double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
151  if( x != 0. || y != 0. )
152  {
153  pairAngleIndex[k].angle = atan2( y, x );
154  }
155  else
156  {
157  pairAngleIndex[k].angle = 0;
158  // it would mean that the cells are touching at a vertex
159  }
160  pairAngleIndex[k].index = k;
161  }
162 
163  // this should be faster than the bubble sort we had before
164  std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
165  // copy now to a new double array
166  double PCopy[10 * MAXEDGES]; // the same dimension as P; very conservative, but faster than
167  // reallocate for a vector
168  for( k = 0; k < nP; k++ ) // index will show where it should go now;
169  {
170  int ck = pairAngleIndex[k].index;
171  PCopy[2 * k] = P[2 * ck];
172  PCopy[2 * k + 1] = P[2 * ck + 1];
173  }
174  // now copy from PCopy over original P location
175  std::copy( PCopy, PCopy + 2 * nP, P );
176 
177  // eliminate duplicates, finally
178 
179  int i = 0, j = 1; // the next one; j may advance faster than i
180  // check the unit
181  // double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
182  // than 1.e-5 cm
183  while( j < nP )
184  {
185  double d2 = dist2( &P[2 * i], &P[2 * j] );
186  if( d2 > epsilon_1 )
187  {
188  i++;
189  P[2 * i] = P[2 * j];
190  P[2 * i + 1] = P[2 * j + 1];
191  }
192  j++;
193  }
194  // test also the last point with the first one (index 0)
195  // the first one could be at -PI; last one could be at +PI, according to atan2 span
196 
197  double d2 = dist2( P, &P[2 * i] ); // check the first and last points (ordered from -pi to +pi)
198  if( d2 > epsilon_1 )
199  {
200  nP = i + 1;
201  }
202  else
203  nP = i; // effectively delete the last point (that would have been the same with first)
204  if( nP == 0 ) nP = 1; // we should be left with at least one point we already tested if nP is 0 originally
205  return 0;
206 }

References moab::angleAndIndex::angle, moab::angleCompare(), dist2(), moab::angleAndIndex::index, and MAXEDGES.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), and moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

◆ spherical_to_cart()

CartVect moab::IntxUtils::spherical_to_cart ( IntxUtils::SphereCoords &  sc)
static

Definition at line 1108 of file IntxUtils.cpp.

1109 {
1110  CartVect res;
1111  res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
1112  res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
1113  res[2] = sc.R * sin( sc.lat ); // z
1114  return res;
1115 }

References moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

Referenced by main().

◆ transform_coordinates()

void moab::IntxUtils::transform_coordinates ( double *  avg_position,
int  projection_type 
)
static

Definition at line 1016 of file IntxUtils.cpp.

1017 {
1018  if( projection_type == 1 )
1019  {
1020  double R =
1021  avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
1022  R = sqrt( R );
1023  double lat = asin( avg_position[2] / R );
1024  double lon = atan2( avg_position[1], avg_position[0] );
1025  avg_position[0] = lon;
1026  avg_position[1] = lat;
1027  avg_position[2] = R;
1028  }
1029  else if( projection_type == 2 ) // gnomonic projection
1030  {
1031  CartVect pos( avg_position );
1032  int gplane;
1033  IntxUtils::decide_gnomonic_plane( pos, gplane );
1034 
1035  IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
1036  avg_position[2] = 0;
1037  IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
1038  }
1039 }

References decide_gnomonic_plane(), gnomonic_projection(), gnomonic_unroll(), and moab::R.

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