MOAB: Mesh Oriented datABase  (version 5.3.0)
IntxUtils.hpp
Go to the documentation of this file.
00001 /*
00002  * IntxUtils.hpp
00003  *
00004  *  Created on: Oct 3, 2012
00005  */
00006 
00007 #ifndef INTXUTILS_HPP_
00008 #define INTXUTILS_HPP_
00009 
00010 #include "moab/CartVect.hpp"
00011 #include "moab/Core.hpp"
00012 
00013 namespace moab
00014 {
00015 
00016 class IntxUtils
00017 {
00018   public:
00019     // vec utilities that could be common between quads on a plane or sphere
00020     static inline double dist2( double* a, double* b )
00021     {
00022         double abx = b[0] - a[0], aby = b[1] - a[1];
00023         return sqrt( abx * abx + aby * aby );
00024     }
00025 
00026     static inline double area2D( double* a, double* b, double* c )
00027     {
00028         // (b-a)x(c-a) / 2
00029         return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
00030     }
00031 
00032     static int borderPointsOfXinY2( double* X, int nX, double* Y, int nY, double* P, int* side, double epsilon_area );
00033 
00034     static int SortAndRemoveDoubles2( double* P, int& nP, double epsilon );
00035     // the marks will show what edges of blue intersect the red
00036 
00037     static ErrorCode EdgeIntersections2( double* blue, int nsBlue, double* red, int nsRed, int* markb, int* markr,
00038                                          double* points, int& nPoints );
00039 
00040     // special one, for intersection between rll (constant latitude)  and cs quads
00041     static ErrorCode EdgeIntxRllCs( double* blue, CartVect* bluec, int* blueEdgeType, int nsBlue, double* red,
00042                                     CartVect* redc, int nsRed, int* markb, int* markr, int plane, double Radius,
00043                                     double* points, int& nPoints );
00044 
00045     // vec utils related to gnomonic projection on a sphere
00046     // vec utils
00047     /*
00048      *
00049      * position on a sphere of radius R
00050      * if plane specified, use it; if not, return the plane, and the point in the plane
00051      * there are 6 planes, numbered 1 to 6
00052      * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
00053      *
00054      * projection on the plane will preserve the orientation, such that a triangle, quad pointing
00055      * outside the sphere will have a positive orientation in the projection plane
00056      * this is similar logic to
00057      * /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90 method:
00058      * function cart2face(cart3D) result(face_no)
00059      */
00060     static void decide_gnomonic_plane( const CartVect& pos, int& oPlane );
00061 
00062     // point on a sphere is projected on one of six planes, decided earlier
00063 
00064     static ErrorCode gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 );
00065 
00066     // given the position on plane (one out of 6), find out the position on sphere
00067     static ErrorCode reverse_gnomonic_projection( const double& c1, const double& c2, double R, int plane,
00068                                                   CartVect& pos );
00069 
00070     // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
00071     static void gnomonic_unroll( double& c1, double& c2, double R, int plane );
00072 
00073     // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
00074     static ErrorCode global_gnomonic_projection( Interface* mb, EntityHandle inSet, double R, bool centers_only,
00075                                                  EntityHandle& outSet );
00076 
00077     static void transform_coordinates( double* avg_position, int projection_type );
00078     /*
00079     *   other methods to convert from spherical coord to cartesian, and back
00080     *   A spherical coordinate triple is (R, lon, lat)
00081     *   should we store it as a CartVect? probably not ...
00082     *   /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90
00083     *
00084         enforce three facts:
00085         ! ==========================================================
00086         ! enforce three facts:
00087         !
00088         ! 1) lon at poles is defined to be zero
00089         !
00090         ! 2) Grid points must be separated by about .01 Meter (on earth)
00091         !    from pole to be considered "not the pole".
00092         !
00093         ! 3) range of lon is { 0<= lon < 2*pi }
00094         !
00095         ! 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
00096         !
00097         ! ==========================================================
00098     */
00099     struct SphereCoords
00100     {
00101         double R, lon, lat;
00102     };
00103 
00104     static SphereCoords cart_to_spherical( CartVect& );
00105 
00106     static CartVect spherical_to_cart( SphereCoords& );
00107 
00108     /*
00109      *  create a mesh used mainly for visualization for now, with nodes corresponding to
00110      *   GL points, a so-called refined mesh, with NP nodes in each direction.
00111      *   input: a range of quads (coarse), and  a desired order (NP is the number of points), so it
00112      *   is order + 1
00113      *
00114      *   output: a set with refined elements; with proper input, it should be pretty
00115      *   similar to a Homme mesh read with ReadNC
00116      */
00117     // ErrorCode SpectralVisuMesh(Interface * mb, Range & input, int NP, EntityHandle & outputSet,
00118     // double tolerance);
00119 
00120     /*
00121      * given an entity set, get all nodes and project them on a sphere with given radius
00122      */
00123     static ErrorCode ScaleToRadius( Interface* mb, EntityHandle set, double R );
00124 
00125     static double distance_on_great_circle( CartVect& p1, CartVect& p2 );
00126 
00127     // break the nonconvex quads into triangles; remove the quad from the set? yes.
00128     // maybe radius is not needed;
00129     //
00130     static ErrorCode enforce_convexity( Interface* mb, EntityHandle set, int rank = 0 );
00131 
00132     // this could be larger than PI, because of orientation; useful for non-convex polygons
00133     static double oriented_spherical_angle( double* A, double* B, double* C );
00134 
00135     // looking at quad connectivity, collapse to triangle if 2 nodes equal
00136     // then delete the old quad
00137     static ErrorCode fix_degenerate_quads( Interface* mb, EntityHandle set );
00138 
00139     // distance along a great circle on a sphere of radius 1
00140     static double distance_on_sphere( double la1, double te1, double la2, double te2 );
00141 
00142     /* End Analytical functions */
00143 
00144     /*
00145      * given 2 arcs AB and CD, compute the unique intersection point, if it exists
00146      *  in between
00147      */
00148     static ErrorCode intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E );
00149     /*
00150      * given 2 arcs AB and CD, compute the intersection points, if it exists
00151      *  AB is a great circle arc
00152      *  CD is a constant latitude arc
00153      */
00154     static ErrorCode intersect_great_circle_arc_with_clat_arc( double* A, double* B, double* C, double* D, double R,
00155                                                                double* E, int& np );
00156 
00157     // ErrorCode  set_edge_type_flag(Interface * mb, EntityHandle sf1);
00158 
00159     static int borderPointsOfCSinRLL( CartVect* redc, double* red2dc, int nsRed, CartVect* bluec, int nsBlue,
00160                                       int* blueEdgeType, double* P, int* side, double epsil );
00161 
00162     // used only by homme
00163     static ErrorCode deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set );
00164 
00165     // used to 'repair' scrip-like meshes
00166     static ErrorCode remove_duplicate_vertices( Interface* mb, EntityHandle file_set, double merge_tol,
00167                                                 std::vector< Tag >& tagList );
00168 };
00169 
00170 class IntxAreaUtils
00171 {
00172   public:
00173     enum AreaMethod
00174     {
00175         lHuiller        = 0,
00176         Girard          = 1,
00177         GaussQuadrature = 2
00178     };
00179 
00180     IntxAreaUtils( AreaMethod p_eAreaMethod = lHuiller ) : m_eAreaMethod( p_eAreaMethod ) {}
00181 
00182     ~IntxAreaUtils() {}
00183 
00184     /*
00185      * utilities to compute area of a polygon on which all edges are arcs of great circles on a
00186      * sphere
00187      */
00188     /*
00189      * this will compute the spherical angle ABC, when A, B, C are on a sphere of radius R
00190      *  the radius will not be needed, usually, just for verification the points are indeed on that
00191      * sphere the center of the sphere is at origin (0,0,0) this angle can be used in Girard's
00192      * theorem to compute the area of a spherical polygon
00193      */
00194     double spherical_angle( double* A, double* B, double* C, double Radius );
00195 
00196     double area_spherical_triangle( double* A, double* B, double* C, double Radius );
00197 
00198     double area_spherical_polygon( double* A, int N, double Radius, int* sign = NULL );
00199 
00200     double area_spherical_element( Interface* mb, EntityHandle elem, double R );
00201 
00202     double area_on_sphere( Interface* mb, EntityHandle set, double R );
00203 
00204     ErrorCode positive_orientation( Interface* mb, EntityHandle set, double R );
00205 
00206   private:
00207     /* lHuiller method for computing area on a spherical triangle */
00208     double area_spherical_triangle_lHuiller( double* ptA, double* ptB, double* ptC, double Radius );
00209 
00210     /* lHuiller method for computing area on a spherical polygon */
00211     double area_spherical_polygon_lHuiller( double* A, int N, double Radius, int* sign = NULL );
00212 
00213     /* Girard method for computing area on a spherical triangle with spherical excess */
00214     double area_spherical_triangle_girard( double* A, double* B, double* C, double Radius );
00215 
00216     /* Girard method for computing area on a spherical polygon with spherical excess */
00217     double area_spherical_polygon_girard( double* A, int N, double Radius );
00218 
00219 #ifdef MOAB_HAVE_TEMPESTREMAP
00220     /* Gauss-quadrature based integration method for computing area on a spherical triangle */
00221     double area_spherical_triangle_GQ( double* ptA, double* ptB, double* ptC, double Radius );
00222 
00223     /* Gauss-quadrature based integration method for computing area on a spherical polygon element
00224      */
00225     double area_spherical_polygon_GQ( double* A, int N, double Radius );
00226 #endif
00227 
00228   private:
00229     /* data */
00230     AreaMethod m_eAreaMethod;
00231 };
00232 
00233 }  // namespace moab
00234 #endif /* INTXUTILS_HPP_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines