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 47 of file IntxUtils.hpp.

Member Function Documentation

◆ area2D()

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

Definition at line 57 of file IntxUtils.hpp.

58  {
59  // (b-a)x(c-a) / 2
60  return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
61  }

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 2375 of file IntxUtils.cpp.

2384 {
2385  int extraPoints = 0;
2386  // first decide the blue z coordinates
2387  CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
2388  for( int i = 0; i < nsBlue; i++ )
2389  {
2390  if( blueEdgeType[i] == 0 )
2391  {
2392  int iP1 = ( i + 1 ) % nsBlue;
2393  if( bluec[i][2] > bluec[iP1][2] )
2394  {
2395  A = bluec[i];
2396  B = bluec[iP1];
2397  C = bluec[( i + 2 ) % nsBlue];
2398  D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
2399  break;
2400  }
2401  }
2402  }
2403  if( nsBlue == 3 && B[2] < 0 )
2404  {
2405  // select D to be C
2406  D = C;
2407  C = B; // B is the south pole then
2408  }
2409  // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
2410  // AB is const longitude, BC and DA constant latitude
2411  // check now each of the red points if they are inside this rectangle
2412  for( int i = 0; i < nsRed; i++ )
2413  {
2414  const CartVect& X = redc[i];
2415  if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
2416  // now decide if it is between the planes OAB and OCD
2417  if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
2418  {
2419  side[i] = 1; //
2420  // it means point X is in the rectangle that we want , on the sphere
2421  // pass the coords 2d
2422  P[extraPoints * 2] = red2dc[2 * i];
2423  P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
2424  extraPoints++;
2425  }
2426  }
2427  return extraPoints;
2428 }

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

This code defines several utility functions for computing edge intersections and performing geometric operations.

  • borderPointsOfXinY2: Computes the border points of a set of points X inside another set of points Y.
  • SortAndRemoveDoubles2: Sorts a set of points P according to their angles and removes duplicate points.
  • EdgeIntersections2: Computes the intersections between the edges of two sets of points blue and red.
  • EdgeIntxRllCs: Computes the intersections between the edges of a set of points blue and a set of points red on a specific plane.

The code also defines some helper structs and functions used by these utility functions. 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 70 of file IntxUtils.cpp.

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

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 1076 of file IntxUtils.cpp.

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

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 414 of file IntxUtils.cpp.

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

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 2430 of file IntxUtils.cpp.

