 |
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 * Intx2Mesh.hpp
3 *
4 * Created on: Oct 2, 2012
5 *
6 */
7
8 #ifndef INTX2MESH_HPP_
9 #define INTX2MESH_HPP_
10
11 #include <iostream>
12 #include <sstream>
13 #include <fstream>
14 #include <map>
15 #include <ctime>
16 #include <cstdlib>
17 #include <cstdio>
18 #include <cstring>
19 #include "moab/Core.hpp"
20 #include "moab/Interface.hpp"
21 #include "moab/Range.hpp"
22 #include "moab/CartVect.hpp"
23 #include "moab/IntxMesh/IntxUtils.hpp"
24
25 // #define ENABLE_DEBUG
26
27 // maximum number of edges on each convex polygon of interest
28 #define MAXEDGES 10
29 #define MAXEDGES2 20 // used for coordinates in plane
30 #define CORRTAGNAME "__correspondent"
31
32 namespace moab
33 {
34
35 #define ERRORR( rval, str ) \
36 if( MB_SUCCESS != ( rval ) ) \
37 { \
38 std::cout << ( str ) << "\n"; \
39 return rval; \
40 }
41
42 #define ERRORV( rval, str ) \
43 if( MB_SUCCESS != ( rval ) ) \
44 { \
45 std::cout << ( str ) << "\n"; \
46 return; \
47 }
48
49 #ifdef MOAB_HAVE_MPI
50 // forward declarations
51 class ParallelComm;
52 class TupleList;
53 #endif
54
55 class Intx2Mesh
56 {
57 public:
58 Intx2Mesh( Interface* mbimpl );
59
60 virtual ~Intx2Mesh();
61
62 /*
63 * computes intersection between 2 sets;
64 * set 2 (mbs2) should be completely covered by set mbs1, as all elements
65 * from set 2 will be decomposed in polygons ; an area verification is done, and
66 * will signal elements from set 2 that are not "recovered"
67 *
68 * the resulting polygons are inserted in output set; total area from output set should be
69 * equal to total area of elements in set 2 (check that!)
70 */
71 ErrorCode intersect_meshes( EntityHandle mbs1, EntityHandle mbs2, EntityHandle& outputSet );
72
73 /*
74 * slower intx, use kd tree only, no adjacency, no adv front
75 */
76 ErrorCode intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet );
77
78 // mark could be (3 or 4, depending on type: ) no, it could go to 10
79 // no, it will be MAXEDGES = 10
80 // this is pure abstract, this needs to be implemented by
81 // all derivations
82 // the max number of intersection points could be 2*MAXEDGES
83 // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES)
84 // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon
85 // with 2*MAXEDGES, at most
86 // will also return the number of nodes of tgt and src elements
87 virtual ErrorCode computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
88 EntityHandle src,
89 double* P,
90 int& nP,
91 double& area,
92 int markb[MAXEDGES],
93 int markr[MAXEDGES],
94 int& nsidesSrc,
95 int& nsidesTgt,
96 bool check_boxes_first = false ) = 0;
97
98 // this is also abstract
99 virtual ErrorCode findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP ) = 0;
100
101 // this is also computing the area of the tgt cell in plane (gnomonic plane for sphere)
102 // setting the local variables:
103 // this will be done once per tgt cell, and once per localQueue for src cells
104 // const EntityHandle * tgtConn;
105 // CartVect redCoords[MAXEDGES];
106 // double redCoords2D[MAXEDGES2]; // these are in plane
107
108 virtual double setup_tgt_cell( EntityHandle tgt, int& nsTgt ) = 0;
109
110 virtual ErrorCode FindMaxEdgesInSet( EntityHandle eset, int& max_edges );
111 virtual ErrorCode FindMaxEdges( EntityHandle set1,
112 EntityHandle set2 ); // this needs to be called before any
113 // covering communication in parallel
114
115 virtual ErrorCode createTags();
116
117 virtual ErrorCode filterByMask(Range & cells);
118
119 ErrorCode DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag );
120
121 void set_error_tolerance( double eps )
122 {
123 epsilon_1 = eps;
124 epsilon_area = eps * sqrt( eps );
125 }
126
127 #ifdef MOAB_HAVE_MPI
128 void set_parallel_comm( moab::ParallelComm* pcomm )
129 {
130 parcomm = pcomm;
131 }
132 #endif
133 // void SetEntityType (EntityType tp) { type=tp;}
134
135 // clean some memory allocated
136 void clean();
137
138 void set_box_error( double berror )
139 {
140 box_error = berror;
141 }
142
143 ErrorCode create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set );
144
145 // in this method, used in parallel, each departure elements are already created, and at their
146 // positions the covering_set is output, will contain the departure cells that cover the euler
147 // set; some of these departure cells might come from different processors so the covering_set
148 // contains some elements from lagr_set and some elements that come from other procs we need to
149 // keep track of what processors "sent" the elements so we know were to send back the info about
150 // the tracers masses
151
152 ErrorCode create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set );
153
154 /* in this method, used in parallel, each cell from first mesh need to be sent to the
155 * tasks whose second mesh bounding boxes intersects bounding box of the cell
156 * this method assumes that the bounding boxes for the second mesh were computed already in a
157 * previous method (called build_processor_euler_boxes)
158 *
159 * the covering_set is output, will cover the second mesh (set) from the task;
160 * Some of the cells in the first mesh will be sent to multiple processors; we keep track of
161 * them using the global id of the vertices and global ids of the cells (we do not want to
162 * create multiple vertices with the same ids). The global id of the cells are needed just for
163 * debugging, the cells cannot come from 2 different tasks, but the vertices can
164 *
165 * Right now, we use crystal router, but an improvement might be to use direct communication
166 * (send_entities) on parallel comm
167 *
168 * param initial_distributed_set (IN) : the initial distribution of the first mesh (set)
169 *
170 * param (OUT) : the covering set in first mesh , which completely covers the second mesh set
171 */
172 #ifdef MOAB_HAVE_MPI
173
174 virtual ErrorCode build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic = true );
175 #endif
176 void correct_polygon( EntityHandle* foundIds, int& nP );
177 #ifdef MOAB_HAVE_MPI
178 // share vertices between the intersection target domains
179 ErrorCode resolve_intersection_sharing();
180 #endif
181 #ifdef ENABLE_DEBUG
182 void enable_debug()
183 {
184 dbg_1 = 1;
185 }
186 void disable_debug()
187 {
188 dbg_1 = 0;
189 }
190 #endif
191
192 #ifdef MOAB_HAVE_TEMPESTREMAP
193 friend class TempestRemapper;
194 #endif
195
196 protected: // so it can be accessed in derived classes, InPlane and OnSphere
197 Interface* mb;
198
199 EntityHandle mbs1;
200 EntityHandle mbs2;
201 Range rs1; // range set 1 (departure set, lagrange set, src set, manufactured set, target mesh)
202 Range rs2; // range set 2 (arrival set, euler set, tgt set, initial set, source mesh)
203
204 EntityHandle outSet; // will contain intersection
205 Tag gid; // global id tag will be used to set the parents of the intersection cell
206
207 // tags used in computation, advancing front
208 Tag TgtFlagTag; // to mark tgt quads already considered
209
210 Range TgtEdges; //
211
212 // tgt parent and src parent tags
213 // these will be on the out mesh
214 Tag tgtParentTag;
215 Tag srcParentTag;
216 Tag countTag;
217
218 Tag srcNeighTag; // will store neighbors for navigating easily in advancing front, for first
219 // mesh (src, target, lagrange)
220 Tag tgtNeighTag; // will store neighbors for navigating easily in advancing front, for second
221 // mesh (tgt, source, euler)
222
223 Tag neighTgtEdgeTag; // will store edge borders for each tgt cell
224
225 Tag orgSendProcTag; /// for coverage mesh, will store the original sender
226 Tag imaskTag; // if it exists, use it for filtering the source or target cells
227
228 // EntityType type; // this will be tri, quad or MBPOLYGON...
229
230 const EntityHandle* tgtConn;
231 const EntityHandle* srcConn;
232 CartVect tgtCoords[MAXEDGES];
233 CartVect srcCoords[MAXEDGES];
234 double tgtCoords2D[MAXEDGES2]; // these are in plane
235 double srcCoords2D[MAXEDGES2]; // these are in plane
236
237 #ifdef ENABLE_DEBUG
238 static int dbg_1;
239 std::ofstream mout_1[6]; // some debug files
240 #endif
241 // for each tgt edge, we keep a vector of extra nodes, coming from intersections
242 // use the index in TgtEdges range
243 // so the extra nodes on each tgt edge are kept track of
244 // only entity handles are in the vector, not the actual coordinates;
245 // actual coordinates are retrieved every time, which could be expensive
246 // maybe we should store the coordinates too, along with entity handles (more memory used,
247 // faster to retrieve)
248 std::vector< std::vector< EntityHandle >* > extraNodesVec;
249
250 double epsilon_1;
251 double epsilon_area;
252
253 std::vector< double > allBoxes;
254 double box_error;
255 /* \brief Local root of the kdtree */
256 EntityHandle localRoot;
257 Range localEnts; // this range is for local elements of interest, euler cells, or "first mesh"
258 unsigned int my_rank;
259
260 #ifdef MOAB_HAVE_MPI
261 ParallelComm* parcomm;
262 TupleList* remote_cells; // not used anymore for communication, just a container
263 TupleList* remote_cells_with_tracers; // these will be used now to update tracers on remote procs
264 std::map< int, EntityHandle > globalID_to_eh; // needed for parallel, mostly
265 #endif
266 int max_edges_1; // maximum number of edges in the lagrange set (first set, src)
267 int max_edges_2; // maximum number of edges in the euler set (second set, tgt)
268 int counting;
269 };
270
271 } /* namespace moab */
272 #endif /* INTX2MESH_HPP_ */
1.9.1.