Mesh Oriented datABase  (version 5.5.1)
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 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 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)
 

Detailed Description

Definition at line 16 of file IntxUtils.hpp.

Member Function Documentation

◆ area2D()

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

Definition at line 26 of file IntxUtils.hpp.

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

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), 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 1918 of file IntxUtils.cpp.

1927 {
1928  int extraPoints = 0;
1929  // first decide the blue z coordinates
1930  CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
1931  for( int i = 0; i < nsBlue; i++ )
1932  {
1933  if( blueEdgeType[i] == 0 )
1934  {
1935  int iP1 = ( i + 1 ) % nsBlue;
1936  if( bluec[i][2] > bluec[iP1][2] )
1937  {
1938  A = bluec[i];
1939  B = bluec[iP1];
1940  C = bluec[( i + 2 ) % nsBlue];
1941  D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
1942  break;
1943  }
1944  }
1945  }
1946  if( nsBlue == 3 && B[2] < 0 )
1947  {
1948  // select D to be C
1949  D = C;
1950  C = B; // B is the south pole then
1951  }
1952  // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
1953  // AB is const longitude, BC and DA constant latitude
1954  // check now each of the red points if they are inside this rectangle
1955  for( int i = 0; i < nsRed; i++ )
1956  {
1957  CartVect& X = redc[i];
1958  if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
1959  // now decide if it is between the planes OAB and OCD
1960  if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
1961  {
1962  side[i] = 1; //
1963  // it means point X is in the rectangle that we want , on the sphere
1964  // pass the coords 2d
1965  P[extraPoints * 2] = red2dc[2 * i];
1966  P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
1967  extraPoints++;
1968  }
1969  }
1970  return extraPoints;
1971 }

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

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

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

789 {
790  SphereCoords res;
791  res.R = cart3d.length();
792  if( res.R < 0 )
793  {
794  res.lon = res.lat = 0.;
795  return res;
796  }
797  res.lat = asin( cart3d[2] / res.R );
798  res.lon = atan2( cart3d[1], cart3d[0] );
799  if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
800  // there are some defines it depends on :(
801  // #if defined __USE_BSD || defined __USE_XOPEN ???
802 
803  return res;
804 }

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

Referenced by distance_on_great_circle().

◆ decide_gnomonic_plane()

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

Definition at line 411 of file IntxUtils.cpp.

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

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), global_gnomonic_projection(), 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 1973 of file IntxUtils.cpp.

1974 {
1975  ReadUtilIface* read_iface;
1976  ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
1977  // create the handle tag for the corresponding element / vertex
1978 
1979  EntityHandle dum = 0;
1980  Tag corrTag = 0; // it will be created here
1981  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
1982 
1983  // give the same global id to new verts and cells created in the lagr(departure) mesh
1984  Tag gid = mb->globalId_tag();
1985 
1986  Range quads;
1987  rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
1988 
1989  Range connecVerts;
1990  rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
1991 
1992  std::map< EntityHandle, EntityHandle > newNodes;
1993 
1994  std::vector< double* > coords;
1995  EntityHandle start_vert, start_elem, *connect;
1996  int num_verts = connecVerts.size();
1997  rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
1998  if( MB_SUCCESS != rval ) return rval;
1999  // fill it up
2000  int i = 0;
2001  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
2002  {
2003  EntityHandle oldV = *vit;
2004  CartVect posi;
2005  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
2006 
2007  int global_id;
2008  rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
2009  EntityHandle new_vert = start_vert + i;
2010  // Cppcheck warning (false positive): variable coords is assigned a value that is never used
2011  coords[0][i] = posi[0];
2012  coords[1][i] = posi[1];
2013  coords[2][i] = posi[2];
2014 
2015  newNodes[oldV] = new_vert;
2016  // set also the correspondent tag :)
2017  rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
2018 
2019  // also the other side
2020  // need to check if we really need this; the new vertex will never need the old vertex
2021  // we have the global id which is the same
2022  rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
2023  // set the global id on the corresponding vertex the same as the initial vertex
2024  rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
2025  }
2026  // now create new quads in order (in a sequence)
2027 
2028  rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
2029  if( MB_SUCCESS != rval ) return rval;
2030  int ie = 0;
2031  for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
2032  {
2033  EntityHandle q = *it;
2034  int nnodes;
2035  const EntityHandle* conn;
2036  rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
2037  int global_id;
2038  rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
2039 
2040  for( int ii = 0; ii < nnodes; ii++ )
2041  {
2042  EntityHandle v1 = conn[ii];
2043  connect[4 * ie + ii] = newNodes[v1];
2044  }
2045  EntityHandle newElement = start_elem + ie;
2046 
2047  // set the corresponding tag; not sure we need this one, from old to new
2048  rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
2049  rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
2050 
2051  // set the global id
2052  rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
2053 
2054  rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
2055  }
2056 
2057  rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
2058 
2059  return MB_SUCCESS;
2060 }

References moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), ErrorCode, 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 20 of file IntxUtils.hpp.

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

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

1264 {
1265  SphereCoords sph1 = cart_to_spherical( p1 );
1266  SphereCoords sph2 = cart_to_spherical( p2 );
1267  // radius should be the same
1268  return sph1.R *
1269  acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
1270 }

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

1510 {
1511  return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
1512 }

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

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

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

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

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

1286 {
1287  Range inputRange;
1288  ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
1289 
1290  Tag corrTag = 0;
1291  EntityHandle dumH = 0;
1292  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
1293  if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
1294 
1295  Tag gidTag;
1296  rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
1297 
1298  std::vector< double > coords;
1299  coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
1300  // we should create a queue with new polygons that need processing for reflex angles
1301  // (obtuse)
1302  std::queue< EntityHandle > newPolys;
1303  int brokenPolys = 0;
1304  Range::iterator eit = inputRange.begin();
1305  while( eit != inputRange.end() || !newPolys.empty() )
1306  {
1307  EntityHandle eh;
1308  if( eit != inputRange.end() )
1309  {
1310  eh = *eit;
1311  ++eit;
1312  }
1313  else
1314  {
1315  eh = newPolys.front();
1316  newPolys.pop();
1317  }
1318  // get the nodes, then the coordinates
1319  const EntityHandle* verts;
1320  int num_nodes;
1321  rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
1322  int nsides = num_nodes;
1323  // account for possible padded polygons
1324  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1325  nsides--;
1326  EntityHandle corrHandle = 0;
1327  if( corrTag )
1328  {
1329  rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
1330  }
1331  int gid = 0;
1332  rval = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
1333  coords.resize( 3 * nsides );
1334  if( nsides < 4 ) continue; // if already triangles, don't bother
1335  // get coordinates
1336  rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
1337  // compute each angle
1338  bool alreadyBroken = false;
1339 
1340  for( int i = 0; i < nsides; i++ )
1341  {
1342  double* A = &coords[3 * i];
1343  double* B = &coords[3 * ( ( i + 1 ) % nsides )];
1344  double* C = &coords[3 * ( ( i + 2 ) % nsides )];
1345  double angle = IntxUtils::oriented_spherical_angle( A, B, C );
1346  if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
1347  {
1348  if( alreadyBroken )
1349  {
1350  mb->list_entities( &eh, 1 );
1351  mb->list_entities( verts, nsides );
1352  double* D = &coords[3 * ( ( i + 3 ) % nsides )];
1353  std::cout << "ABC: " << angle << " \n";
1354  std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
1355  std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
1356  std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
1357  std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
1358 
1359  return MB_FAILURE;
1360  }
1361  // the bad angle is at i+1;
1362  // create 1 triangle and one polygon; add the polygon to the input range, so
1363  // it will be processed too
1364  // also, add both to the set :) and remove the original polygon from the set
1365  // break the next triangle, even though not optimal
1366  // so create the triangle i+1, i+2, i+3; remove i+2 from original list
1367  // even though not optimal in general, it is good enough.
1368  EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
1369  verts[( i + 3 ) % nsides] };
1370  // create a polygon with num_nodes-1 vertices, and connectivity
1371  // verts[i+1], verts[i+3], (all except i+2)
1372  std::vector< EntityHandle > conn( nsides - 1 );
1373  for( int j = 1; j < nsides; j++ )
1374  {
1375  conn[j - 1] = verts[( i + j + 2 ) % nsides];
1376  }
1377  EntityHandle newElement;
1378  rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
1379 
1380  rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
1381  if( corrTag )
1382  {
1383  rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
1384  }
1385  rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
1386  if( nsides == 4 )
1387  {
1388  // create another triangle
1389  rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
1390  }
1391  else
1392  {
1393  // create another polygon, and add it to the inputRange
1394  rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
1395  newPolys.push( newElement ); // because it has less number of edges, the
1396  // reverse should work to find it.
1397  }
1398  rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
1399  if( corrTag )
1400  {
1401  rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
1402  }
1403  rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
1404  rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
1405  brokenPolys++;
1406  alreadyBroken = true; // get out of the loop, element is broken
1407  }
1408  }
1409  }
1410  if( brokenPolys > 0 )
1411  {
1412  std::cout << "on local process " << my_rank << ", " << brokenPolys
1413  << " concave polygons were decomposed in convex ones \n";
1414 #ifdef VERBOSE
1415  std::stringstream fff;
1416  fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
1417  rval = mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 );MB_CHK_ERR( rval );
1418  std::cout << "wrote new file set: " << fff.str() << "\n";
1419 #endif
1420  }
1421  return MB_SUCCESS;
1422 }

References moab::angle(), moab::Range::begin(), CORRTAGNAME, moab::Range::end(), ErrorCode, MAXEDGES, mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TAG_NOT_FOUND, MB_TYPE_HANDLE, MB_TYPE_INTEGER, 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 1426 of file IntxUtils.cpp.

1427 {
1428  Range quads;
1429  ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
1430  Tag gid;
1431  gid = mb->globalId_tag();
1432  for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
1433  {
1434  EntityHandle quad = *qit;
1435  const EntityHandle* conn4 = NULL;
1436  int num_nodes = 0;
1437  rval = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
1438  for( int i = 0; i < num_nodes; i++ )
1439  {
1440  int next_node_index = ( i + 1 ) % num_nodes;
1441  if( conn4[i] == conn4[next_node_index] )
1442  {
1443  // form a triangle and delete the quad
1444  // first get the global id, to set it on triangle later
1445  int global_id = 0;
1446  rval = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
1447  int i2 = ( i + 2 ) % num_nodes;
1448  int i3 = ( i + 3 ) % num_nodes;
1449  EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
1450  EntityHandle tri;
1451  rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
1452  mb->add_entities( set, &tri, 1 );
1453  mb->remove_entities( set, &quad, 1 );
1454  mb->delete_entities( &quad, 1 );
1455  rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
1456  }
1457  }
1458  }
1459  return MB_SUCCESS;
1460 }

References moab::Range::begin(), moab::Range::end(), ErrorCode, 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 611 of file IntxUtils.cpp.

