 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
93 static ErrorCode global_gnomonic_projection( Interface* mb,
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 */
121 struct SphereCoords
122 {
123 double R, lon, lat;
124 };
125
126 static SphereCoords cart_to_spherical( CartVect& );
127
128 static CartVect spherical_to_cart( SphereCoords& );
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
159 static ErrorCode fix_degenerate_quads( Interface* mb, EntityHandle set );
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 */
176 static ErrorCode intersect_great_circle_arc_with_clat_arc( double* A,
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
200 static ErrorCode remove_duplicate_vertices( Interface* mb,
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 // used now to compute maximum diagonal for a range of cells
207 static ErrorCode max_diagonal(Interface* mb, Range cells, int max_edges, double & diagonal);
208 };
209
210 class IntxAreaUtils
211 {
212 public:
213 enum AreaMethod
214 {
215 lHuiller = 0,
216 Girard = 1,
217 GaussQuadrature = 2
218 };
219
220 IntxAreaUtils( AreaMethod p_eAreaMethod = lHuiller ) : m_eAreaMethod( p_eAreaMethod ) {}
221
222 ~IntxAreaUtils() {}
223
224 /*
225 * utilities to compute area of a polygon on which all edges are arcs of great circles on a
226 * sphere
227 */
228 /*
229 * this will compute the spherical angle ABC, when A, B, C are on a sphere of radius R
230 * the radius will not be needed, usually, just for verification the points are indeed on that
231 * sphere the center of the sphere is at origin (0,0,0) this angle can be used in Girard's
232 * theorem to compute the area of a spherical polygon
233 */
234 double spherical_angle( const double* A, const double* B, const double* C, double Radius );
235
236 double area_spherical_triangle( const double* A, const double* B, const double* C, double Radius );
237
238 double area_spherical_polygon( const double* A, int N, double Radius, int* sign = NULL );
239
240 double area_spherical_element( Interface* mb, EntityHandle elem, double R );
241
242 double area_on_sphere( Interface* mb, EntityHandle set, double R );
243
244 ErrorCode positive_orientation( Interface* mb, EntityHandle set, double R );
245
246 private:
247 /* lHuiller method for computing area on a spherical triangle */
248 double area_spherical_triangle_lHuiller( const double* ptA, const double* ptB, const double* ptC, double Radius );
249
250 /* lHuiller method for computing area on a spherical polygon */
251 double area_spherical_polygon_lHuiller( const double* A, int N, double Radius, int* sign = NULL );
252
253 /* Girard method for computing area on a spherical triangle with spherical excess */
254 double area_spherical_triangle_girard( const double* A, const double* B, const double* C, double Radius );
255
256 /* Girard method for computing area on a spherical polygon with spherical excess */
257 double area_spherical_polygon_girard( const double* A, int N, double Radius );
258
259 #ifdef MOAB_HAVE_TEMPESTREMAP
260 /* Gauss-quadrature based integration method for computing area on a spherical triangle */
261 double area_spherical_triangle_GQ( const double* ptA, const double* ptB, const double* ptC );
262
263 /* Gauss-quadrature based integration method for computing area on a spherical polygon element
264 */
265 double area_spherical_polygon_GQ( const double* A, int N );
266 #endif
267
268 private:
269 /* data */
270 AreaMethod m_eAreaMethod;
271 };
272
273 } // namespace moab
274 #endif /* INTXUTILS_HPP_ */
1.9.1.