IntxUtils.hpp

Go to the documentation of this file.

/*
 * IntxUtils.hpp
 *
 * Created on: Oct 3, 2012
 */
 
#ifndef INTXUTILS_HPP_
#define INTXUTILS_HPP_
 
#include "moab/CartVect.hpp"
#include "moab/Core.hpp"
 
namespace moab
{
 
class IntxUtils
{
 public:
 // vec utilities that could be common between quads on a plane or sphere
 static inline double dist2( double* a, double* b )
 {
 double abx = b[0] - a[0], aby = b[1] - a[1];
 return sqrt( abx * abx + aby * aby );
 }
 
 static inline double area2D( double* a, double* b, double* c )
 {
 // (b-a)x(c-a) / 2
 return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
 }
 
 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 );
 // the marks will show what edges of blue intersect the red
 
 static ErrorCode EdgeIntersections2( double* blue,
 int nsBlue,
 double* red,
 int nsRed,
 int* markb,
 int* markr,
 double* points,
 int& nPoints );
 
 // special one, for intersection between rll (constant latitude) and cs quads
 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 );
 
 // vec utils related to gnomonic projection on a sphere
 // vec utils
 /*
 *
 * position on a sphere of radius R
 * if plane specified, use it; if not, return the plane, and the point in the plane
 * there are 6 planes, numbered 1 to 6
 * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
 *
 * projection on the plane will preserve the orientation, such that a triangle, quad pointing
 * outside the sphere will have a positive orientation in the projection plane
 * this is similar logic to
 * /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90 method:
 * function cart2face(cart3D) result(face_no)
 */
 static void decide_gnomonic_plane( const CartVect& pos, int& oPlane );
 
 // point on a sphere is projected on one of six planes, decided earlier
 
 static ErrorCode gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 );
 
 // given the position on plane (one out of 6), find out the position on sphere
 static ErrorCode reverse_gnomonic_projection( const double& c1,
 const double& c2,
 double R,
 int plane,
 CartVect& pos );
 
 // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
 static void gnomonic_unroll( double& c1, double& c2, double R, int plane );
 
 // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
 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 );
 /*
 * other methods to convert from spherical coord to cartesian, and back
 * A spherical coordinate triple is (R, lon, lat)
 * should we store it as a CartVect? probably not ...
 * /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90
 *
 enforce three facts:
 ! ==========================================================
 ! enforce three facts:
 !
 ! 1) lon at poles is defined to be zero
 !
 ! 2) Grid points must be separated by about .01 Meter (on earth)
 ! from pole to be considered "not the pole".
 !
 ! 3) range of lon is { 0<= lon < 2*pi }
 !
 ! 4) range of lat is from -pi/2 to pi/2; -pi/2 or pi/2 are the poles, so there the lon is 0
 !
 ! ==========================================================
 */
 struct SphereCoords
 {
 double R, lon, lat;
 };
 
 static SphereCoords cart_to_spherical( CartVect& );
 
 static CartVect spherical_to_cart( SphereCoords& );
 
 /*
 * create a mesh used mainly for visualization for now, with nodes corresponding to
 * GL points, a so-called refined mesh, with NP nodes in each direction.
 * input: a range of quads (coarse), and a desired order (NP is the number of points), so it
 * is order + 1
 *
 * output: a set with refined elements; with proper input, it should be pretty
 * similar to a Homme mesh read with ReadNC
 */
 // ErrorCode SpectralVisuMesh(Interface * mb, Range & input, int NP, EntityHandle & outputSet,
 // double tolerance);
 
 /*
 * given an entity set, get all nodes and project them on a sphere with given radius
 */
 static ErrorCode ScaleToRadius( Interface* mb, EntityHandle set, double R );
 
 static double distance_on_great_circle( CartVect& p1, CartVect& p2 );
 
 // break the nonconvex quads into triangles; remove the quad from the set? yes.
 // maybe radius is not needed;
 //
 static ErrorCode enforce_convexity( Interface* mb, EntityHandle set, int rank = 0 );
 
 // this could be larger than PI, because of orientation; useful for non-convex polygons
 static double oriented_spherical_angle( const double* A, const double* B, const double* C );
 
 // looking at quad connectivity, collapse to triangle if 2 nodes equal
 // then delete the old quad
 static ErrorCode fix_degenerate_quads( Interface* mb, EntityHandle set );
 
 // distance along a great circle on a sphere of radius 1
 static double distance_on_sphere( double la1, double te1, double la2, double te2 );
 
 /* End Analytical functions */
 
 /*
 * given 2 arcs AB and CD, compute the unique intersection point, if it exists
 * in between
 */
 static ErrorCode intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E );
 /*
 * given 2 arcs AB and CD, compute the intersection points, if it exists
 * AB is a great circle arc
 * CD is a constant latitude arc
 */
 static ErrorCode intersect_great_circle_arc_with_clat_arc( double* A,
 double* B,
 double* C,
 double* D,
 double R,
 double* E,
 int& np );
 
 // ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1);
 
 static int borderPointsOfCSinRLL( CartVect* redc,
 double* red2dc,
 int nsRed,
 CartVect* bluec,
 int nsBlue,
 int* blueEdgeType,
 double* P,
 int* side,
 double epsil );
 
 // used only by homme
 static ErrorCode deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set );
 
 // used to 'repair' scrip-like meshes
 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 );
};
 