616 {
617  std::string parTagName( "PARALLEL_PARTITION" );
618  Tag part_tag;
619  Range partSets;
620  ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
621  if( MB_SUCCESS == rval && part_tag != 0 )
622  {
623  rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
624  }
625  rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
626  // Get all entities of dimension 2
627  Range inputRange; // get
628  rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
629  rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
630 
631  std::map< EntityHandle, int > partsAssign;
632  std::map< int, EntityHandle > newPartSets;
633  if( !partSets.empty() )
634  {
635  // get all cells, and assign parts
636  for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
637  {
638  EntityHandle pSet = *setIt;
639  Range ents;
640  rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
641  int val;
642  rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
643  // create a new set with the same part id tag, in the outSet
644  EntityHandle newPartSet;
645  rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
646  rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
647  newPartSets[val] = newPartSet;
648  rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
649  for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
650  {
651  partsAssign[*it] = val;
652  }
653  }
654  }
655 
656  if( centers_only )
657  {
658  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
659  {
660  CartVect center;
661  EntityHandle cell = *it;
662  rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
663  int plane;
664  decide_gnomonic_plane( center, plane );
665  double c[3];
666  c[2] = 0.;
667  gnomonic_projection( center, R, plane, c[0], c[1] );
668 
669  gnomonic_unroll( c[0], c[1], R, plane );
670 
671  EntityHandle vertex;
672  rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
673  rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
674  }
675  }
676  else
677  {
678  // distribute the cells to 6 planes, based on the center
679  Range subranges[6];
680  for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
681  {
682  CartVect center;
683  EntityHandle cell = *it;
684  rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
685  int plane;
686  decide_gnomonic_plane( center, plane );
687  subranges[plane - 1].insert( cell );
688  }
689  for( int i = 1; i <= 6; i++ )
690  {
691  Range verts;
692  rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
693  std::map< EntityHandle, EntityHandle > corr;
694  for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
695  {
696  CartVect vect;
697  EntityHandle v = *vt;
698  rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
699  double c[3];
700  c[2] = 0.;
701  gnomonic_projection( vect, R, i, c[0], c[1] );
702  gnomonic_unroll( c[0], c[1], R, i );
703  EntityHandle vertex;
704  rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
705  corr[v] = vertex; // for new connectivity
706  }
707  EntityHandle new_conn[20]; // max edges in 2d ?
708  for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
709  {
710  EntityHandle eh = *eit;
711  const EntityHandle* conn = NULL;
712  int num_nodes;
713  rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
714  // build a new vertex array
715  for( int j = 0; j < num_nodes; j++ )
716  new_conn[j] = corr[conn[j]];
717  EntityType type = mb->type_from_handle( eh );
718  EntityHandle newCell;
719  rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
720  rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
721  std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
722  if( mit != partsAssign.end() )
723  {
724  int val = mit->second;
725  rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
726  }
727  }
728  }
729  }
730 
731  return MB_SUCCESS;
732 }

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_SUCCESS, MBENTITYSET, MESHSET_SET, moab::R, ScaleToRadius(), and moab::Interface::UNION.

◆ gnomonic_projection()

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

Definition at line 456 of file IntxUtils.cpp.

457 {
458  double alfa = 1.; // the new point will be on line alfa*pos
459 
460  switch( plane )
461  {
462  case 1: {
463  // the plane with x = R; x>y, x>z
464  // c1->y, c2->z
465  alfa = R / pos[0];
466  c1 = alfa * pos[1];
467  c2 = alfa * pos[2];
468  break;
469  }
470  case 2: {
471  // y = R -> zx
472  alfa = R / pos[1];
473  c1 = alfa * pos[2];
474  c2 = alfa * pos[0];
475  break;
476  }
477  case 3: {
478  // x=-R, -> yz
479  alfa = -R / pos[0];
480  c1 = -alfa * pos[1]; // the sign is to preserve orientation
481  c2 = alfa * pos[2];
482  break;
483  }
484  case 4: {
485  // y = -R
486  alfa = -R / pos[1];
487  c1 = -alfa * pos[2]; // the sign is to preserve orientation
488  c2 = alfa * pos[0];
489  break;
490  }
491  case 5: {
492  // z = -R
493  alfa = -R / pos[2];
494  c1 = -alfa * pos[0]; // the sign is to preserve orientation
495  c2 = alfa * pos[1];
496  break;
497  }
498  case 6: {
499  alfa = R / pos[2];
500  c1 = alfa * pos[0];
501  c2 = alfa * pos[1];
502  break;
503  }
504  default:
505  return MB_FAILURE; // error
506  }
507 
508  return MB_SUCCESS; // no error
509 }

References MB_SUCCESS, and moab::R.

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

◆ gnomonic_unroll()

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

Definition at line 574 of file IntxUtils.cpp.

