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