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