Branch data Line data Source code
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 <time.h>
16 : : #include <stdlib.h>
17 : : #include <stdio.h>
18 : : #include <string.h>
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, EntityHandle src, double* P, int& nP,
88 : : double& area, int markb[MAXEDGES], int markr[MAXEDGES],
89 : : int& nsidesSrc, int& nsidesTgt,
90 : : bool check_boxes_first = false ) = 0;
91 : :
92 : : // this is also abstract
93 : : virtual ErrorCode findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP ) = 0;
94 : :
95 : : // this is also computing the area of the tgt cell in plane (gnomonic plane for sphere)
96 : : // setting the local variables:
97 : : // this will be done once per tgt cell, and once per localQueue for src cells
98 : : // const EntityHandle * tgtConn;
99 : : // CartVect redCoords[MAXEDGES];
100 : : // double redCoords2D[MAXEDGES2]; // these are in plane
101 : :
102 : : virtual double setup_tgt_cell( EntityHandle tgt, int& nsTgt ) = 0;
103 : :
104 : : virtual ErrorCode FindMaxEdgesInSet( EntityHandle eset, int& max_edges );
105 : : virtual ErrorCode FindMaxEdges( EntityHandle set1,
106 : : EntityHandle set2 ); // this needs to be called before any
107 : : // covering communication in parallel
108 : :
109 : : virtual ErrorCode createTags();
110 : :
111 : : ErrorCode DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag );
112 : :
113 : 3 : void set_error_tolerance( double eps )
114 : : {
115 : 3 : epsilon_1 = eps;
116 : 3 : epsilon_area = eps * eps / 2;
117 : 3 : }
118 : :
119 : : #ifdef MOAB_HAVE_MPI
120 : 1 : void set_parallel_comm( moab::ParallelComm* pcomm )
121 : : {
122 : 1 : parcomm = pcomm;
123 : 1 : }
124 : : #endif
125 : : // void SetEntityType (EntityType tp) { type=tp;}
126 : :
127 : : // clean some memory allocated
128 : : void clean();
129 : :
130 : 1 : void set_box_error( double berror )
131 : : {
132 : 1 : box_error = berror;
133 : 1 : }
134 : :
135 : : ErrorCode create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set );
136 : :
137 : : // in this method, used in parallel, each departure elements are already created, and at their
138 : : // positions the covering_set is output, will contain the departure cells that cover the euler
139 : : // set; some of these departure cells might come from different processors so the covering_set
140 : : // contains some elements from lagr_set and some elements that come from other procs we need to
141 : : // keep track of what processors "sent" the elements so we know were to send back the info about
142 : : // the tracers masses
143 : :
144 : : ErrorCode create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set );
145 : :
146 : : /* in this method, used in parallel, each cell from first mesh need to be sent to the
147 : : * tasks whose second mesh bounding boxes intersects bounding box of the cell
148 : : * this method assumes that the bounding boxes for the second mesh were computed already in a
149 : : * previous method (called build_processor_euler_boxes)
150 : : *
151 : : * the covering_set is output, will cover the second mesh (set) from the task;
152 : : * Some of the cells in the first mesh will be sent to multiple processors; we keep track of
153 : : * them using the global id of the vertices and global ids of the cells (we do not want to
154 : : * create multiple vertices with the same ids). The global id of the cells are needed just for
155 : : * debugging, the cells cannot come from 2 different tasks, but the vertices can
156 : : *
157 : : * Right now, we use crystal router, but an improvement might be to use direct communication
158 : : * (send_entities) on parallel comm
159 : : *
160 : : * param initial_distributed_set (IN) : the initial distribution of the first mesh (set)
161 : : *
162 : : * param (OUT) : the covering set in first mesh , which completely covers the second mesh set
163 : : */
164 : : #ifdef MOAB_HAVE_MPI
165 : :
166 : : virtual ErrorCode build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts );
167 : : #endif
168 : : void correct_polygon( EntityHandle* foundIds, int& nP );
169 : : #ifdef MOAB_HAVE_MPI
170 : : // share vertices between the intersection target domains
171 : : ErrorCode resolve_intersection_sharing();
172 : : #endif
173 : : #ifdef ENABLE_DEBUG
174 : : void enable_debug()
175 : : {
176 : : dbg_1 = 1;
177 : : }
178 : : void disable_debug()
179 : : {
180 : : dbg_1 = 0;
181 : : }
182 : : #endif
183 : :
184 : : #ifdef MOAB_HAVE_TEMPESTREMAP
185 : : friend class TempestRemapper;
186 : : #endif
187 : :
188 : : protected: // so it can be accessed in derived classes, InPlane and OnSphere
189 : : Interface* mb;
190 : :
191 : : EntityHandle mbs1;
192 : : EntityHandle mbs2;
193 : : Range rs1; // range set 1 (departure set, lagrange set, src set, manufactured set, target mesh)
194 : : Range rs2; // range set 2 (arrival set, euler set, tgt set, initial set, source mesh)
195 : :
196 : : EntityHandle outSet; // will contain intersection
197 : : Tag gid; // global id tag will be used to set the parents of the intersection cell
198 : :
199 : : // tags used in computation, advancing front
200 : : Tag TgtFlagTag; // to mark tgt quads already considered
201 : :
202 : : Range TgtEdges; //
203 : :
204 : : // tgt parent and src parent tags
205 : : // these will be on the out mesh
206 : : Tag tgtParentTag;
207 : : Tag srcParentTag;
208 : : Tag countTag;
209 : :
210 : : Tag srcNeighTag; // will store neighbors for navigating easily in advancing front, for first
211 : : // mesh (src, target, lagrange)
212 : : Tag tgtNeighTag; // will store neighbors for navigating easily in advancing front, for second
213 : : // mesh (tgt, source, euler)
214 : :
215 : : Tag neighTgtEdgeTag; // will store edge borders for each tgt cell
216 : :
217 : : Tag orgSendProcTag; /// for coverage mesh, will store the original sender
218 : :
219 : : // EntityType type; // this will be tri, quad or MBPOLYGON...
220 : :
221 : : const EntityHandle* tgtConn;
222 : : const EntityHandle* srcConn;
223 : : CartVect tgtCoords[MAXEDGES];
224 : : CartVect srcCoords[MAXEDGES];
225 : : double tgtCoords2D[MAXEDGES2]; // these are in plane
226 : : double srcCoords2D[MAXEDGES2]; // these are in plane
227 : :
228 : : #ifdef ENABLE_DEBUG
229 : : static int dbg_1;
230 : : std::ofstream mout_1[6]; // some debug files
231 : : #endif
232 : : // for each tgt edge, we keep a vector of extra nodes, coming from intersections
233 : : // use the index in TgtEdges range
234 : : // so the extra nodes on each tgt edge are kept track of
235 : : // only entity handles are in the vector, not the actual coordinates;
236 : : // actual coordinates are retrieved every time, which could be expensive
237 : : // maybe we should store the coordinates too, along with entity handles (more memory used,
238 : : // faster to retrieve)
239 : : std::vector< std::vector< EntityHandle >* > extraNodesVec;
240 : :
241 : : double epsilon_1;
242 : : double epsilon_area;
243 : :
244 : : std::vector< double > allBoxes;
245 : : double box_error;
246 : : /* \brief Local root of the kdtree */
247 : : EntityHandle localRoot;
248 : : Range localEnts; // this range is for local elements of interest, euler cells, or "first mesh"
249 : : unsigned int my_rank;
250 : :
251 : : #ifdef MOAB_HAVE_MPI
252 : : ParallelComm* parcomm;
253 : : TupleList* remote_cells; // not used anymore for communication, just a container
254 : : TupleList* remote_cells_with_tracers; // these will be used now to update tracers on remote procs
255 : : std::map< int, EntityHandle > globalID_to_eh; // needed for parallel, mostly
256 : : #endif
257 : : int max_edges_1; // maximum number of edges in the lagrange set (first set, src)
258 : : int max_edges_2; // maximum number of edges in the euler set (second set, tgt)
259 : : int counting;
260 : : };
261 : :
262 : : } /* namespace moab */
263 : : #endif /* INTX2MESH_HPP_ */
|