2431 {
2432  ReadUtilIface* read_iface;
2433  MB_CHK_ERR( mb->query_interface( read_iface ) );
2434  // create the handle tag for the corresponding element / vertex
2435 
2436  EntityHandle dum = 0;
2437  Tag corrTag = 0; // it will be created here
2438  MB_CHK_ERR( mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum ) );
2439 
2440  // give the same global id to new verts and cells created in the lagr(departure) mesh
2441  Tag gid = mb->globalId_tag();
2442 
2443  Range quads;
2444  MB_CHK_ERR( mb->get_entities_by_type( source_set, MBQUAD, quads ) );
2445 
2446  Range connecVerts;
2447  MB_CHK_ERR( mb->get_connectivity( quads, connecVerts ) );
2448 
2449  std::map< EntityHandle, EntityHandle > newNodes;
2450 
2451  std::vector< double* > coords;
2452  EntityHandle start_vert, start_elem, *connect;
2453  int num_verts = connecVerts.size();
2454  MB_CHK_ERR( read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords ) );
2455 
2456  // fill it up
2457  int i = 0;
2458  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
2459  {
2460  EntityHandle oldV = *vit;
2461  CartVect posi;
2462  MB_CHK_ERR( mb->get_coords( &oldV, 1, &( posi[0] ) ) );
2463 
2464  int global_id;
2465  MB_CHK_ERR( mb->tag_get_data( gid, &oldV, 1, &global_id ) );
2466  EntityHandle new_vert = start_vert + i;
2467  // Cppcheck warning (false positive): variable coords is assigned a value that is never used
2468  coords[0][i] = posi[0];
2469  coords[1][i] = posi[1];
2470  coords[2][i] = posi[2];
2471 
2472  newNodes[oldV] = new_vert;
2473  // set also the correspondent tag :)
2474  MB_CHK_ERR( mb->tag_set_data( corrTag, &oldV, 1, &new_vert ) );
2475 
2476  // also the other side
2477  // need to check if we really need this; the new vertex will never need the old vertex
2478  // we have the global id which is the same
2479  MB_CHK_ERR( mb->tag_set_data( corrTag, &new_vert, 1, &oldV ) );
2480  // set the global id on the corresponding vertex the same as the initial vertex
2481  MB_CHK_ERR( mb->tag_set_data( gid, &new_vert, 1, &global_id ) );
2482  }
2483  // now create new quads in order (in a sequence)
2484 
2485  MB_CHK_ERR( read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect ) );
2486 
2487  int ie = 0;
2488  for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
2489  {
2490  EntityHandle q = *it;
2491  int nnodes;
2492  const EntityHandle* conn;
2493  MB_CHK_ERR( mb->get_connectivity( q, conn, nnodes ) );
2494  int global_id;
2495  MB_CHK_ERR( mb->tag_get_data( gid, &q, 1, &global_id ) );
2496 
2497  for( int ii = 0; ii < nnodes; ii++ )
2498  {
2499  EntityHandle v1 = conn[ii];
2500  connect[4 * ie + ii] = newNodes[v1];
2501  }
2502  EntityHandle newElement = start_elem + ie;
2503 
2504  // set the corresponding tag; not sure we need this one, from old to new
2505  MB_CHK_ERR( mb->tag_set_data( corrTag, &q, 1, &newElement ) );
2506  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &q ) );
2507 
2508  // set the global id
2509  MB_CHK_ERR( mb->tag_set_data( gid, &newElement, 1, &global_id ) );
2510 
2511  MB_CHK_ERR( mb->add_entities( dest_set, &newElement, 1 ) );
2512  }
2513 
2514  MB_CHK_ERR( read_iface->update_adjacencies( start_elem, quads.size(), 4, connect ) );
2515 
2516  return MB_SUCCESS;
2517 }

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 51 of file IntxUtils.hpp.

52  {
53  double abx = b[0] - a[0], aby = b[1] - a[1];
54  return sqrt( abx * abx + aby * aby );
55  }

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 1657 of file IntxUtils.cpp.

1658 {
1659  SphereCoords sph1 = cart_to_spherical( p1 );
1660  SphereCoords sph2 = cart_to_spherical( p2 );
1661  // radius should be the same
1662  return sph1.R *
1663  acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
1664 }

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 1966 of file IntxUtils.cpp.

1967 {
1968  return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
1969 }

◆ 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 225 of file IntxUtils.cpp.

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

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 309 of file IntxUtils.cpp.

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

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 1679 of file IntxUtils.cpp.