575 {
576  double tmp;
577  switch( plane )
578  {
579  case 1:
580  break; // nothing
581  case 2: // rotate + 90
582  tmp = c1;
583  c1 = -c2;
584  c2 = tmp;
585  c1 = c1 + 2 * R;
586  break;
587  case 3:
588  c1 = c1 + 4 * R;
589  break;
590  case 4: // rotate with -90 x-> -y; y -> x
591 
592  tmp = c1;
593  c1 = c2;
594  c2 = -tmp;
595  c1 = c1 - 2 * R;
596  break;
597  case 5: // South Pole
598  // rotate 180 then move to (-2, -2)
599  c1 = -c1 - 2. * R;
600  c2 = -c2 - 2. * R;
601  break;
602  ;
603  case 6: // North Pole
604  c1 = c1 - 2. * R;
605  c2 = c2 + 2. * R;
606  break;
607  }
608  return;
609 }

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

1612 {
1613  const double distTol = R * 1.e-6;
1614  const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
1615  np = 0; // number of points in intersection
1616  CartVect a( A ), b( B ), c( C ), d( D );
1617  // check input first
1618  double R2 = R * R;
1619  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1620  fabs( d.length_squared() - R2 ) >
1621  10 * Tolerance )
1622  return MB_FAILURE;
1623 
1624  if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
1625  if( ( c - d ).length_squared() < Tolerance ) // edges are too short
1626  return MB_FAILURE;
1627 
1628  // CD is the const latitude arc
1629  if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
1630  return MB_FAILURE;
1631 
1632  if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
1633 
1634  // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
1635  // great circle arc AB normal to the AB circle:
1636  CartVect n1 = a * b; // the normal to the great circle arc (circle)
1637  // solve the system of equations:
1638  /*
1639  * n1%(x, y, z) = 0 // on the great circle
1640  * z = C[2];
1641  * x^2+y^2+z^2 = R^2
1642  */
1643  double z = C[2];
1644  if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
1645  {
1646  // it is the Equator; check if the const lat edge is Equator too
1647  if( fabs( C[2] ) > distTol )
1648  {
1649  return MB_FAILURE; // no intx, too far from Eq
1650  }
1651  else
1652  {
1653  // all points are on the equator
1654  //
1655  CartVect cd = c * d;
1656  // by convention, c<d, positive is from c to d
1657  // is a or b between c , d?
1658  CartVect ca = c * a;
1659  CartVect ad = a * d;
1660  CartVect cb = c * b;
1661  CartVect bd = b * d;
1662  bool agtc = ( ca % cd >= -Tolerance ); // a>c?
1663  bool dgta = ( ad % cd >= -Tolerance ); // d>a?
1664  bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
1665  bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
1666  if( agtc )
1667  {
1668  if( dgta )
1669  {
1670  // a is for sure a point
1671  E[0] = a[0];
1672  E[1] = a[1];
1673  E[2] = a[2];
1674  np++;
1675  if( bgtc )
1676  {
1677  if( dgtb )
1678  {
1679  // b is also in between c and d
1680  E[3] = b[0];
1681  E[4] = b[1];
1682  E[5] = b[2];
1683  np++;
1684  }
1685  else
1686  {
1687  // then order is c a d b, intx is ad
1688  E[3] = d[0];
1689  E[4] = d[1];
1690  E[5] = d[2];
1691  np++;
1692  }
1693  }
1694  else
1695  {
1696  // b is less than c, so b c a d, intx is ac
1697  E[3] = c[0];
1698  E[4] = c[1];
1699  E[5] = c[2];
1700  np++; // what if E[0] is E[3]?
1701  }
1702  }
1703  else // c < d < a
1704  {
1705  if( dgtb ) // d is for sure in
1706  {
1707  E[0] = d[0];
1708  E[1] = d[1];
1709  E[2] = d[2];
1710  np++;
1711  if( bgtc ) // c<b<d<a
1712  {
1713  // another point is b
1714  E[3] = b[0];
1715  E[4] = b[1];
1716  E[5] = b[2];
1717  np++;
1718  }
1719  else // b<c<d<a
1720  {
1721  // another point is c
1722  E[3] = c[0];
1723  E[4] = c[1];
1724  E[5] = c[2];
1725  np++;
1726  }
1727  }
1728  else
1729  {
1730  // nothing, order is c, d < a, b
1731  }
1732  }
1733  }
1734  else // a < c < d
1735  {
1736  if( bgtc )
1737  {
1738  // c is for sure in
1739  E[0] = c[0];
1740  E[1] = c[1];
1741  E[2] = c[2];
1742  np++;
1743  if( dgtb )
1744  {
1745  // a < c < b < d; second point is b
1746  E[3] = b[0];
1747  E[4] = b[1];
1748  E[5] = b[2];
1749  np++;
1750  }
1751  else
1752  {
1753  // a < c < d < b; second point is d
1754  E[3] = d[0];
1755  E[4] = d[1];
1756  E[5] = d[2];
1757  np++;
1758  }
1759  }
1760  else // a, b < c < d
1761  {
1762  // nothing
1763  }
1764  }
1765  }
1766  // for the 2 points selected, see if it is only one?
1767  // no problem, maybe it will be collapsed later anyway
1768  if( np > 0 ) return MB_SUCCESS;
1769  return MB_FAILURE; // no intersection
1770  }
1771  {
1772  if( fabs( n1[0] ) <= fabs( n1[1] ) )
1773  {
1774  // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
1775  // (u+v*x)^2+x^2=R2-z^2
1776  // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
1777  // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
1778  // x1,2 =
1779  double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
1780  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1781  double delta = b1 * b1 - 4 * a1 * c1;
1782  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1783  if( delta > Tolerance ) // 2 solutions possible
1784  {
1785  double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1786  double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1787  double y1 = u + v * x1;
1788  double y2 = u + v * x2;
1789  if( verify( a, b, c, d, x1, y1, z ) )
1790  {
1791  E[0] = x1;
1792  E[1] = y1;
1793  E[2] = z;
1794  np++;
1795  }
1796  if( verify( a, b, c, d, x2, y2, z ) )
1797  {
1798  E[3 * np + 0] = x2;
1799  E[3 * np + 1] = y2;
1800  E[3 * np + 2] = z;
1801  np++;
1802  }
1803  }
1804  else
1805  {
1806  // one solution
1807  double x1 = -b1 / 2 / a1;
1808  double y1 = u + v * x1;
1809  if( verify( a, b, c, d, x1, y1, z ) )
1810  {
1811  E[0] = x1;
1812  E[1] = y1;
1813  E[2] = z;
1814  np++;
1815  }
1816  }
1817  }
1818  else
1819  {
1820  // resolve eq in y, reverse
1821  // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
1822  // (u+v*y)^2+y^2 -R2+z^2 =0
1823  // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
1824  //
1825  // x1,2 =
1826  double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
1827  double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1828  double delta = b1 * b1 - 4 * a1 * c1;
1829  if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1830  if( delta > Tolerance ) // 2 solutions possible
1831  {
1832  double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1833  double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1834  double x1 = u + v * y1;
1835  double x2 = u + v * y2;
1836  if( verify( a, b, c, d, x1, y1, z ) )
1837  {
1838  E[0] = x1;
1839  E[1] = y1;
1840  E[2] = z;
1841  np++;
1842  }
1843  if( verify( a, b, c, d, x2, y2, z ) )
1844  {
1845  E[3 * np + 0] = x2;
1846  E[3 * np + 1] = y2;
1847  E[3 * np + 2] = z;
1848  np++;
1849  }
1850  }
1851  else
1852  {
1853  // one solution
1854  double y1 = -b1 / 2 / a1;
1855  double x1 = u + v * y1;
1856  if( verify( a, b, c, d, x1, y1, z ) )
1857  {
1858  E[0] = x1;
1859  E[1] = y1;
1860  E[2] = z;
1861  np++;
1862  }
1863  }
1864  }
1865  }
1866 
1867  if( np <= 0 ) return MB_FAILURE;
1868  return MB_SUCCESS;
1869 }

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

