Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
IntxUtils.hpp
Go to the documentation of this file.
1 /*
2  * IntxUtils.hpp
3  *
4  * Created on: Oct 3, 2012
5  */
6 
7 #ifndef INTXUTILS_HPP_
8 #define INTXUTILS_HPP_
9 
10 #include "moab/CartVect.hpp"
11 #include "moab/Core.hpp"
12 namespace moab
13 {
14 
15 class IntxUtils
16 {
17  public:
18  // vec utilities that could be common between quads on a plane or sphere
19  static inline double dist2( double* a, double* b )
20  {
21  double abx = b[0] - a[0], aby = b[1] - a[1];
22  return sqrt( abx * abx + aby * aby );
23  }
24 
25  static inline double area2D( double* a, double* b, double* c )
26  {
27  // (b-a)x(c-a) / 2
28  return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
29  }
30 
31  static int borderPointsOfXinY2( double* X, int nX, double* Y, int nY, double* P, int* side, double epsilon_area );
32 
33  static int SortAndRemoveDoubles2( double* P, int& nP, double epsilon );
34  // the marks will show what edges of blue intersect the red
35 
36  static ErrorCode EdgeIntersections2( double* blue,
37  int nsBlue,
38  double* red,
39  int nsRed,
40  int* markb,
41  int* markr,
42  double* points,
43  int& nPoints );
44 
45  // special one, for intersection between rll (constant latitude) and cs quads
46  static ErrorCode EdgeIntxRllCs( double* blue,
47  CartVect* bluec,
48  int* blueEdgeType,
49  int nsBlue,
50  double* red,
51  CartVect* redc,
52  int nsRed,
53  int* markb,
54  int* markr,
55  int plane,
56  double Radius,
57  double* points,
58  int& nPoints );
59 
60  // vec utils related to gnomonic projection on a sphere
61  // vec utils
62  /*
63  *
64  * position on a sphere of radius R
65  * if plane specified, use it; if not, return the plane, and the point in the plane
66  * there are 6 planes, numbered 1 to 6
67  * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
68  *
69  * projection on the plane will preserve the orientation, such that a triangle, quad pointing
70  * outside the sphere will have a positive orientation in the projection plane
71  * this is similar logic to
72  * /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90 method:
73  * function cart2face(cart3D) result(face_no)
74  */
75  static void decide_gnomonic_plane( const CartVect& pos, int& oPlane );
76 
77  // point on a sphere is projected on one of six planes, decided earlier
78 
79  static ErrorCode gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 );
80 
81  // point on a sphere is projected on plane decided by gnomonic center P, and u, v computed in a previous method
82  // these points are in a list axis[3], P, u and v
83 
84  static ErrorCode gnomonic_projection_generalized( const CartVect& pos,
85  const CartVect axis[3],
86  double& c1,
87  double& c2 );
88 
89  // given a mesh on a hemisphere, and a point P that defines the hemisphere, project the mesh
90  // on a plane tangent at P (gnomonic plane at P)
91  static ErrorCode global_gnomonic_projection_general( Interface* mb,
92  EntityHandle inSet,
93  CartVect P,
94  EntityHandle& outSet );
95  // given the position on plane (one out of 6), find out the position on sphere
96  static ErrorCode reverse_gnomonic_projection( const double& c1,
97  const double& c2,
98  double R,
99  int plane,
100  CartVect& pos );
101 
102  // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
103  static void gnomonic_unroll( double& c1, double& c2, double R, int plane );
104 
105  // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
106  static ErrorCode global_gnomonic_projection( Interface* mb,
107  EntityHandle inSet,
108  double R,
109  bool centers_only,
110  EntityHandle& outSet );
111 
112  // given a point P, |P| > 0, find u,v directions for the gnomonic projection plane perpendicular to OP
113  // the result will be a cartesian coordinate system centered at P, and 2 unit vectors at P that
114  // define the gnomonic plane orientation around P, such as a positive oriented triangle on a sphere will
115  // remain positive oriented after projection in the plane.
116  // all vertices in this projection will be on one hemisphere to the side of P
117  // the orientation of u and v uniquely determine the gnomonic plane orientation, and u x v should be in OP direction
118  // also, u and v are in plane, and perpendicular to each other
119  static ErrorCode gnomonic_projection_plane_at_point( CartVect P, CartVect& u, CartVect& v );
120 
121  static void transform_coordinates( double* avg_position, int projection_type );
122  /*
123  * other methods to convert from spherical coord to cartesian, and back
124  * A spherical coordinate triple is (R, lon, lat)
125  * should we store it as a CartVect? probably not ...
126  * /cesm1_0_4/models/atm/cam/src/dynamics/homme/share/coordinate_systems_mod.F90
127  *
128  enforce three facts:
129  ! ==========================================================
130  ! enforce three facts:
131  !
132  ! 1) lon at poles is defined to be zero
133  !
134  ! 2) Grid points must be separated by about .01 Meter (on earth)
135  ! from pole to be considered "not the pole".
136  !
137  ! 3) range of lon is { 0<= lon < 2*pi }
138  !
139  ! 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
140  !
141  ! ==========================================================
142  */
143  struct SphereCoords
144  {
145  double R, lon, lat;
146  };
147 
148  static SphereCoords cart_to_spherical( CartVect& );
149 
150  static CartVect spherical_to_cart( SphereCoords& );
151 
152  /*
153  * create a mesh used mainly for visualization for now, with nodes corresponding to
154  * GL points, a so-called refined mesh, with NP nodes in each direction.
155  * input: a range of quads (coarse), and a desired order (NP is the number of points), so it
156  * is order + 1
157  *
158  * output: a set with refined elements; with proper input, it should be pretty
159  * similar to a Homme mesh read with ReadNC
160  */
161  // ErrorCode SpectralVisuMesh(Interface * mb, Range & input, int NP, EntityHandle & outputSet,
162  // double tolerance);
163 
164  /*
165  * given an entity set, get all nodes and project them on a sphere with given radius
166  */
167  static ErrorCode ScaleToRadius( Interface* mb, EntityHandle set, double R );
168 
169  static double distance_on_great_circle( CartVect& p1, CartVect& p2 );
170 
171  // break the nonconvex quads into triangles; remove the quad from the set? yes.
172  // maybe radius is not needed;
173  //
174  static ErrorCode enforce_convexity( Interface* mb, EntityHandle set, int rank = 0 );
175 
176  // this could be larger than PI, because of orientation; useful for non-convex polygons
177  static double oriented_spherical_angle( const double* A, const double* B, const double* C );
178 
179  // looking at quad connectivity, collapse to triangle if 2 nodes equal
180  // then delete the old quad
181  static ErrorCode fix_degenerate_quads( Interface* mb, EntityHandle set );
182 
183  // distance along a great circle on a sphere of radius 1
184  static double distance_on_sphere( double la1, double te1, double la2, double te2 );
185 
186  /* End Analytical functions */
187 
188  /*
189  * given 2 arcs AB and CD, compute the unique intersection point, if it exists
190  * in between
191  */
192  static ErrorCode intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E );
193  /*
194  * given 2 arcs AB and CD, compute the intersection points, if it exists
195  * AB is a great circle arc
196  * CD is a constant latitude arc
197  */
198  static ErrorCode intersect_great_circle_arc_with_clat_arc( double* A,
199  double* B,
200  double* C,
201  double* D,
202  double R,
203  double* E,
204  int& np );
205 
206  // ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1);
207 
208  static int borderPointsOfCSinRLL( CartVect* redc,
209  double* red2dc,
210  int nsRed,
211  CartVect* bluec,
212  int nsBlue,
213  int* blueEdgeType,
214  double* P,
215  int* side,
216  double epsil );
217 
218  // used only by homme
219  static ErrorCode deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set );
220 
221  // used to 'repair' scrip-like meshes
222  static ErrorCode remove_duplicate_vertices( Interface* mb,
223  EntityHandle file_set,
224  double merge_tol,
225  std::vector< Tag >& tagList );
226 
227  static ErrorCode remove_padded_vertices( Interface* mb, EntityHandle file_set, std::vector< Tag >& tagList );
228  // used now to compute maximum diagonal for a range of cells
229  static ErrorCode max_diagonal( Interface* mb, Range cells, int max_edges, double& diagonal );
230 };
231 
232 class IntxAreaUtils
233 {
234  public:
235  enum AreaMethod
236  {
237  lHuiller = 0,
238  Girard = 1,
239  GaussQuadrature = 2
240  };
241 
242  IntxAreaUtils( AreaMethod p_eAreaMethod = lHuiller ) : m_eAreaMethod( p_eAreaMethod ) {}
243 
244  ~IntxAreaUtils() {}
245 
246  /*
247  * utilities to compute area of a polygon on which all edges are arcs of great circles on a
248  * sphere
249  */
250  /*
251  * this will compute the spherical angle ABC, when A, B, C are on a sphere of radius R
252  * the radius will not be needed, usually, just for verification the points are indeed on that
253  * sphere the center of the sphere is at origin (0,0,0) this angle can be used in Girard's
254  * theorem to compute the area of a spherical polygon
255  */
256  double spherical_angle( const double* A, const double* B, const double* C, double Radius );
257 
258  double area_spherical_triangle( const double* A, const double* B, const double* C, double Radius );
259 
260  double area_spherical_polygon( const double* A, int N, double Radius, int* sign = NULL );
261 
262  double area_spherical_element( Interface* mb, EntityHandle elem, double R );
263 
264  double area_on_sphere( Interface* mb, EntityHandle set, double R );
265 
266  ErrorCode positive_orientation( Interface* mb, EntityHandle set, double R );
267 
268  private:
269  /* lHuiller method for computing area on a spherical triangle */
270  double area_spherical_triangle_lHuiller( const double* ptA, const double* ptB, const double* ptC, double Radius );
271 
272  /* lHuiller method for computing area on a spherical polygon */
273  double area_spherical_polygon_lHuiller( const double* A, int N, double Radius, int* sign = NULL );
274 
275  /* Girard method for computing area on a spherical triangle with spherical excess */
276  double area_spherical_triangle_girard( const double* A, const double* B, const double* C, double Radius );
277 
278  /* Girard method for computing area on a spherical polygon with spherical excess */
279  double area_spherical_polygon_girard( const double* A, int N, double Radius );
280 
281 #ifdef MOAB_HAVE_TEMPESTREMAP
282  /* Gauss-quadrature based integration method for computing area on a spherical triangle */
283  double area_spherical_triangle_GQ( const double* ptA, const double* ptB, const double* ptC );
284 
285  /* Gauss-quadrature based integration method for computing area on a spherical polygon element
286  */
287  double area_spherical_polygon_GQ( const double* A, int N );
288 #endif
289 
290  private:
291  /* data */
292  AreaMethod m_eAreaMethod;
293 };
294 
295 } // namespace moab
296 #endif /* INTXUTILS_HPP_ */