1680 {
1681  Range inputRange;
1682  MB_CHK_ERR( mb->get_entities_by_dimension( lset, 2, inputRange ) );
1683 
1684  Tag corrTag = nullptr;
1685  EntityHandle dumH = 0;
1686  // no need to check return error
1687  mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
1688 
1689  Tag gidTag = mb->globalId_tag();
1690 
1691  std::vector< double > coords;
1692  coords.resize( 3 * MAXEDGES ); // initial guess only; resized per cell below
1693  // we should create a queue with new polygons that need processing for reflex angles
1694  // (obtuse)
1695  std::queue< EntityHandle > newPolys;
1696  int brokenPolys = 0;
1697  Range::iterator eit = inputRange.begin();
1698  while( eit != inputRange.end() || !newPolys.empty() )
1699  {
1700  EntityHandle eh;
1701  if( eit != inputRange.end() )
1702  {
1703  eh = *eit;
1704  ++eit;
1705  }
1706  else
1707  {
1708  eh = newPolys.front();
1709  newPolys.pop();
1710  }
1711  // get the nodes, then the coordinates
1712  const EntityHandle* verts;
1713  int num_nodes;
1714  MB_CHK_ERR( mb->get_connectivity( eh, verts, num_nodes ) );
1715  int nsides = num_nodes;
1716  // account for possible padded polygons
1717  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1718  nsides--;
1719  EntityHandle corrHandle = 0;
1720  if( corrTag )
1721  {
1722  MB_CHK_ERR( mb->tag_get_data( corrTag, &eh, 1, &corrHandle ) );
1723  }
1724  int gid = 0;
1725  MB_CHK_ERR( mb->tag_get_data( gidTag, &eh, 1, &gid ) );
1726  coords.resize( 3 * nsides );
1727  if( nsides < 4 ) continue; // if already triangles, don't bother
1728  // get coordinates
1729  MB_CHK_ERR( mb->get_coords( verts, nsides, &coords[0] ) );
1730  // compute each angle
1731  bool alreadyBroken = false;
1732 
1733  for( int i = 0; i < nsides; i++ )
1734  {
1735  double* A = &coords[3 * i];
1736  double* B = &coords[3 * ( ( i + 1 ) % nsides )];
1737  double* C = &coords[3 * ( ( i + 2 ) % nsides )];
1738  double angle = IntxUtils::oriented_spherical_angle( A, B, C );
1739  if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
1740  {
1741  if( alreadyBroken )
1742  {
1743  mb->list_entities( &eh, 1 );
1744  mb->list_entities( verts, nsides );
1745  double* D = &coords[3 * ( ( i + 3 ) % nsides )];
1746  std::cout << "ABC: " << angle << " \n";
1747  std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
1748  std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
1749  std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
1750  std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
1751 
1752  return MB_FAILURE;
1753  }
1754  // the bad angle is at i+1;
1755  // create 1 triangle and one polygon; add the polygon to the input range, so
1756  // it will be processed too
1757  // also, add both to the set :) and remove the original polygon from the set
1758  // break the next triangle, even though not optimal
1759  // so create the triangle i+1, i+2, i+3; remove i+2 from original list
1760  // even though not optimal in general, it is good enough.
1761  EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
1762  verts[( i + 3 ) % nsides] };
1763  // create a polygon with num_nodes-1 vertices, and connectivity
1764  // verts[i+1], verts[i+3], (all except i+2)
1765  std::vector< EntityHandle > conn( nsides - 1 );
1766  for( int j = 1; j < nsides; j++ )
1767  {
1768  conn[j - 1] = verts[( i + j + 2 ) % nsides];
1769  }
1770  EntityHandle newElement;
1771  MB_CHK_ERR( mb->create_element( MBTRI, conn3, 3, newElement ) );
1772 
1773  MB_CHK_ERR( mb->add_entities( lset, &newElement, 1 ) );
1774  if( corrTag )
1775  {
1776  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &corrHandle ) );
1777  }
1778  MB_CHK_ERR( mb->tag_set_data( gidTag, &newElement, 1, &gid ) );
1779  if( nsides == 4 )
1780  {
1781  // create another triangle
1782  MB_CHK_ERR( mb->create_element( MBTRI, &conn[0], 3, newElement ) );
1783  }
1784  else
1785  {
1786  // create another polygon, and add it to the inputRange
1787  MB_CHK_ERR( mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement ) );
1788  newPolys.push( newElement ); // because it has less number of edges, the
1789  // reverse should work to find it.
1790  }
1791  MB_CHK_ERR( mb->add_entities( lset, &newElement, 1 ) );
1792  if( corrTag )
1793  {
1794  MB_CHK_ERR( mb->tag_set_data( corrTag, &newElement, 1, &corrHandle ) );
1795  }
1796  MB_CHK_ERR( mb->tag_set_data( gidTag, &newElement, 1, &gid ) );
1797  MB_CHK_ERR( mb->remove_entities( lset, &eh, 1 ) );
1798  brokenPolys++;
1799  alreadyBroken = true; // get out of the loop, element is broken
1800  }
1801  }
1802  }
1803  if( brokenPolys > 0 )
1804  {
1805  std::cout << "on local process " << my_rank << ", " << brokenPolys
1806  << " concave polygons were decomposed in convex ones \n";
1807 #ifdef VERBOSE
1808  std::stringstream fff;
1809  fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
1810  MB_CHK_ERR( mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 ) );
1811  std::cout << "wrote new file set: " << fff.str() << "\n";
1812 #endif
1813  }
1814  return MB_SUCCESS;
1815 }

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 1819 of file IntxUtils.cpp.