1519 {
1520  // first verify A, B, C, D are on the same sphere
1521  double R2 = R * R;
1522  const double Tolerance = 1.e-12 * R2;
1523 
1524  CartVect a( A ), b( B ), c( C ), d( D );
1525 
1526  if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
1527  fabs( d.length_squared() - R2 ) >
1528  10 * Tolerance )
1529  return MB_FAILURE;
1530 
1531  CartVect n1 = a * b;
1532  if( n1.length_squared() < Tolerance ) return MB_FAILURE;
1533 
1534  CartVect n2 = c * d;
1535  if( n2.length_squared() < Tolerance ) return MB_FAILURE;
1536  CartVect n3 = n1 * n2;
1537  n3.normalize();
1538 
1539  n3 = R * n3;
1540  // the intersection is either n3 or -n3
1541  CartVect n4 = a * n3, n5 = n3 * b;
1542  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1543  {
1544  // n3 is good for ab, see if it is good for cd
1545  n4 = c * n3;
1546  n5 = n3 * d;
1547  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1548  {
1549  E[0] = n3[0];
1550  E[1] = n3[1];
1551  E[2] = n3[2];
1552  }
1553  else
1554  return MB_FAILURE;
1555  }
1556  else
1557  {
1558  // try -n3
1559  n3 = -n3;
1560  n4 = a * n3, n5 = n3 * b;
1561  if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
1562  {
1563  // n3 is good for ab, see if it is good for cd
1564  n4 = c * n3;
1565  n5 = n3 * d;
1566  if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
1567  {
1568  E[0] = n3[0];
1569  E[1] = n3[1];
1570  E[2] = n3[2];
1571  }
1572  else
1573  return MB_FAILURE;
1574  }
1575  else
1576  return MB_FAILURE;
1577  }
1578 
1579  return MB_SUCCESS;
1580 }

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

◆ oriented_spherical_angle()

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

Definition at line 876 of file IntxUtils.cpp.

877 {
878  // assume the same radius, sphere at origin
879  CartVect a( A ), b( B ), c( C );
880  CartVect normalOAB = a * b;
881  CartVect normalOCB = c * b;
882  CartVect orient = ( c - b ) * ( a - b );
883  double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
884  if( ang != ang )
885  {
886  // signal of a nan
887  std::cout << a << " " << b << " " << c << "\n";
888  std::cout << ang << "\n";
889  }
890  if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
891 
892  return ang;
893 }

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

2066 {
2067  Range verts;
2068  ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
2069  rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
2070 
2071  MergeMesh mm( mb );
2072 
2073  // remove the vertices from the set, before merging
2074 
2075  rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
2076 
2077  // now correct vertices that are repeated in polygons
2078  rval = remove_padded_vertices( mb, file_set, tagList );
2079  return MB_SUCCESS;
2080 }

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

◆ remove_padded_vertices()

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

Definition at line 2082 of file IntxUtils.cpp.

2083 {
2084 
2085  // now correct vertices that are repeated in polygons
2086  Range cells;
2087  ErrorCode rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
2088 
2089  Range verts;
2090  rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
2091 
2092  Range modifiedCells; // will be deleted at the end; keep the gid
2093  Range newCells;
2094 
2095  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
2096  {
2097  EntityHandle cell = *cit;
2098  const EntityHandle* connec = NULL;
2099  int num_verts = 0;
2100  rval = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
2101 
2102  std::vector< EntityHandle > newConnec;
2103  newConnec.push_back( connec[0] ); // at least one vertex
2104  int index = 0;
2105  int new_size = 1;
2106  while( index < num_verts - 2 )
2107  {
2108  int next_index = ( index + 1 );
2109  if( connec[next_index] != newConnec[new_size - 1] )
2110  {
2111  newConnec.push_back( connec[next_index] );
2112  new_size++;
2113  }
2114  index++;
2115  }
2116  // add the last one only if different from previous and first node
2117  if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
2118  {
2119  newConnec.push_back( connec[num_verts - 1] );
2120  new_size++;
2121  }
2122  if( new_size < num_verts )
2123  {
2124  // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
2125  modifiedCells.insert( cell );
2126  // create a new cell with type triangle, quad or polygon
2127  EntityType type = MBTRI;
2128  if( new_size == 3 )
2129  type = MBTRI;
2130  else if( new_size == 4 )
2131  type = MBQUAD;
2132  else if( new_size > 4 )
2133  type = MBPOLYGON;
2134 
2135  // create new cell
2136  EntityHandle newCell;
2137  rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
2138  // set the old id to the new element
2139  newCells.insert( newCell );
2140  double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
2141  for( size_t i = 0; i < tagList.size(); i++ )
2142  {
2143  rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
2144  rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
2145  }
2146  }
2147  }
2148 
2149  rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
2150  rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
2151  rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
2152  rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
2153 
2154  return MB_SUCCESS;
2155 }

