|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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) |
| static double | oriented_spherical_angle (double *A, double *B, 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) |
Definition at line 16 of file IntxUtils.hpp.
| static double moab::IntxUtils::area2D | ( | double * | a, |
| double * | b, | ||
| double * | c | ||
| ) | [inline, static] |
Definition at line 26 of file IntxUtils.hpp.
Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::IntxRllCssphere::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxAreaUtils::positive_orientation(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), and moab::Intx2MeshOnSphere::setup_tgt_cell().
{
// (b-a)x(c-a) / 2
return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
}
| 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 1733 of file IntxUtils.cpp.
Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().
{
int extraPoints = 0;
// first decide the blue z coordinates
CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
for( int i = 0; i < nsBlue; i++ )
{
if( blueEdgeType[i] == 0 )
{
int iP1 = ( i + 1 ) % nsBlue;
if( bluec[i][2] > bluec[iP1][2] )
{
A = bluec[i];
B = bluec[iP1];
C = bluec[( i + 2 ) % nsBlue];
D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
break;
}
}
}
if( nsBlue == 3 && B[2] < 0 )
{
// select D to be C
D = C;
C = B; // B is the south pole then
}
// so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
// AB is const longitude, BC and DA constant latitude
// check now each of the red points if they are inside this rectangle
for( int i = 0; i < nsRed; i++ )
{
CartVect& X = redc[i];
if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
// now decide if it is between the planes OAB and OCD
if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
{
side[i] = 1; //
// it means point X is in the rectangle that we want , on the sphere
// pass the coords 2d
P[extraPoints * 2] = red2dc[2 * i];
P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
extraPoints++;
}
}
return extraPoints;
}
| int moab::IntxUtils::borderPointsOfXinY2 | ( | double * | X, |
| int | nX, | ||
| double * | Y, | ||
| int | nY, | ||
| double * | P, | ||
| int * | side, | ||
| double | epsilon_area | ||
| ) | [static] |
Definition at line 44 of file IntxUtils.cpp.
Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().
{
// 2 triangles, 3 corners, is the corner of X in Y?
// Y must have a positive area
/*
*/
int extraPoint = 0;
for( int i = 0; i < nX; i++ )
{
// compute double the area of all nY triangles formed by a side of Y and a corner of X; if
// one is negative, stop (negative means it is outside; X and Y are all oriented such that
// they are positive oriented;
// if one area is negative, it means it is outside the convex region, for sure)
double* A = X + 2 * i;
int inside = 1;
for( int j = 0; j < nY; j++ )
{
double* B = Y + 2 * j;
int j1 = ( j + 1 ) % nY;
double* C = Y + 2 * j1; // no copy of data
double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
if( area2 < -epsilon_area )
{
inside = 0;
break;
}
}
if( inside )
{
side[i] = 1; // so vertex i of X is inside the convex region formed by Y
// so side has nX dimension (first array)
P[extraPoint * 2] = A[0];
P[extraPoint * 2 + 1] = A[1];
extraPoint++;
}
}
return extraPoint;
}
| IntxUtils::SphereCoords moab::IntxUtils::cart_to_spherical | ( | CartVect & | cart3d | ) | [static] |
Definition at line 726 of file IntxUtils.cpp.
References moab::IntxUtils::SphereCoords::lat, moab::CartVect::length(), moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.
Referenced by add_field_value(), IntxUtilsCSLAM::departure_point_case1(), departure_point_swirl(), departure_point_swirl_rot(), distance_on_great_circle(), manufacture_lagrange_mesh_on_sphere(), set_density(), and IntxUtilsCSLAM::velocity_case1().
{
SphereCoords res;
res.R = cart3d.length();
if( res.R < 0 )
{
res.lon = res.lat = 0.;
return res;
}
res.lat = asin( cart3d[2] / res.R );
res.lon = atan2( cart3d[1], cart3d[0] );
if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
// there are some defines it depends on :(
// #if defined __USE_BSD || defined __USE_XOPEN ???
return res;
}
| void moab::IntxUtils::decide_gnomonic_plane | ( | const CartVect & | pos, |
| int & | oPlane | ||
| ) | [static] |
Definition at line 349 of file IntxUtils.cpp.
Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), get_gnomonic_plane(), global_gnomonic_projection(), moab::IntxRllCssphere::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().
{
// decide plane, based on max x, y, z
if( fabs( pos[0] ) < fabs( pos[1] ) )
{
if( fabs( pos[2] ) < fabs( pos[1] ) )
{
// pos[1] is biggest
if( pos[1] > 0 )
plane = 2;
else
plane = 4;
}
else
{
// pos[2] is biggest
if( pos[2] < 0 )
plane = 5;
else
plane = 6;
}
}
else
{
if( fabs( pos[2] ) < fabs( pos[0] ) )
{
// pos[0] is the greatest
if( pos[0] > 0 )
plane = 1;
else
plane = 3;
}
else
{
// pos[2] is biggest
if( pos[2] < 0 )
plane = 5;
else
plane = 6;
}
}
return;
}
| ErrorCode moab::IntxUtils::deep_copy_set_with_quads | ( | Interface * | mb, |
| EntityHandle | source_set, | ||
| EntityHandle | dest_set | ||
| ) | [static] |
Definition at line 1788 of file IntxUtils.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::ReadUtilIface::get_element_connect(), moab::Interface::get_entities_by_type(), moab::ReadUtilIface::get_node_coords(), moab::Interface::globalId_tag(), ie, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, MBQUAD, moab::Interface::query_interface(), moab::Range::size(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::ReadUtilIface::update_adjacencies().
{
ReadUtilIface* read_iface;
ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
// create the handle tag for the corresponding element / vertex
EntityHandle dum = 0;
Tag corrTag = 0; // it will be created here
rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
// give the same global id to new verts and cells created in the lagr(departure) mesh
Tag gid = mb->globalId_tag();
Range quads;
rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
Range connecVerts;
rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
std::map< EntityHandle, EntityHandle > newNodes;
std::vector< double* > coords;
EntityHandle start_vert, start_elem, *connect;
int num_verts = connecVerts.size();
rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
if( MB_SUCCESS != rval ) return rval;
// fill it up
int i = 0;
for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
{
EntityHandle oldV = *vit;
CartVect posi;
rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
int global_id;
rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
EntityHandle new_vert = start_vert + i;
// Cppcheck warning (false positive): variable coords is assigned a value that is never used
coords[0][i] = posi[0];
coords[1][i] = posi[1];
coords[2][i] = posi[2];
newNodes[oldV] = new_vert;
// set also the correspondent tag :)
rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
// also the other side
// need to check if we really need this; the new vertex will never need the old vertex
// we have the global id which is the same
rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
// set the global id on the corresponding vertex the same as the initial vertex
rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
}
// now create new quads in order (in a sequence)
rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
if( MB_SUCCESS != rval ) return rval;
int ie = 0;
for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
{
EntityHandle q = *it;
int nnodes;
const EntityHandle* conn;
rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
int global_id;
rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
for( int ii = 0; ii < nnodes; ii++ )
{
EntityHandle v1 = conn[ii];
connect[4 * ie + ii] = newNodes[v1];
}
EntityHandle newElement = start_elem + ie;
// set the corresponding tag; not sure we need this one, from old to new
rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
// set the global id
rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
}
rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
return MB_SUCCESS;
}
| static double moab::IntxUtils::dist2 | ( | double * | a, |
| double * | b | ||
| ) | [inline, static] |
Definition at line 20 of file IntxUtils.hpp.
Referenced by moab::Intx2MeshInPlane::findNodes(), moab::IntxRllCssphere::findNodes(), moab::Intx2MeshOnSphere::findNodes(), and SortAndRemoveDoubles2().
{
double abx = b[0] - a[0], aby = b[1] - a[1];
return sqrt( abx * abx + aby * aby );
}
| double moab::IntxUtils::distance_on_great_circle | ( | CartVect & | p1, |
| CartVect & | p2 | ||
| ) | [static] |
Definition at line 1075 of file IntxUtils.cpp.
References cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.
{
SphereCoords sph1 = cart_to_spherical( p1 );
SphereCoords sph2 = cart_to_spherical( p2 );
// radius should be the same
return sph1.R *
acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
}
| double moab::IntxUtils::distance_on_sphere | ( | double | la1, |
| double | te1, | ||
| double | la2, | ||
| double | te2 | ||
| ) | [static] |
Definition at line 1324 of file IntxUtils.cpp.
Referenced by IntxUtilsCSLAM::quasi_smooth_field(), and IntxUtilsCSLAM::slotted_cylinder_field().
{
return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
}
| ErrorCode moab::IntxUtils::EdgeIntersections2 | ( | double * | blue, |
| int | nsBlue, | ||
| double * | red, | ||
| int | nsRed, | ||
| int * | markb, | ||
| int * | markr, | ||
| double * | points, | ||
| int & | nPoints | ||
| ) | [static] |
Definition at line 179 of file IntxUtils.cpp.
References MAXEDGES, and MB_SUCCESS.
Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().
{
/* EDGEINTERSECTIONS computes edge intersections of two elements
[P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
and blue ( stored column wise )
(point coordinates are stored column-wise, in counter clock
order) the points P where their edges intersect. In addition,
in n the indices of which neighbors of red are also intersecting
with blue are given.
*/
// points is an array with enough slots (24 * 2 doubles)
nPoints = 0;
for( int i = 0; i < MAXEDGES; i++ )
{
markb[i] = markr[i] = 0;
}
for( int i = 0; i < nsBlue; i++ )
{
for( int j = 0; j < nsRed; j++ )
{
double b[2];
double a[2][2]; // 2*2
int iPlus1 = ( i + 1 ) % nsBlue;
int jPlus1 = ( j + 1 ) % nsRed;
for( int k = 0; k < 2; k++ )
{
b[k] = red[2 * j + k] - blue[2 * i + k];
// row k of a: a(k, 0), a(k, 1)
a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
}
double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
if( fabs( delta ) > 1.e-14 ) // this is close to machine epsilon
{
// not parallel
double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
{
// the intersection is good
for( int k = 0; k < 2; k++ )
{
points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
}
markb[i] = 1; // so neighbor number i of blue will be considered too.
markr[j] = 1; // this will be used in advancing red around blue quad
nPoints++;
}
}
// the case delta ~ 0. will be considered by the interior points logic
}
}
return MB_SUCCESS;
}
| 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 | Radius, | ||
| double * | points, | ||
| int & | nPoints | ||
| ) | [static] |
Definition at line 244 of file IntxUtils.cpp.
References moab::CartVect::array(), moab::E, gnomonic_projection(), intersect_great_circle_arc_with_clat_arc(), MB_SUCCESS, np, and moab::R.
Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().
{
// if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
// everything else the same (except if there are 2 points resulting, which is rare)
for( int i = 0; i < 4; i++ )
{ // always at most 4 , so maybe don't bother
markb[i] = markr[i] = 0;
}
for( int i = 0; i < nsBlue; i++ )
{
int iPlus1 = ( i + 1 ) % nsBlue;
if( blueEdgeType[i] == 0 ) // old style, just 2d
{
for( int j = 0; j < nsRed; j++ )
{
double b[2];
double a[2][2]; // 2*2
int jPlus1 = ( j + 1 ) % nsRed;
for( int k = 0; k < 2; k++ )
{
b[k] = red[2 * j + k] - blue[2 * i + k];
// row k of a: a(k, 0), a(k, 1)
a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
}
double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
if( fabs( delta ) > 1.e-14 )
{
// not parallel
double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
{
// the intersection is good
for( int k = 0; k < 2; k++ )
{
points[2 * nPoints + k] =
blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
}
markb[i] = 1; // so neighbor number i of blue will be considered too.
markr[j] = 1; // this will be used in advancing red around blue quad
nPoints++;
}
} // if the edges are too "parallel", skip them
}
}
else // edge type is 1, so use 3d intersection
{
CartVect& C = bluec[i];
CartVect& D = bluec[iPlus1];
for( int j = 0; j < nsRed; j++ )
{
int jPlus1 = ( j + 1 ) % nsRed; // nsRed is just 4, forget about it, usually
CartVect& A = redc[j];
CartVect& B = redc[jPlus1];
int np = 0;
double E[9];
intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
if( np == 0 ) continue;
if( np >= 2 )
{
std::cout << "intersection with 2 points :" << A << B << C << D << "\n";
}
for( int k = 0; k < np; k++ )
{
gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
points[2 * nPoints + 1] );
nPoints++;
}
markb[i] = 1; // so neighbor number i of blue will be considered too.
markr[j] = 1; // this will be used in advancing red around blue quad
}
}
}
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::enforce_convexity | ( | Interface * | mb, |
| EntityHandle | set, | ||
| int | rank = 0 |
||
| ) | [static] |
Definition at line 1087 of file IntxUtils.cpp.
References moab::Interface::add_entities(), moab::angle(), moab::Range::begin(), CORRTAGNAME, moab::Interface::create_element(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::id_from_handle(), moab::Interface::list_entities(), MAXEDGES, MB_CHK_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TAG_NOT_FOUND, MB_TYPE_HANDLE, MB_TYPE_INTEGER, MBPOLYGON, MBTRI, oriented_spherical_angle(), moab::Interface::remove_entities(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::Interface::write_file().
Referenced by main(), test_intx_mpas(), and update_tracer().
{
// look at each polygon; compute all angles; if one is reflex, break that angle with
// the next triangle; put the 2 new polys in the set;
// still look at the next poly
// replace it with 2 triangles, and remove from set;
// it should work for all polygons / tested first for case 1, with dt 0.5 (too much deformation)
// get all entities of dimension 2
// then get the connectivity, etc
Range inputRange;
ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
Tag corrTag = 0;
EntityHandle dumH = 0;
rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
Tag gidTag;
rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
std::vector< double > coords;
coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
// we should create a queue with new polygons that need processing for reflex angles
// (obtuse)
std::queue< EntityHandle > newPolys;
int brokenPolys = 0;
Range::iterator eit = inputRange.begin();
while( eit != inputRange.end() || !newPolys.empty() )
{
EntityHandle eh;
if( eit != inputRange.end() )
{
eh = *eit;
++eit;
}
else
{
eh = newPolys.front();
newPolys.pop();
}
// get the nodes, then the coordinates
const EntityHandle* verts;
int num_nodes;
rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
int nsides = num_nodes;
// account for possible padded polygons
while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
nsides--;
EntityHandle corrHandle = 0;
if( corrTag )
{
rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
}
int gid = 0;
rval = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
coords.resize( 3 * nsides );
if( nsides < 4 ) continue; // if already triangles, don't bother
// get coordinates
rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
// compute each angle
bool alreadyBroken = false;
for( int i = 0; i < nsides; i++ )
{
double* A = &coords[3 * i];
double* B = &coords[3 * ( ( i + 1 ) % nsides )];
double* C = &coords[3 * ( ( i + 2 ) % nsides )];
double angle = IntxUtils::oriented_spherical_angle( A, B, C );
if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
{
if( alreadyBroken )
{
mb->list_entities( &eh, 1 );
mb->list_entities( verts, nsides );
double* D = &coords[3 * ( ( i + 3 ) % nsides )];
std::cout << "ABC: " << angle << " \n";
std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
return MB_FAILURE;
}
// the bad angle is at i+1;
// create 1 triangle and one polygon; add the polygon to the input range, so
// it will be processed too
// also, add both to the set :) and remove the original polygon from the set
// break the next triangle, even though not optimal
// so create the triangle i+1, i+2, i+3; remove i+2 from original list
// even though not optimal in general, it is good enough.
EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
verts[( i + 3 ) % nsides] };
// create a polygon with num_nodes-1 vertices, and connectivity
// verts[i+1], verts[i+3], (all except i+2)
std::vector< EntityHandle > conn( nsides - 1 );
for( int j = 1; j < nsides; j++ )
{
conn[j - 1] = verts[( i + j + 2 ) % nsides];
}
EntityHandle newElement;
rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
if( corrTag )
{
rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
}
rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
if( nsides == 4 )
{
// create another triangle
rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
}
else
{
// create another polygon, and add it to the inputRange
rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
newPolys.push( newElement ); // because it has less number of edges, the
// reverse should work to find it.
}
rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
if( corrTag )
{
rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
}
rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
brokenPolys++;
alreadyBroken = true; // get out of the loop, element is broken
}
}
}
if( brokenPolys > 0 )
{
std::cout << "on local process " << my_rank << ", " << brokenPolys
<< " concave polygons were decomposed in convex ones \n";
#ifdef VERBOSE
std::stringstream fff;
fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
rval = mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 );MB_CHK_ERR( rval );
std::cout << "wrote new file set: " << fff.str() << "\n";
#endif
}
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::fix_degenerate_quads | ( | Interface * | mb, |
| EntityHandle | set | ||
| ) | [static] |
Definition at line 1236 of file IntxUtils.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_element(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_type(), moab::Interface::globalId_tag(), MB_CHK_ERR, MB_SUCCESS, MBQUAD, MBTRI, moab::Interface::remove_entities(), moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().
Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().
{
Range quads;
ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
Tag gid;
gid = mb->globalId_tag();
for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
{
EntityHandle quad = *qit;
const EntityHandle* conn4 = NULL;
int num_nodes = 0;
rval = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
for( int i = 0; i < num_nodes; i++ )
{
int next_node_index = ( i + 1 ) % num_nodes;
if( conn4[i] == conn4[next_node_index] )
{
// form a triangle and delete the quad
// first get the global id, to set it on triangle later
int global_id = 0;
rval = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
int i2 = ( i + 2 ) % num_nodes;
int i3 = ( i + 3 ) % num_nodes;
EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
EntityHandle tri;
rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
mb->add_entities( set, &tri, 1 );
mb->remove_entities( set, &quad, 1 );
mb->delete_entities( &quad, 1 );
rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
}
}
}
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::global_gnomonic_projection | ( | Interface * | mb, |
| EntityHandle | inSet, | ||
| double | R, | ||
| bool | centers_only, | ||
| EntityHandle & | outSet | ||
| ) | [static] |
Definition at line 549 of file IntxUtils.cpp.
References moab::Interface::add_entities(), moab::CartVect::array(), moab::Range::begin(), center(), moab::Interface::create_element(), moab::Interface::create_meshset(), moab::Interface::create_vertex(), decide_gnomonic_plane(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_handle(), moab::Interface::get_entities_by_type_and_tag(), gnomonic_projection(), gnomonic_unroll(), moab::Range::insert(), MB_CHK_ERR, MB_SUCCESS, MBENTITYSET, MESHSET_SET, ScaleToRadius(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), moab::Interface::type_from_handle(), and moab::Interface::UNION.
Referenced by main().
{
std::string parTagName( "PARALLEL_PARTITION" );
Tag part_tag;
Range partSets;
ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
if( MB_SUCCESS == rval && part_tag != 0 )
{
rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
}
rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
// Get all entities of dimension 2
Range inputRange; // get
rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
std::map< EntityHandle, int > partsAssign;
std::map< int, EntityHandle > newPartSets;
if( !partSets.empty() )
{
// get all cells, and assign parts
for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
{
EntityHandle pSet = *setIt;
Range ents;
rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
int val;
rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
// create a new set with the same part id tag, in the outSet
EntityHandle newPartSet;
rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
newPartSets[val] = newPartSet;
rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
{
partsAssign[*it] = val;
}
}
}
if( centers_only )
{
for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
{
CartVect center;
EntityHandle cell = *it;
rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
int plane;
decide_gnomonic_plane( center, plane );
double c[3];
c[2] = 0.;
gnomonic_projection( center, R, plane, c[0], c[1] );
gnomonic_unroll( c[0], c[1], R, plane );
EntityHandle vertex;
rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
}
}
else
{
// distribute the cells to 6 planes, based on the center
Range subranges[6];
for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
{
CartVect center;
EntityHandle cell = *it;
rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
int plane;
decide_gnomonic_plane( center, plane );
subranges[plane - 1].insert( cell );
}
for( int i = 1; i <= 6; i++ )
{
Range verts;
rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
std::map< EntityHandle, EntityHandle > corr;
for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
{
CartVect vect;
EntityHandle v = *vt;
rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
double c[3];
c[2] = 0.;
gnomonic_projection( vect, R, i, c[0], c[1] );
gnomonic_unroll( c[0], c[1], R, i );
EntityHandle vertex;
rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
corr[v] = vertex; // for new connectivity
}
EntityHandle new_conn[20]; // max edges in 2d ?
for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
{
EntityHandle eh = *eit;
const EntityHandle* conn = NULL;
int num_nodes;
rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
// build a new vertex array
for( int j = 0; j < num_nodes; j++ )
new_conn[j] = corr[conn[j]];
EntityType type = mb->type_from_handle( eh );
EntityHandle newCell;
rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
if( mit != partsAssign.end() )
{
int val = mit->second;
rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
}
}
}
}
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::gnomonic_projection | ( | const CartVect & | pos, |
| double | R, | ||
| int | plane, | ||
| double & | c1, | ||
| double & | c2 | ||
| ) | [static] |
Definition at line 394 of file IntxUtils.cpp.
References MB_SUCCESS.
Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), EdgeIntxRllCs(), get_barycenters(), get_intersection_weights(), get_linear_reconstruction(), global_gnomonic_projection(), moab::IntxRllCssphere::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), test_linear_reconstruction(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().
{
double alfa = 1.; // the new point will be on line alfa*pos
switch( plane )
{
case 1: {
// the plane with x = R; x>y, x>z
// c1->y, c2->z
alfa = R / pos[0];
c1 = alfa * pos[1];
c2 = alfa * pos[2];
break;
}
case 2: {
// y = R -> zx
alfa = R / pos[1];
c1 = alfa * pos[2];
c2 = alfa * pos[0];
break;
}
case 3: {
// x=-R, -> yz
alfa = -R / pos[0];
c1 = -alfa * pos[1]; // the sign is to preserve orientation
c2 = alfa * pos[2];
break;
}
case 4: {
// y = -R
alfa = -R / pos[1];
c1 = -alfa * pos[2]; // the sign is to preserve orientation
c2 = alfa * pos[0];
break;
}
case 5: {
// z = -R
alfa = -R / pos[2];
c1 = -alfa * pos[0]; // the sign is to preserve orientation
c2 = alfa * pos[1];
break;
}
case 6: {
alfa = R / pos[2];
c1 = alfa * pos[0];
c2 = alfa * pos[1];
break;
}
default:
return MB_FAILURE; // error
}
return MB_SUCCESS; // no error
}
| void moab::IntxUtils::gnomonic_unroll | ( | double & | c1, |
| double & | c2, | ||
| double | R, | ||
| int | plane | ||
| ) | [static] |
Definition at line 512 of file IntxUtils.cpp.
References moab::R.
Referenced by global_gnomonic_projection(), and transform_coordinates().
{
double tmp;
switch( plane )
{
case 1:
break; // nothing
case 2: // rotate + 90
tmp = c1;
c1 = -c2;
c2 = tmp;
c1 = c1 + 2 * R;
break;
case 3:
c1 = c1 + 4 * R;
break;
case 4: // rotate with -90 x-> -y; y -> x
tmp = c1;
c1 = c2;
c2 = -tmp;
c1 = c1 - 2 * R;
break;
case 5: // South Pole
// rotate 180 then move to (-2, -2)
c1 = -c1 - 2. * R;
c2 = -c2 - 2. * R;
break;
;
case 6: // North Pole
c1 = c1 - 2. * R;
c2 = c2 + 2. * R;
break;
}
return;
}
| 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 1420 of file IntxUtils.cpp.
References moab::CartVect::length_squared(), MB_SUCCESS, moab::R, and moab::verify().
Referenced by EdgeIntxRllCs(), and test_great_arc_clat_intx().
{
const double distTol = R * 1.e-6;
const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
np = 0; // number of points in intersection
CartVect a( A ), b( B ), c( C ), d( D );
// check input first
double R2 = R * R;
if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
fabs( d.length_squared() - R2 ) >
10 * Tolerance )
return MB_FAILURE;
if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
if( ( c - d ).length_squared() < Tolerance ) // edges are too short
return MB_FAILURE;
// CD is the const latitude arc
if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
return MB_FAILURE;
if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
// find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
// great circle arc AB normal to the AB circle:
CartVect n1 = a * b; // the normal to the great circle arc (circle)
// solve the system of equations:
/*
* n1%(x, y, z) = 0 // on the great circle
* z = C[2];
* x^2+y^2+z^2 = R^2
*/
double z = C[2];
if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
{
// it is the Equator; check if the const lat edge is Equator too
if( fabs( C[2] ) > distTol )
{
return MB_FAILURE; // no intx, too far from Eq
}
else
{
// all points are on the equator
//
CartVect cd = c * d;
// by convention, c<d, positive is from c to d
// is a or b between c , d?
CartVect ca = c * a;
CartVect ad = a * d;
CartVect cb = c * b;
CartVect bd = b * d;
bool agtc = ( ca % cd >= -Tolerance ); // a>c?
bool dgta = ( ad % cd >= -Tolerance ); // d>a?
bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
if( agtc )
{
if( dgta )
{
// a is for sure a point
E[0] = a[0];
E[1] = a[1];
E[2] = a[2];
np++;
if( bgtc )
{
if( dgtb )
{
// b is also in between c and d
E[3] = b[0];
E[4] = b[1];
E[5] = b[2];
np++;
}
else
{
// then order is c a d b, intx is ad
E[3] = d[0];
E[4] = d[1];
E[5] = d[2];
np++;
}
}
else
{
// b is less than c, so b c a d, intx is ac
E[3] = c[0];
E[4] = c[1];
E[5] = c[2];
np++; // what if E[0] is E[3]?
}
}
else // c < d < a
{
if( dgtb ) // d is for sure in
{
E[0] = d[0];
E[1] = d[1];
E[2] = d[2];
np++;
if( bgtc ) // c<b<d<a
{
// another point is b
E[3] = b[0];
E[4] = b[1];
E[5] = b[2];
np++;
}
else // b<c<d<a
{
// another point is c
E[3] = c[0];
E[4] = c[1];
E[5] = c[2];
np++;
}
}
else
{
// nothing, order is c, d < a, b
}
}
}
else // a < c < d
{
if( bgtc )
{
// c is for sure in
E[0] = c[0];
E[1] = c[1];
E[2] = c[2];
np++;
if( dgtb )
{
// a < c < b < d; second point is b
E[3] = b[0];
E[4] = b[1];
E[5] = b[2];
np++;
}
else
{
// a < c < d < b; second point is d
E[3] = d[0];
E[4] = d[1];
E[5] = d[2];
np++;
}
}
else // a, b < c < d
{
// nothing
}
}
}
// for the 2 points selected, see if it is only one?
// no problem, maybe it will be collapsed later anyway
if( np > 0 ) return MB_SUCCESS;
return MB_FAILURE; // no intersection
}
{
if( fabs( n1[0] ) <= fabs( n1[1] ) )
{
// resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
// (u+v*x)^2+x^2=R2-z^2
// (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
// delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
// x1,2 =
double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
double delta = b1 * b1 - 4 * a1 * c1;
if( delta < -Tolerance ) return MB_FAILURE; // no intersection
if( delta > Tolerance ) // 2 solutions possible
{
double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
double y1 = u + v * x1;
double y2 = u + v * x2;
if( verify( a, b, c, d, x1, y1, z ) )
{
E[0] = x1;
E[1] = y1;
E[2] = z;
np++;
}
if( verify( a, b, c, d, x2, y2, z ) )
{
E[3 * np + 0] = x2;
E[3 * np + 1] = y2;
E[3 * np + 2] = z;
np++;
}
}
else
{
// one solution
double x1 = -b1 / 2 / a1;
double y1 = u + v * x1;
if( verify( a, b, c, d, x1, y1, z ) )
{
E[0] = x1;
E[1] = y1;
E[2] = z;
np++;
}
}
}
else
{
// resolve eq in y, reverse
// n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
// (u+v*y)^2+y^2 -R2+z^2 =0
// (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
//
// x1,2 =
double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
double delta = b1 * b1 - 4 * a1 * c1;
if( delta < -Tolerance ) return MB_FAILURE; // no intersection
if( delta > Tolerance ) // 2 solutions possible
{
double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
double x1 = u + v * y1;
double x2 = u + v * y2;
if( verify( a, b, c, d, x1, y1, z ) )
{
E[0] = x1;
E[1] = y1;
E[2] = z;
np++;
}
if( verify( a, b, c, d, x2, y2, z ) )
{
E[3 * np + 0] = x2;
E[3 * np + 1] = y2;
E[3 * np + 2] = z;
np++;
}
}
else
{
// one solution
double y1 = -b1 / 2 / a1;
double x1 = u + v * y1;
if( verify( a, b, c, d, x1, y1, z ) )
{
E[0] = x1;
E[1] = y1;
E[2] = z;
np++;
}
}
}
}
if( np <= 0 ) return MB_FAILURE;
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::intersect_great_circle_arcs | ( | double * | A, |
| double * | B, | ||
| double * | C, | ||
| double * | D, | ||
| double | R, | ||
| double * | E | ||
| ) | [static] |
Definition at line 1333 of file IntxUtils.cpp.
References moab::CartVect::length_squared(), MB_SUCCESS, moab::CartVect::normalize(), and moab::R.
Referenced by test_great_arc_intx().
{
// first verify A, B, C, D are on the same sphere
double R2 = R * R;
const double Tolerance = 1.e-12 * R2;
CartVect a( A ), b( B ), c( C ), d( D );
if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
fabs( d.length_squared() - R2 ) >
10 * Tolerance )
return MB_FAILURE;
CartVect n1 = a * b;
if( n1.length_squared() < Tolerance ) return MB_FAILURE;
CartVect n2 = c * d;
if( n2.length_squared() < Tolerance ) return MB_FAILURE;
CartVect n3 = n1 * n2;
n3.normalize();
n3 = R * n3;
// the intersection is either n3 or -n3
CartVect n4 = a * n3, n5 = n3 * b;
if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
{
// n3 is good for ab, see if it is good for cd
n4 = c * n3;
n5 = n3 * d;
if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
{
E[0] = n3[0];
E[1] = n3[1];
E[2] = n3[2];
}
else
return MB_FAILURE;
}
else
{
// try -n3
n3 = -n3;
n4 = a * n3, n5 = n3 * b;
if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
{
// n3 is good for ab, see if it is good for cd
n4 = c * n3;
n5 = n3 * d;
if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
{
E[0] = n3[0];
E[1] = n3[1];
E[2] = n3[2];
}
else
return MB_FAILURE;
}
else
return MB_FAILURE;
}
return MB_SUCCESS;
}
| double moab::IntxUtils::oriented_spherical_angle | ( | double * | A, |
| double * | B, | ||
| double * | C | ||
| ) | [static] |
Definition at line 814 of file IntxUtils.cpp.
References moab::angle().
Referenced by moab::IntxAreaUtils::area_spherical_polygon_girard(), and enforce_convexity().
{
// assume the same radius, sphere at origin
CartVect a( A ), b( B ), c( C );
CartVect normalOAB = a * b;
CartVect normalOCB = c * b;
CartVect orient = ( c - b ) * ( a - b );
double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
if( ang != ang )
{
// signal of a nan
std::cout << a << " " << b << " " << c << "\n";
std::cout << ang << "\n";
}
if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
return ang;
}
| ErrorCode moab::IntxUtils::remove_duplicate_vertices | ( | Interface * | mb, |
| EntityHandle | file_set, | ||
| double | merge_tol, | ||
| std::vector< Tag > & | tagList | ||
| ) | [static] |
Definition at line 1877 of file IntxUtils.cpp.
References ErrorCode, moab::Interface::get_entities_by_dimension(), MB_CHK_ERR, MB_SUCCESS, moab::MergeMesh::merge_all(), moab::Interface::remove_entities(), and remove_padded_vertices().
{
Range verts;
ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
MergeMesh mm( mb );
// remove the vertices from the set, before merging
rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
// now correct vertices that are repeated in polygons
rval = remove_padded_vertices( mb, file_set, tagList );
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::remove_padded_vertices | ( | Interface * | mb, |
| EntityHandle | file_set, | ||
| std::vector< Tag > & | tagList | ||
| ) | [static] |
Definition at line 1897 of file IntxUtils.cpp.
References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_element(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), moab::Range::insert(), MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, MBQUAD, MBTRI, moab::Interface::remove_entities(), moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().
Referenced by moab::NCHelperDomain::create_mesh(), moab::NCHelperScrip::create_mesh(), and remove_duplicate_vertices().
{
// now correct vertices that are repeated in polygons
Range cells;
ErrorCode rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
Range verts;
rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
Range modifiedCells; // will be deleted at the end; keep the gid
Range newCells;
for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
{
EntityHandle cell = *cit;
const EntityHandle* connec = NULL;
int num_verts = 0;
rval = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
std::vector< EntityHandle > newConnec;
newConnec.push_back( connec[0] ); // at least one vertex
int index = 0;
int new_size = 1;
while( index < num_verts - 2 )
{
int next_index = ( index + 1 );
if( connec[next_index] != newConnec[new_size - 1] )
{
newConnec.push_back( connec[next_index] );
new_size++;
}
index++;
}
// add the last one only if different from previous and first node
if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
{
newConnec.push_back( connec[num_verts - 1] );
new_size++;
}
if( new_size < num_verts )
{
// cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
modifiedCells.insert( cell );
// create a new cell with type triangle, quad or polygon
EntityType type = MBTRI;
if( new_size == 3 )
type = MBTRI;
else if( new_size == 4 )
type = MBQUAD;
else if( new_size > 4 )
type = MBPOLYGON;
// create new cell
EntityHandle newCell;
rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
// set the old id to the new element
newCells.insert( newCell );
double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
for( size_t i = 0; i < tagList.size(); i++ )
{
rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
}
}
}
rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
return MB_SUCCESS;
}
| ErrorCode moab::IntxUtils::reverse_gnomonic_projection | ( | const double & | c1, |
| const double & | c2, | ||
| double | R, | ||
| int | plane, | ||
| CartVect & | pos | ||
| ) | [static] |
Definition at line 450 of file IntxUtils.cpp.
References MB_SUCCESS, and moab::R.
Referenced by moab::IntxRllCssphere::findNodes(), moab::Intx2MeshOnSphere::findNodes(), get_barycenters(), and moab::BoundBox::update_box_spherical_elem().
{
// the new point will be on line beta*pos
double len = sqrt( c1 * c1 + c2 * c2 + R * R );
double beta = R / len; // it is less than 1, in general
switch( plane )
{
case 1: {
// the plane with x = R; x>y, x>z
// c1->y, c2->z
pos[0] = beta * R;
pos[1] = c1 * beta;
pos[2] = c2 * beta;
break;
}
case 2: {
// y = R -> zx
pos[1] = R * beta;
pos[2] = c1 * beta;
pos[0] = c2 * beta;
break;
}
case 3: {
// x=-R, -> yz
pos[0] = -R * beta;
pos[1] = -c1 * beta; // the sign is to preserve orientation
pos[2] = c2 * beta;
break;
}
case 4: {
// y = -R
pos[1] = -R * beta;
pos[2] = -c1 * beta; // the sign is to preserve orientation
pos[0] = c2 * beta;
break;
}
case 5: {
// z = -R
pos[2] = -R * beta;
pos[0] = -c1 * beta; // the sign is to preserve orientation
pos[1] = c2 * beta;
break;
}
case 6: {
pos[2] = R * beta;
pos[0] = c1 * beta;
pos[1] = c2 * beta;
break;
}
default:
return MB_FAILURE; // error
}
return MB_SUCCESS; // no error
}
| ErrorCode moab::IntxUtils::ScaleToRadius | ( | Interface * | mb, |
| EntityHandle | set, | ||
| double | R | ||
| ) | [static] |
Definition at line 771 of file IntxUtils.cpp.
References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::CartVect::length(), MB_SUCCESS, MBVERTEX, and moab::Interface::set_coords().
Referenced by CreateTempestMesh(), global_gnomonic_projection(), and main().
{
Range nodes;
ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
if( rval != moab::MB_SUCCESS ) return rval;
// one by one, get the node and project it on the sphere, with a radius given
// the center of the sphere is at 0,0,0
for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
{
EntityHandle nd = *nit;
CartVect pos;
rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
if( rval != moab::MB_SUCCESS ) return rval;
double len = pos.length();
if( len == 0. ) return MB_FAILURE;
pos = R / len * pos;
rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
if( rval != moab::MB_SUCCESS ) return rval;
}
return MB_SUCCESS;
}
| int moab::IntxUtils::SortAndRemoveDoubles2 | ( | double * | P, |
| int & | nP, | ||
| double | epsilon | ||
| ) | [static] |
Definition at line 98 of file IntxUtils.cpp.
References moab::angleAndIndex::angle, moab::angleCompare(), dist2(), moab::angleAndIndex::index, and MAXEDGES.
Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().
{
if( nP < 2 ) return 0; // nothing to do
// center of gravity for the points
double c[2] = { 0., 0. };
int k = 0;
for( k = 0; k < nP; k++ )
{
c[0] += P[2 * k];
c[1] += P[2 * k + 1];
}
c[0] /= nP;
c[1] /= nP;
// how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
// intersection points
struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
for( k = 0; k < nP; k++ )
{
double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
if( x != 0. || y != 0. )
{
pairAngleIndex[k].angle = atan2( y, x );
}
else
{
pairAngleIndex[k].angle = 0;
// it would mean that the cells are touching at a vertex
}
pairAngleIndex[k].index = k;
}
// this should be faster than the bubble sort we had before
std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
// copy now to a new double array
double PCopy[10 * MAXEDGES]; // the same dimension as P; very conservative, but faster than
// reallocate for a vector
for( k = 0; k < nP; k++ ) // index will show where it should go now;
{
int ck = pairAngleIndex[k].index;
PCopy[2 * k] = P[2 * ck];
PCopy[2 * k + 1] = P[2 * ck + 1];
}
// now copy from PCopy over original P location
std::copy( PCopy, PCopy + 2 * nP, P );
// eliminate duplicates, finally
int i = 0, j = 1; // the next one; j may advance faster than i
// check the unit
// double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
// than 1.e-5 cm
while( j < nP )
{
double d2 = dist2( &P[2 * i], &P[2 * j] );
if( d2 > epsilon_1 )
{
i++;
P[2 * i] = P[2 * j];
P[2 * i + 1] = P[2 * j + 1];
}
j++;
}
// test also the last point with the first one (index 0)
// the first one could be at -PI; last one could be at +PI, according to atan2 span
double d2 = dist2( P, &P[2 * i] ); // check the first and last points (ordered from -pi to +pi)
if( d2 > epsilon_1 )
{
nP = i + 1;
}
else
nP = i; // effectively delete the last point (that would have been the same with first)
if( nP == 0 ) nP = 1; // we should be left with at least one point we already tested if nP is 0 originally
return 0;
}
| CartVect moab::IntxUtils::spherical_to_cart | ( | IntxUtils::SphereCoords & | sc | ) | [static] |
Definition at line 762 of file IntxUtils.cpp.
References moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.
Referenced by add_field_value(), IntxUtilsCSLAM::departure_point_case1(), departure_point_swirl(), departure_point_swirl_rot(), main(), set_density(), and IntxUtilsCSLAM::smooth_field().
{
CartVect res;
res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
res[2] = sc.R * sin( sc.lat ); // z
return res;
}
| void moab::IntxUtils::transform_coordinates | ( | double * | avg_position, |
| int | projection_type | ||
| ) | [static] |
Definition at line 671 of file IntxUtils.cpp.
References decide_gnomonic_plane(), gnomonic_projection(), gnomonic_unroll(), and moab::R.
{
if( projection_type == 1 )
{
double R =
avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
R = sqrt( R );
double lat = asin( avg_position[2] / R );
double lon = atan2( avg_position[1], avg_position[0] );
avg_position[0] = lon;
avg_position[1] = lat;
avg_position[2] = R;
}
else if( projection_type == 2 ) // gnomonic projection
{
CartVect pos( avg_position );
int gplane;
IntxUtils::decide_gnomonic_plane( pos, gplane );
IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
avg_position[2] = 0;
IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
}
}
1.7.6.1