1820 {
1821  Range quads;
1822  MB_CHK_ERR( mb->get_entities_by_type( set, MBQUAD, quads ) );
1823 
1824  Tag gid = mb->globalId_tag();
1825  for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
1826  {
1827  EntityHandle quad = *qit;
1828  const EntityHandle* conn4 = nullptr;
1829  int num_nodes = 0;
1830  MB_CHK_ERR( mb->get_connectivity( quad, conn4, num_nodes ) );
1831  for( int i = 0; i < num_nodes; i++ )
1832  {
1833  int next_node_index = ( i + 1 ) % num_nodes;
1834  if( conn4[i] == conn4[next_node_index] )
1835  {
1836  // form a triangle and delete the quad
1837  // first get the global id, to set it on triangle later
1838  int global_id = 0;
1839  MB_CHK_ERR( mb->tag_get_data( gid, &quad, 1, &global_id ) );
1840  int i2 = ( i + 2 ) % num_nodes;
1841  int i3 = ( i + 3 ) % num_nodes;
1842  EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
1843  EntityHandle tri;
1844  MB_CHK_ERR( mb->create_element( MBTRI, conn3, 3, tri ) );
1845  MB_CHK_ERR( mb->add_entities( set, &tri, 1 ) );
1846  MB_CHK_ERR( mb->remove_entities( set, &quad, 1 ) );
1847  MB_CHK_ERR( mb->delete_entities( &quad, 1 ) );
1848  MB_CHK_ERR( mb->tag_set_data( gid, &tri, 1, &global_id ) );
1849  }
1850  }
1851  }
1852  return MB_SUCCESS;
1853 }

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

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), anonymous_namespace{mbtempest.cpp}::convertAndWriteMOABMesh(), 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 844 of file IntxUtils.cpp.

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

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 720 of file IntxUtils.cpp.

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

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 563 of file IntxUtils.cpp.

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

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 546 of file IntxUtils.cpp.

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

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 458 of file IntxUtils.cpp.

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

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 681 of file IntxUtils.cpp.

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

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 2062 of file IntxUtils.cpp.

2069 {
2070  const double distTol = R * 1.e-6;
2071  const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
2072  np = 0; // number of points in intersection
2073  CartVect a( A ), b( B ), c( C ), d( D );
2074  // check input first
2075  double R2 = R * R;
2076  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
2077  fabs( d.length_squared() - R2 ) >
2078  10 * Tolerance )
2079  return MB_FAILURE;
2080 
2081  if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
2082  if( ( c - d ).length_squared() < Tolerance ) // edges are too short
2083  return MB_FAILURE;
2084 
2085  // CD is the const latitude arc
2086  if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
2087  return MB_FAILURE;
2088 
2089  if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
2090 
2091  // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
2092  // great circle arc AB normal to the AB circle:
2093  CartVect n1 = a * b; // the normal to the great circle arc (circle)
2094  // solve the system of equations:
2095  /*
2096  * n1%(x, y, z) = 0 // on the great circle
2097  * z = C[2];
2098  * x^2+y^2+z^2 = R^2
2099  */
2100  double z = C[2];
2101  if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
2102  {
2103  // it is the Equator; check if the const lat edge is Equator too
2104  if( fabs( C[2] ) > distTol )
2105  {
2106  return MB_FAILURE; // no intx, too far from Eq
2107  }
2108  else
2109  {
2110  // all points are on the equator
2111  //
2112  CartVect cd = c * d;
2113  // by convention, c<d, positive is from c to d
2114  // is a or b between c , d?
2115  CartVect ca = c * a;
2116  CartVect ad = a * d;
2117  CartVect cb = c * b;
2118  CartVect bd = b * d;
2119  bool agtc = ( ca % cd >= -Tolerance ); // a>c?
2120  bool dgta = ( ad % cd >= -Tolerance ); // d>a?
2121  bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
2122  bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
2123  if( agtc )
2124  {
2125  if( dgta )
2126  {
2127  // a is for sure a point
2128  E[0] = a[0];
2129  E[1] = a[1];
2130  E[2] = a[2];
2131  np++;
2132  if( bgtc )
2133  {
2134  if( dgtb )
2135  {
2136  // b is also in between c and d
2137  E[3] = b[0];
2138  E[4] = b[1];
2139  E[5] = b[2];
2140  np++;
2141  }
2142  else
2143  {
2144  // then order is c a d b, intx is ad
2145  E[3] = d[0];
2146  E[4] = d[1];
2147  E[5] = d[2];
2148  np++;
2149  }
2150  }
2151  else
2152  {
2153  // b is less than c, so b c a d, intx is ac
2154  E[3] = c[0];
2155  E[4] = c[1];
2156  E[5] = c[2];
2157  np++; // what if E[0] is E[3]?
2158  }
2159  }
2160  else // c < d < a
2161  {
2162  if( dgtb ) // d is for sure in
2163  {
2164  E[0] = d[0];
2165  E[1] = d[1];
2166  E[2] = d[2];
2167  np++;
2168  if( bgtc ) // c<b<d<a
2169  {
2170  // another point is b
2171  E[3] = b[0];
2172  E[4] = b[1];
2173  E[5] = b[2];
2174  np++;
2175  }
2176  else // b<c<d<a
2177  {
2178  // another point is c
2179  E[3] = c[0];
2180  E[4] = c[1];
2181  E[5] = c[2];
2182  np++;
2183  }
2184  }
2185  else
2186  {
2187  // nothing, order is c, d < a, b
2188  }
2189  }
2190  }
2191  else // a < c < d
2192  {
2193  if( bgtc )
2194  {
2195  // c is for sure in
2196  E[0] = c[0];
2197  E[1] = c[1];
2198  E[2] = c[2];
2199  np++;
2200  if( dgtb )
2201  {
2202  // a < c < b < d; second point is b
2203  E[3] = b[0];
2204  E[4] = b[1];
2205  E[5] = b[2];
2206  np++;
2207  }
2208  else
2209  {
2210  // a < c < d < b; second point is d
2211  E[3] = d[0];
2212  E[4] = d[1];
2213  E[5] = d[2];
2214  np++;
2215  }
2216  }
2217  else // a, b < c < d
2218  {
2219  // nothing
2220  }
2221  }
2222  }
2223  // for the 2 points selected, see if it is only one?
2224  // no problem, maybe it will be collapsed later anyway
2225  if( np > 0 ) return MB_SUCCESS;
2226  return MB_FAILURE; // no intersection
2227  }
2228  {
2229  if( fabs( n1[0] ) <= fabs( n1[1] ) )
2230  {
2231  // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
2232  // (u+v*x)^2+x^2=R2-z^2
2233  // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
2234  // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
2235  // x1,2 =
2236  double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
2237  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
2238  double delta = b1 * b1 - 4 * a1 * c1;
2239  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
2240  if( delta > Tolerance ) // 2 solutions possible
2241  {
2242  double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
2243  double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
2244  double y1 = u + v * x1;
2245  double y2 = u + v * x2;
2246  if( verify( a, b, c, d, x1, y1, z ) )
2247  {
2248  E[0] = x1;
2249  E[1] = y1;
2250  E[2] = z;
2251  np++;
2252  }
2253  if( verify( a, b, c, d, x2, y2, z ) )
2254  {
2255  E[3 * np + 0] = x2;
2256  E[3 * np + 1] = y2;
2257  E[3 * np + 2] = z;
2258  np++;
2259  }
2260  }
2261  else
2262  {
2263  // one solution
2264  double x1 = -b1 / 2 / a1;
2265  double y1 = u + v * x1;
2266  if( verify( a, b, c, d, x1, y1, z ) )
2267  {
2268  E[0] = x1;
2269  E[1] = y1;
2270  E[2] = z;
2271  np++;
2272  }
2273  }
2274  }
2275  else
2276  {
2277  // resolve eq in y, reverse
2278  // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
2279  // (u+v*y)^2+y^2 -R2+z^2 =0
2280  // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
2281  //
2282  // x1,2 =
2283  double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
2284  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
2285  double delta = b1 * b1 - 4 * a1 * c1;
2286  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
2287  if( delta > Tolerance ) // 2 solutions possible
2288  {
2289  double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
2290  double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
2291  double x1 = u + v * y1;
2292  double x2 = u + v * y2;
2293  if( verify( a, b, c, d, x1, y1, z ) )
2294  {
2295  E[0] = x1;
2296  E[1] = y1;
2297  E[2] = z;
2298  np++;
2299  }
2300  if( verify( a, b, c, d, x2, y2, z ) )
2301  {
2302  E[3 * np + 0] = x2;
2303  E[3 * np + 1] = y2;
2304  E[3 * np + 2] = z;
2305  np++;
2306  }
2307  }
2308  else
2309  {
2310  // one solution
2311  double y1 = -b1 / 2 / a1;
2312  double x1 = u + v * y1;
2313  if( verify( a, b, c, d, x1, y1, z ) )
2314  {
2315  E[0] = x1;
2316  E[1] = y1;
2317  E[2] = z;
2318  np++;
2319  }
2320  }
2321  }
2322  }
2323 
2324  if( np <= 0 ) return MB_FAILURE;
2325  return MB_SUCCESS;
2326 }

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 1975 of file IntxUtils.cpp.