References moab::Range::begin(), moab::Range::end(), ErrorCode, 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 512 of file IntxUtils.cpp.

517 {
518 
519  // the new point will be on line beta*pos
520  double len = sqrt( c1 * c1 + c2 * c2 + R * R );
521  double beta = R / len; // it is less than 1, in general
522 
523  switch( plane )
524  {
525  case 1: {
526  // the plane with x = R; x>y, x>z
527  // c1->y, c2->z
528  pos[0] = beta * R;
529  pos[1] = c1 * beta;
530  pos[2] = c2 * beta;
531  break;
532  }
533  case 2: {
534  // y = R -> zx
535  pos[1] = R * beta;
536  pos[2] = c1 * beta;
537  pos[0] = c2 * beta;
538  break;
539  }
540  case 3: {
541  // x=-R, -> yz
542  pos[0] = -R * beta;
543  pos[1] = -c1 * beta; // the sign is to preserve orientation
544  pos[2] = c2 * beta;
545  break;
546  }
547  case 4: {
548  // y = -R
549  pos[1] = -R * beta;
550  pos[2] = -c1 * beta; // the sign is to preserve orientation
551  pos[0] = c2 * beta;
552  break;
553  }
554  case 5: {
555  // z = -R
556  pos[2] = -R * beta;
557  pos[0] = -c1 * beta; // the sign is to preserve orientation
558  pos[1] = c2 * beta;
559  break;
560  }
561  case 6: {
562  pos[2] = R * beta;
563  pos[0] = c1 * beta;
564  pos[1] = c2 * beta;
565  break;
566  }
567  default:
568  return MB_FAILURE; // error
569  }
570 
571  return MB_SUCCESS; // no error
572 }

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

834 {
835  Range nodes;
836  ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
837  if( rval != moab::MB_SUCCESS ) return rval;
838 
839  // one by one, get the node and project it on the sphere, with a radius given
840  // the center of the sphere is at 0,0,0
841  for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
842  {
843  EntityHandle nd = *nit;
844  CartVect pos;
845  rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
846  if( rval != moab::MB_SUCCESS ) return rval;
847  double len = pos.length();
848  if( len == 0. ) return MB_FAILURE;
849  pos = R / len * pos;
850  rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
851  if( rval != moab::MB_SUCCESS ) return rval;
852  }
853  return MB_SUCCESS;
854 }

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

Referenced by CreateTempestMesh(), global_gnomonic_projection(), 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 130 of file IntxUtils.cpp.

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

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

825 {
826  CartVect res;
827  res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
828  res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
829  res[2] = sc.R * sin( sc.lat ); // z
830  return res;
831 }

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

734 {
735  if( projection_type == 1 )
736  {
737  double R =
738  avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
739  R = sqrt( R );
740  double lat = asin( avg_position[2] / R );
741  double lon = atan2( avg_position[1], avg_position[0] );
742  avg_position[0] = lon;
743  avg_position[1] = lat;
744  avg_position[2] = R;
745  }
746  else if( projection_type == 2 ) // gnomonic projection
747  {
748  CartVect pos( avg_position );
749  int gplane;
750  IntxUtils::decide_gnomonic_plane( pos, gplane );
751 
752  IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
753  avg_position[2] = 0;
754  IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
755  }
756 }

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

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