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