 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#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) |
| static ErrorCode | max_diagonal (Interface *mb, Range cells, int max_edges, double &diagonal) |
Definition at line 16 of file IntxUtils.hpp.
|
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().
|
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().
|
static |
Computes the border points of X in Y2.
| X | The array of points representing X. |
| nX | The number of points in X. |
| Y | The array of points representing Y. |
| nY | The number of points in Y. |
| P | The array to store the border points of X in Y2. |
| side | The array to store the side information for each point in X. |
| epsilon_area | The epsilon value for area comparison. |
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().
|
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().
|
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().
|
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().
|
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().
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.
|
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 }
|
static |
Computes the edge intersections of two elements.
| blue | The array of points representing the blue element. |
| nsBlue | The number of points in the blue element. |
| red | The array of points representing the red element. |
| nsRed | The number of points in the red element. |
| markb | The array to mark the intersecting edges of the blue element. |
| markr | The array to mark the intersecting edges of the red element. |
| points | The array to store the intersection points. |
| nPoints | The number of intersection points found. |
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().
|
static |
Computes the edge intersections between a RLL and CS quad.
Note: Special function.
| blue | The array of points representing the blue element. |
| bluec | The array of Cartesian coordinates for the blue element. |
| blueEdgeType | The array of edge types for the blue element. |
| nsBlue | The number of points in the blue element. |
| red | The array of points representing the red element. |
| redc | The array of Cartesian coordinates for the red element. |
| nsRed | The number of points in the red element. |
| markb | The array to mark the intersecting edges of the blue element. |
| markr | The array to mark the intersecting edges of the red element. |
| plane | The plane of intersection. |
| R | The radius of the sphere. |
| points | The array to store the intersection points. |
| nPoints | The number of intersection points found. |
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().
|
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.
| mb | The interface to the MOAB instance. |
| lset | The handle of the input set containing the polygons. |
| my_rank | The rank of the local process. |
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().
|
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().
|
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.
|
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().
|
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().
|
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().
|
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.
|
static |
Definition at line 2157 of file IntxUtils.cpp.
2158 {
2159 diagonal = 0.0;
2160 std::vector< CartVect > coords( max_edges ); // maximum number of nodes in a cell? hard coded ?
2161 for( auto it = cells.begin(); it != cells.end(); ++it )
2162 {
2163 // get the connectivity, then the coordinates
2164 EntityHandle cell = *it;
2165 const EntityHandle* connec = NULL;
2166 int num_verts = 0;
2167 MB_CHK_SET_ERR( mb->get_connectivity( cell, connec, num_verts ), "Failed to get connectivity" );
2168 MB_CHK_SET_ERR( mb->get_coords( connec, num_verts, &( coords[0][0] ) ), "Failed to get coordinates" );
2169 // compute the max diagonal, in a double loop
2170 for( int i = 0; i < num_verts - 1; i++ )
2171 {
2172 for( int j = i + 1; j < num_verts; j++ )
2173 {
2174 double len_sq = ( coords[i] - coords[j] ).length_squared();
2175 if( len_sq > diagonal ) diagonal = len_sq;
2176 }
2177 }
2178 }
2179
2180 // return the root of the diagonal
2181 diagonal = std::sqrt( diagonal );
2182 return MB_SUCCESS;
2183 }
References moab::Range::begin(), moab::Range::end(), length_squared(), mb, MB_CHK_SET_ERR, and MB_SUCCESS.
|
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().
|
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().
Referenced by iMOAB_MergeVertices().
|
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().
|
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().
|
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().
|
static |
Sorts and removes duplicate points in the given array.
Note: nP might be modified too, we will remove duplicates if found
| P | The array of points to be sorted and checked for duplicates. |
| nP | The number of points in P. |
| epsilon_1 | The epsilon value for distance comparison. |
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().
|
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().
|
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.
1.9.1.