class IntxAreaUtils
{
 public:
 enum AreaMethod
 {
 lHuiller = 0,
 Girard = 1,
 GaussQuadrature = 2
 };
 
 IntxAreaUtils( AreaMethod p_eAreaMethod = lHuiller ) : m_eAreaMethod( p_eAreaMethod ) {}
 
 ~IntxAreaUtils() {}
 
 /*
 * utilities to compute area of a polygon on which all edges are arcs of great circles on a
 * sphere
 */
 /*
 * this will compute the spherical angle ABC, when A, B, C are on a sphere of radius R
 * the radius will not be needed, usually, just for verification the points are indeed on that
 * sphere the center of the sphere is at origin (0,0,0) this angle can be used in Girard's
 * theorem to compute the area of a spherical polygon
 */
 double spherical_angle( const double* A, const double* B, const double* C, double Radius );
 
 double area_spherical_triangle( const double* A, const double* B, const double* C, double Radius );
 
 double area_spherical_polygon( const double* A, int N, double Radius, int* sign = NULL );
 
 double area_spherical_element( Interface* mb, EntityHandle elem, double R );
 
 double area_on_sphere( Interface* mb, EntityHandle set, double R );
 
 ErrorCode positive_orientation( Interface* mb, EntityHandle set, double R );
 
 private:
 /* lHuiller method for computing area on a spherical triangle */
 double area_spherical_triangle_lHuiller( const double* ptA, const double* ptB, const double* ptC, double Radius );
 
 /* lHuiller method for computing area on a spherical polygon */
 double area_spherical_polygon_lHuiller( const double* A, int N, double Radius, int* sign = NULL );
 
 /* Girard method for computing area on a spherical triangle with spherical excess */
 double area_spherical_triangle_girard( const double* A, const double* B, const double* C, double Radius );
 
 /* Girard method for computing area on a spherical polygon with spherical excess */
 double area_spherical_polygon_girard( const double* A, int N, double Radius );
 
#ifdef MOAB_HAVE_TEMPESTREMAP
 /* Gauss-quadrature based integration method for computing area on a spherical triangle */
 double area_spherical_triangle_GQ( const double* ptA, const double* ptB, const double* ptC );
 
 /* Gauss-quadrature based integration method for computing area on a spherical polygon element
 */
 double area_spherical_polygon_GQ( const double* A, int N );
#endif
 
 private:
 /* data */
 AreaMethod m_eAreaMethod;
};
 
} // namespace moab
#endif /* INTXUTILS_HPP_ */