1976 {
1977  // first verify A, B, C, D are on the same sphere
1978  double R2 = R * R;
1979  const double Tolerance = 1.e-12 * R2;
1980 
1981  CartVect a( A ), b( B ), c( C ), d( D );
1982 
1983  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1984  fabs( d.length_squared() - R2 ) >
1985  10 * Tolerance )
1986  return MB_FAILURE;
1987 
1988  CartVect n1 = a * b;
1989  if( n1.length_squared() < Tolerance ) return MB_FAILURE;
1990 
1991  CartVect n2 = c * d;
1992  if( n2.length_squared() < Tolerance ) return MB_FAILURE;
1993  CartVect n3 = n1 * n2;
1994  n3.normalize();
1995 
1996  n3 = R * n3;
1997  // the intersection is either n3 or -n3
1998  CartVect n4 = a * n3, n5 = n3 * b;
1999  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
2000  {
2001  // n3 is good for ab, see if it is good for cd
2002  n4 = c * n3;
2003  n5 = n3 * d;
2004  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
2005  {
2006  E[0] = n3[0];
2007  E[1] = n3[1];
2008  E[2] = n3[2];
2009  }
2010  else
2011  return MB_FAILURE;
2012  }
2013  else
2014  {
2015  // try -n3
2016  n3 = -n3;
2017  n4 = a * n3, n5 = n3 * b;
2018  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
2019  {
2020  // n3 is good for ab, see if it is good for cd
2021  n4 = c * n3;
2022  n5 = n3 * d;
2023  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
2024  {
2025  E[0] = n3[0];
2026  E[1] = n3[1];
2027  E[2] = n3[2];
2028  }
2029  else
2030  return MB_FAILURE;
2031  }
2032  else
2033  return MB_FAILURE;
2034  }
2035 
2036  return MB_SUCCESS;
2037 }

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 2619 of file IntxUtils.cpp.

2620 {
2621  diagonal = 0.0;
2622  std::vector< CartVect > coords( max_edges ); // maximum number of nodes in a cell? hard coded ?
2623  for( auto it = cells.begin(); it != cells.end(); ++it )
2624  {
2625  // get the connectivity, then the coordinates
2626  EntityHandle cell = *it;
2627  const EntityHandle* connec = nullptr;
2628  int num_verts = 0;
2629  MB_CHK_SET_ERR( mb->get_connectivity( cell, connec, num_verts ), "Failed to get connectivity" );
2630  MB_CHK_SET_ERR( mb->get_coords( connec, num_verts, &( coords[0][0] ) ), "Failed to get coordinates" );
2631  // compute the max diagonal, in a double loop
2632  for( int i = 0; i < num_verts - 1; i++ )
2633  {
2634  for( int j = i + 1; j < num_verts; j++ )
2635  {
2636  double len_sq = ( coords[i] - coords[j] ).length_squared();
2637  if( len_sq > diagonal ) diagonal = len_sq;
2638  }
2639  }
2640  }
2641 
2642  // return the root of the diagonal
2643  diagonal = std::sqrt( diagonal );
2644  return MB_SUCCESS;
2645 }

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 1164 of file IntxUtils.cpp.

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

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 2519 of file IntxUtils.cpp.

2523 {
2524  Range verts;
2525  MB_CHK_ERR( mb->get_entities_by_dimension( file_set, 0, verts ) );
2526  MB_CHK_ERR( mb->remove_entities( file_set, verts ) );
2527 
2528  MergeMesh mm( mb );
2529 
2530  // remove the vertices from the set, before merging
2531 
2532  MB_CHK_ERR( mm.merge_all( file_set, merge_tol ) );
2533 
2534  // now correct vertices that are repeated in polygons
2535  MB_CHK_ERR( remove_padded_vertices( mb, file_set, tagList ) );
2536  return MB_SUCCESS;
2537 }

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 2539 of file IntxUtils.cpp.

2540 {
2541 
2542  // now correct vertices that are repeated in polygons
2543  Range cells;
2544  MB_CHK_ERR( mb->get_entities_by_dimension( file_set, 2, cells ) );
2545 
2546  Range verts;
2547  MB_CHK_ERR( mb->get_connectivity( cells, verts ) );
2548 
2549  Range modifiedCells; // will be deleted at the end; keep the gid
2550  Range newCells;
2551 
2552  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
2553  {
2554  EntityHandle cell = *cit;
2555  const EntityHandle* connec = nullptr;
2556  int num_verts = 0;
2557  MB_CHK_SET_ERR( mb->get_connectivity( cell, connec, num_verts ), "Failed to get connectivity" );
2558 
2559  std::vector< EntityHandle > newConnec;
2560  newConnec.push_back( connec[0] ); // at least one vertex
2561  int index = 0;
2562  int new_size = 1;
2563  while( index < num_verts - 2 )
2564  {
2565  int next_index = ( index + 1 );
2566  if( connec[next_index] != newConnec[new_size - 1] )
2567  {
2568  newConnec.push_back( connec[next_index] );
2569  new_size++;
2570  }
2571  index++;
2572  }
2573  // add the last one only if different from previous and first node
2574  if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
2575  {
2576  newConnec.push_back( connec[num_verts - 1] );
2577  new_size++;
2578  }
2579  if( new_size < num_verts && new_size >= 3 )
2580  {
2581  // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
2582  modifiedCells.insert( cell );
2583  // create a new cell with type triangle, quad or polygon
2584  EntityType type = MBTRI;
2585  if( new_size == 3 )
2586  type = MBTRI;
2587  else if( new_size == 4 )
2588  type = MBQUAD;
2589  else if( new_size > 4 )
2590  type = MBPOLYGON;
2591 
2592  // create new cell
2593  EntityHandle newCell;
2594  MB_CHK_SET_ERR( mb->create_element( type, &newConnec[0], new_size, newCell ), "Failed to create new cell" );
2595  // set the old id to the new element
2596  newCells.insert( newCell );
2597  double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
2598  for( size_t i = 0; i < tagList.size(); i++ )
2599  {
2600  MB_CHK_SET_ERR( mb->tag_get_data( tagList[i], &cell, 1, &value ),
2601  "Failed to get tag value" );
2602  MB_CHK_SET_ERR( mb->tag_set_data( tagList[i], &newCell, 1, &value ),
2603  "Failed to set tag value on new cell" );
2604  }
2605  }
2606  // new_size < 3: degenerate cell collapses below a valid polygon (e.g. SCRIP-style
2607  // padded placeholder with only 2 distinct vertices). Leave the original cell in
2608  // place — its padded connectivity faithfully represents what the input file stored.
2609  }
2610 
2611  MB_CHK_SET_ERR( mb->remove_entities( file_set, modifiedCells ), "Failed to remove old cells from file set" );
2612  MB_CHK_SET_ERR( mb->delete_entities( modifiedCells ), "Failed to delete old cells" );
2613  MB_CHK_SET_ERR( mb->add_entities( file_set, newCells ), "Failed to add new cells to file set" );
2614  MB_CHK_SET_ERR( mb->add_entities( file_set, verts ), "Failed to add verts to the file set" );
2615 
2616  return MB_SUCCESS;
2617 }

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 619 of file IntxUtils.cpp.

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

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 1121 of file IntxUtils.cpp.

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

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

Referenced by 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 133 of file IntxUtils.cpp.

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

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 1112 of file IntxUtils.cpp.

1113 {
1114  CartVect res;
1115  res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
1116  res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
1117  res[2] = sc.R * sin( sc.lat ); // z
1118  return res;
1119 }

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 1020 of file IntxUtils.cpp.

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

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

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