Branch data Line data Source code
1 : : /*
2 : : * Intx2MeshOnSphere.cpp
3 : : *
4 : : * Created on: Oct 3, 2012
5 : : */
6 : :
7 : : #ifdef WIN32 /* windows */
8 : : #define _USE_MATH_DEFINES // For M_PI
9 : : #endif
10 : : #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
11 : : #include "moab/IntxMesh/IntxUtils.hpp"
12 : : #include "moab/GeomUtil.hpp"
13 : : #ifdef MOAB_HAVE_MPI
14 : : #include "moab/ParallelComm.hpp"
15 : : #endif
16 : : #include "MBTagConventions.hpp"
17 : :
18 : : // #define ENABLE_DEBUG
19 : : //#define CHECK_CONVEXITY
20 : : namespace moab
21 : : {
22 : :
23 : 1 : Intx2MeshOnSphere::Intx2MeshOnSphere( Interface* mbimpl, IntxAreaUtils::AreaMethod amethod )
24 : 1 : : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
25 : : {
26 : 1 : }
27 : :
28 [ - + ]: 2 : Intx2MeshOnSphere::~Intx2MeshOnSphere() {}
29 : :
30 : : /*
31 : : * return also the area for robustness verification
32 : : */
33 : 434 : double Intx2MeshOnSphere::setup_tgt_cell( EntityHandle tgt, int& nsTgt )
34 : : {
35 : :
36 : : // get coordinates of the target quad, to decide the gnomonic plane
37 : 434 : double cellArea = 0;
38 : :
39 : : int num_nodes;
40 [ + - ][ - + ]: 434 : ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
[ # # ][ # # ]
41 : :
42 : 434 : nsTgt = num_nodes;
43 : : // account for possible padded polygons
44 [ - + ][ # # ]: 434 : while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
45 : 0 : nsTgt--;
46 : :
47 : : // CartVect coords[4];
48 [ + - ][ + - ]: 434 : rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
[ - + ][ # # ]
[ # # ]
49 : :
50 : 434 : CartVect middle = tgtCoords[0];
51 [ + + ]: 1736 : for( int i = 1; i < nsTgt; i++ )
52 [ + - ]: 1302 : middle += tgtCoords[i];
53 [ + - ]: 434 : middle = 1. / nsTgt * middle;
54 : :
55 [ + - ]: 434 : IntxUtils::decide_gnomonic_plane( middle, plane ); // output the plane
56 [ + + ]: 2170 : for( int j = 0; j < nsTgt; j++ )
57 : : {
58 : : // populate coords in the plane for intersection
59 : : // they should be oriented correctly, positively
60 [ + - ][ - + ]: 1736 : rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
[ # # ][ # # ]
61 : : }
62 : :
63 [ + + ]: 1302 : for( int j = 1; j < nsTgt - 1; j++ )
64 [ + - ]: 868 : cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
65 : :
66 : : // take target coords in order and compute area in plane
67 : 434 : return cellArea;
68 : : }
69 : :
70 : : /* the elements are convex for sure, then do a gnomonic projection of both,
71 : : * compute intersection in the plane, then go back to the sphere for the points
72 : : * */
73 : 724 : ErrorCode Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt, EntityHandle src, double* P,
74 : : int& nP, double& area, int markb[MAXEDGES],
75 : : int markr[MAXEDGES], int& nsBlue, int& nsTgt,
76 : : bool check_boxes_first )
77 : : {
78 : : // the area will be used from now on, to see how well we fill the target cell with polygons
79 : : // the points will be at most 40; they will describe a convex patch, after the points will be
80 : : // ordered and collapsed (eliminate doubles)
81 : :
82 : : // CartVect srccoords[4];
83 : 724 : int num_nodes = 0;
84 [ + - ][ - + ]: 724 : ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
85 : 724 : nsBlue = num_nodes;
86 : : // account for possible padded polygons
87 [ - + ][ # # ]: 724 : while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
88 : 0 : nsBlue--;
89 [ + - ][ + - ]: 724 : rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
90 : :
91 : 724 : area = 0.;
92 : 724 : nP = 0; // number of intersection points we are marking the boundary of src!
93 [ + + ]: 724 : if( check_boxes_first )
94 : : {
95 : : // look at the boxes formed with vertices; if they are far away, return false early
96 : : // make sure the target is setup already
97 [ + - ]: 3 : setup_tgt_cell( tgt, nsTgt ); // we do not need area here
98 : : // use here gnomonic plane (plane) to see where source is
99 [ + - ]: 3 : bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
100 : : int planeb;
101 [ + - ][ + - ]: 3 : CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
[ + - ]
102 [ + - ]: 3 : IntxUtils::decide_gnomonic_plane( mid3, planeb );
103 [ + + ][ + - ]: 3 : if( !overlap3d && ( plane != planeb ) ) // plane was set at setup_tgt_cell
104 : 1 : return MB_SUCCESS; // no error, but no intersection, decide early to get out
105 : : // if same plane, still check for gnomonic plane in 2d
106 : : // if no overlap in 2d, get out
107 [ - + ][ # # ]: 2 : if( !overlap3d && plane == planeb ) // CHECK 2D too
108 : : {
109 [ # # ]: 0 : for( int j = 0; j < nsBlue; j++ )
110 : : {
111 : 0 : rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j],
112 [ # # ][ # # ]: 0 : srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
[ # # ][ # # ]
113 : : }
114 [ # # ]: 0 : bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
115 [ # # ]: 2 : if( !overlap2d ) return MB_SUCCESS; // we are sure they are not overlapping in 2d , either
116 : : }
117 : : }
118 : : #ifdef ENABLE_DEBUG
119 : : if( dbg_1 )
120 : : {
121 : : std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
122 : : for( int j = 0; j < nsTgt; j++ )
123 : : {
124 : : std::cout << tgtCoords[j] << "\n";
125 : : }
126 : : std::cout << "src " << mb->id_from_handle( src ) << "\n";
127 : : for( int j = 0; j < nsBlue; j++ )
128 : : {
129 : : std::cout << srcCoords[j] << "\n";
130 : : }
131 : : mb->list_entities( &tgt, 1 );
132 : : mb->list_entities( &src, 1 );
133 : : }
134 : : #endif
135 : :
136 [ + + ]: 3615 : for( int j = 0; j < nsBlue; j++ )
137 : : {
138 [ + - ][ - + ]: 2892 : rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
[ # # ][ # # ]
139 : : }
140 : :
141 : : #ifdef ENABLE_DEBUG
142 : : if( dbg_1 )
143 : : {
144 : : std::cout << "gnomonic plane: " << plane << "\n";
145 : : std::cout << " target src\n";
146 : : for( int j = 0; j < nsTgt; j++ )
147 : : {
148 : : std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
149 : : }
150 : : for( int j = 0; j < nsBlue; j++ )
151 : : {
152 : : std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
153 : : }
154 : : }
155 : : #endif
156 : :
157 [ + - ][ - + ]: 723 : rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
[ # # ][ # # ]
158 : :
159 : 723 : int side[MAXEDGES] = { 0 }; // this refers to what side? source or tgt?
160 : : int extraPoints =
161 [ + - ]: 723 : IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
162 [ + + ]: 723 : if( extraPoints >= 1 )
163 : : {
164 [ + + ]: 1360 : for( int k = 0; k < nsBlue; k++ )
165 : : {
166 [ + + ]: 1088 : if( side[k] )
167 : : {
168 : : // this means that vertex k of source is inside convex tgt; mark edges k-1 and k in
169 : : // src,
170 : : // as being "intersected" by tgt; (even though they might not be intersected by
171 : : // other edges, the fact that their apex is inside, is good enough)
172 : 272 : markb[k] = 1;
173 : 272 : markb[( k + nsBlue - 1 ) % nsBlue] =
174 : 272 : 1; // it is the previous edge, actually, but instead of doing -1, it is
175 : : // better to do modulo +3 (modulo 4)
176 : : // null side b for next call
177 : 272 : side[k] = 0;
178 : : }
179 : : }
180 : : }
181 : 723 : nP += extraPoints;
182 : :
183 : : extraPoints =
184 [ + - ]: 723 : IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
185 [ + + ]: 723 : if( extraPoints >= 1 )
186 : : {
187 [ + + ]: 3535 : for( int k = 0; k < nsTgt; k++ )
188 : : {
189 [ + + ]: 2828 : if( side[k] )
190 : : {
191 : : // this is to mark that target edges k-1 and k are intersecting src
192 : 1350 : markr[k] = 1;
193 : 1350 : markr[( k + nsTgt - 1 ) % nsTgt] =
194 : 1350 : 1; // it is the previous edge, actually, but instead of doing -1, it is
195 : : // better to do modulo +3 (modulo 4)
196 : : // null side b for next call
197 : : }
198 : : }
199 : : }
200 : 723 : nP += extraPoints;
201 : :
202 : : // now sort and orient the points in P, such that they are forming a convex polygon
203 : : // this will be the foundation of our new mesh
204 : : // this works if the polygons are convex
205 [ + - ]: 723 : IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 ); // nP should be at most 8 in the end ?
206 : : // if there are more than 3 points, some area will be positive
207 : :
208 [ + + ]: 723 : if( nP >= 3 )
209 : : {
210 [ + + ]: 2200 : for( int k = 1; k < nP - 1; k++ )
211 [ + - ]: 1478 : area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
212 : : #ifdef CHECK_CONVEXITY
213 : : // each edge should be large enough that we can compute angles between edges
214 : : for( int k = 0; k < nP; k++ )
215 : : {
216 : : int k1 = ( k + 1 ) % nP;
217 : : int k2 = ( k1 + 1 ) % nP;
218 : : double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
219 : : if( orientedArea < 0 )
220 : : {
221 : : std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
222 : : << " " << src << " \n";
223 : : }
224 : : }
225 : : #endif
226 : : }
227 : :
228 : 724 : return MB_SUCCESS; // no error
229 : : }
230 : :
231 : : // this method will also construct the triangles/quads/polygons in the new mesh
232 : : // if we accept planar polygons, we just save them
233 : : // also, we could just create new vertices every time, and merge only in the end;
234 : : // could be too expensive, and the tolerance for merging could be an
235 : : // interesting topic
236 : 506 : ErrorCode Intx2MeshOnSphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsBlue, double* iP, int nP )
237 : : {
238 : : #ifdef ENABLE_DEBUG
239 : : // first of all, check against target and source vertices
240 : : //
241 : : if( dbg_1 )
242 : : {
243 : : std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
244 : : << "\n";
245 : : for( int n = 0; n < nP; n++ )
246 : : std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
247 : : }
248 : : #endif
249 : :
250 : : // get the edges for the target triangle; the extra points will be on those edges, saved as
251 : : // lists (unordered)
252 : :
253 : : // first get the list of edges adjacent to the target cell
254 : : // use the neighTgtEdgeTag
255 : : EntityHandle adjTgtEdges[MAXEDGES];
256 [ + - ][ - + ]: 506 : ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge target tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
257 : : // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
258 : : // potatoes some of them will be handles to the initial vertices from source or target meshes
259 : :
260 [ + - ]: 506 : std::vector< EntityHandle > foundIds;
261 [ + - ]: 506 : foundIds.resize( nP );
262 : : #ifdef CHECK_CONVEXITY
263 : : int npBefore1 = nP;
264 : : int oldNodes = 0;
265 : : int otherIntx = 0;
266 : : #endif
267 [ + + ]: 2538 : for( int i = 0; i < nP; i++ )
268 : : {
269 : 2032 : double* pp = &iP[2 * i]; // iP+2*i
270 : : // project the point back on the sphere
271 [ + - ]: 2032 : CartVect pos;
272 [ + - ]: 2032 : IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
273 : 2032 : int found = 0;
274 : : // first, are they on vertices from target or src?
275 : : // priority is the target mesh (mb2?)
276 : 2032 : int j = 0;
277 : 2032 : EntityHandle outNode = (EntityHandle)0;
278 [ + + ][ + + ]: 8864 : for( j = 0; j < nsTgt && !found; j++ )
279 : : {
280 : : // int node = tgtTri.v[j];
281 [ + - ]: 6832 : double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
282 [ + + ]: 6832 : if( d2 < epsilon_1 )
283 : : {
284 : :
285 [ + - ]: 864 : foundIds[i] = tgtConn[j]; // no new node
286 : 864 : found = 1;
287 : : #ifdef CHECK_CONVEXITY
288 : : oldNodes++;
289 : : #endif
290 : : #ifdef ENABLE_DEBUG
291 : : if( dbg_1 )
292 : : std::cout << " target node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
293 : : << " 2d coords:" << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << " d2: " << d2
294 : : << " \n";
295 : : #endif
296 : : }
297 : : }
298 : :
299 [ + + ][ + + ]: 6380 : for( j = 0; j < nsBlue && !found; j++ )
300 : : {
301 : : // int node = srcTri.v[j];
302 [ + - ]: 4348 : double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
303 [ + + ]: 4348 : if( d2 < epsilon_1 )
304 : : {
305 : : // suspect is srcConn[j] corresponding in mbOut
306 : :
307 [ + - ]: 216 : foundIds[i] = srcConn[j]; // no new node
308 : 216 : found = 1;
309 : : #ifdef CHECK_CONVEXITY
310 : : oldNodes++;
311 : : #endif
312 : : #ifdef ENABLE_DEBUG
313 : : if( dbg_1 )
314 : : std::cout << " source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
315 : : << " \n";
316 : : #endif
317 : : }
318 : : }
319 : :
320 [ + + ]: 2032 : if( !found )
321 : : {
322 : : // find the edge it belongs, first, on the red element
323 : : // look at the minimum area, not at the first below some tolerance
324 : 952 : double minArea = 1.e+38;
325 : 952 : int index_min = -1;
326 [ + + ]: 4760 : for( j = 0; j < nsTgt; j++ )
327 : : {
328 : 3808 : int j1 = ( j + 1 ) % nsTgt;
329 [ + - ]: 3808 : double area = fabs( IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp ) );
330 : : // how to check if pp is between redCoords2D[j] and redCoords2D[j1] ?
331 : : // they should form a straight line; the sign should be -1
332 [ + - ]: 3808 : double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
333 [ + - ]: 3808 : IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
334 [ + - ]: 3808 : IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
335 [ + + ][ + + ]: 3808 : if( area < minArea && checkx < 2 * epsilon_1 ) // round off error or not?
336 : : {
337 : 952 : index_min = j;
338 : 952 : minArea = area;
339 : : }
340 : : }
341 [ + - ]: 952 : if( minArea < epsilon_1 / 2 ) // we found the smallest area, so we think we found the
342 : : // target edge it belongs
343 : : {
344 : : // found the edge; now find if there is a point in the list here
345 : : // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
346 [ + - ]: 952 : int indx = TgtEdges.index( adjTgtEdges[index_min] );
347 [ - + ]: 952 : if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
348 : : {
349 [ # # ][ # # ]: 0 : std::cerr << " error in adjacent target edge: " << mb->id_from_handle( adjTgtEdges[index_min] )
[ # # ]
350 [ # # ]: 0 : << "\n";
351 : 0 : return MB_FAILURE;
352 : : }
353 [ + - ]: 952 : std::vector< EntityHandle >* expts = extraNodesVec[indx];
354 : : // if the points pp is between extra points, then just give that id
355 : : // if not, create a new point, (check the id)
356 : : // get the coordinates of the extra points so far
357 : 952 : int nbExtraNodesSoFar = expts->size();
358 [ + + ]: 952 : if( nbExtraNodesSoFar > 0 )
359 : : {
360 [ + - ]: 729 : std::vector< CartVect > coords1;
361 [ + - ]: 729 : coords1.resize( nbExtraNodesSoFar );
362 [ + - ][ + - ]: 729 : mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
[ + - ][ + - ]
363 : : // std::list<int>::iterator it;
364 [ + + ][ + + ]: 1503 : for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
365 : : {
366 : : // int pnt = *it;
367 [ + - ][ + - ]: 774 : double d2 = ( pos - coords1[k] ).length();
[ + - ]
368 [ + + ]: 774 : if( d2 < 2 * epsilon_1 ) // is this below machine precision?
369 : : {
370 : 714 : found = 1;
371 [ + - ][ + - ]: 714 : foundIds[i] = ( *expts )[k];
372 : : #ifdef CHECK_CONVEXITY
373 : : otherIntx++;
374 : : #endif
375 : : }
376 : 729 : }
377 : : }
378 [ + + ]: 952 : if( !found )
379 : : {
380 : : // create a new point in 2d (at the intersection)
381 : : // foundIds[i] = m_num2dPoints;
382 : : // expts.push_back(m_num2dPoints);
383 : : // need to create a new node in mbOut
384 : : // this will be on the edge, and it will be added to the local list
385 [ + - ][ + - ]: 238 : rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
386 [ + - ]: 238 : ( *expts ).push_back( outNode );
387 : : // CID 181168; avoid leak storage error
388 [ + - ][ - + ]: 238 : rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
389 [ + - ]: 238 : foundIds[i] = outNode;
390 : 952 : found = 1;
391 : : }
392 : : }
393 : : }
394 [ - + ]: 2032 : if( !found )
395 : : {
396 [ # # ]: 0 : std::cout << " target quad: ";
397 [ # # ]: 0 : for( int j1 = 0; j1 < nsTgt; j1++ )
398 : : {
399 [ # # ][ # # ]: 0 : std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
[ # # ][ # # ]
400 : : }
401 [ # # ][ # # ]: 0 : std::cout << " a point pp is not on a target quad " << *pp << " " << pp[1] << " target quad "
[ # # ][ # # ]
[ # # ]
402 [ # # ][ # # ]: 0 : << mb->id_from_handle( tgt ) << " \n";
[ # # ]
403 : 0 : return MB_FAILURE;
404 : : }
405 : : }
406 : : #ifdef ENABLE_DEBUG
407 : : if( dbg_1 )
408 : : {
409 : : std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
410 : : for( int i1 = 0; i1 < nP; i1++ )
411 : : std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
412 : : }
413 : : #endif
414 : : // first, find out if we have nodes collapsed; shrink them
415 : : // we may have to reduce nP
416 : : // it is possible that some nodes are collapsed after intersection only
417 : : // nodes will always be in order (convex intersection)
418 : : #ifdef CHECK_CONVEXITY
419 : : int npBefore2 = nP;
420 : : #endif
421 [ + - ][ + - ]: 506 : correct_polygon( &foundIds[0], nP );
422 : : // now we can build the triangles, from P array, with foundIds
423 : : // we will put them in the out set
424 [ + - ]: 506 : if( nP >= 3 )
425 : : {
426 : : EntityHandle polyNew;
427 [ + - ][ + - ]: 506 : rval = mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
428 [ + - ][ - + ]: 506 : rval = mb->add_entities( outSet, &polyNew, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
429 : :
430 : : // tag it with the global ids from target and source elements
431 : : int globalID;
432 [ + - ][ - + ]: 506 : rval = mb->tag_get_data( gid, &src, 1, &globalID );MB_CHK_ERR( rval );
[ # # ][ # # ]
433 [ + - ][ - + ]: 506 : rval = mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
[ # # ][ # # ]
434 : : // if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
435 : : // " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";
436 [ + - ][ - + ]: 506 : rval = mb->tag_get_data( gid, &tgt, 1, &globalID );MB_CHK_ERR( rval );
[ # # ][ # # ]
437 [ + - ][ - + ]: 506 : rval = mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
[ # # ][ # # ]
438 : : // if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<
439 : : // " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";
440 : :
441 : 506 : counting++;
442 [ + - ][ - + ]: 506 : rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
[ # # ][ # # ]
443 [ - + ]: 506 : if( orgSendProcTag )
444 : : {
445 : 0 : int org_proc = -1;
446 [ # # ][ # # ]: 506 : rval = mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc );MB_CHK_ERR( rval );
[ # # ][ # # ]
447 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc );MB_CHK_ERR( rval ); // yet another tag
[ # # ][ # # ]
448 : : }
449 : : #ifdef CHECK_CONVEXITY
450 : : // each edge should be large enough that we can compute angles between edges
451 : : std::vector< double > coords;
452 : : coords.resize( 3 * nP );
453 : : rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
454 : : std::vector< CartVect > posi( nP );
455 : : rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
456 : :
457 : : for( int k = 0; k < nP; k++ )
458 : : {
459 : : int k1 = ( k + 1 ) % nP;
460 : : int k2 = ( k1 + 1 ) % nP;
461 : : double orientedArea =
462 : : area_spherical_triangle_lHuiller( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
463 : : if( orientedArea < 0 )
464 : : {
465 : : std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
466 : : for( int i = 0; i < nP; i++ )
467 : : {
468 : : int nexti = ( i + 1 ) % nP;
469 : : double lengthEdge = ( posi[i] - posi[nexti] ).length();
470 : : std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n";
471 : : }
472 : : std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n";
473 : :
474 : : std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
475 : : << " target, src:" << tgt << " " << src << " \n";
476 : : }
477 : : }
478 : : #endif
479 : :
480 : : #ifdef ENABLE_DEBUG
481 : : if( dbg_1 )
482 : : {
483 : : std::cout << "Counting: " << counting << "\n";
484 : : std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
485 : : for( int i1 = 0; i1 < nP; i1++ )
486 : : std::cout << " " << mb->id_from_handle( foundIds[i1] );
487 : : std::cout << " plane: " << plane << "\n";
488 : : std::vector< CartVect > posi( nP );
489 : : mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
490 : : for( int i1 = 0; i1 < nP; i1++ )
491 : : std::cout << foundIds[i1] << " " << posi[i1] << "\n";
492 : :
493 : : std::stringstream fff;
494 : : fff << "file0" << counting << ".vtk";
495 : : rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
496 : : }
497 : : #endif
498 : : }
499 : : // else {
500 : : // std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";
501 : : // }
502 : : // disable_debug();
503 : 506 : return MB_SUCCESS;
504 : : }
505 : :
506 : 0 : ErrorCode Intx2MeshOnSphere::update_tracer_data( EntityHandle out_set, Tag& tagElem, Tag& tagArea )
507 : : {
508 : 0 : EntityHandle dum = 0;
509 : :
510 : : Tag corrTag;
511 : : ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE,
512 [ # # ]: 0 : &dum ); // it should have been created
513 [ # # ][ # # ]: 0 : MB_CHK_SET_ERR( rval, "can't get correlation tag" );
[ # # ][ # # ]
[ # # ][ # # ]
514 : :
515 : : // get all polygons out of out_set; then see where are they coming from
516 [ # # ]: 0 : Range polys;
517 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
518 : :
519 : : // rs2 is the target range, arrival; rs1 is src, departure;
520 : : // there is a connection between rs1 and rs2, through the corrTag
521 : : // corrTag is __correlation
522 : : // basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly);
523 : : // also, mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly);
524 : : // we start from rs2 existing, then we have to update something
525 : :
526 : : // tagElem will have multiple tracers
527 : 0 : int numTracers = 0;
528 [ # # ][ # # ]: 0 : rval = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
529 [ # # ][ # # ]: 0 : if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
[ # # ][ # # ]
[ # # ][ # # ]
530 : :
531 [ # # ][ # # ]: 0 : std::vector< double > currentVals( rs2.size() * numTracers );
532 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( tagElem, rs2, ¤tVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
533 : :
534 : : // create new tuple list for tracers to other processors, from remote_cells
535 : : #ifdef MOAB_HAVE_MPI
536 [ # # ]: 0 : if( remote_cells )
537 : : {
538 [ # # ]: 0 : int n = remote_cells->get_n();
539 [ # # ]: 0 : if( n > 0 )
540 : : {
541 [ # # ][ # # ]: 0 : remote_cells_with_tracers = new TupleList();
542 : : remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
543 [ # # ]: 0 : n ); // tracers are in these tuples
544 [ # # ]: 0 : remote_cells_with_tracers->enableWriteAccess();
545 [ # # ]: 0 : for( int i = 0; i < n; i++ )
546 : : {
547 : 0 : remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i];
548 : 0 : remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
549 : : // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
550 : 0 : remote_cells_with_tracers->vul_wr[i] =
551 : 0 : remote_cells->vul_wr[i]; // this is the corresponding target cell (arrival)
552 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
553 : 0 : remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0; // initialize tracers to be transported
554 [ # # ]: 0 : remote_cells_with_tracers->inc_n();
555 : : }
556 : : }
557 [ # # ]: 0 : delete remote_cells;
558 : 0 : remote_cells = NULL;
559 : : }
560 : : #endif
561 : : // for each polygon, we have 2 indices: target and source parents
562 : : // we need index source to update index tgt?
563 [ # # ]: 0 : std::vector< double > newValues( rs2.size() * numTracers,
564 [ # # ]: 0 : 0. ); // initialize with 0 all of them
565 : : // area of the polygon * conc on target (old) current quantity
566 : : // finally, divide by the area of the tgt
567 : 0 : double check_intx_area = 0.;
568 [ # # ]: 0 : moab::IntxAreaUtils intxAreas( this->areaMethod ); // use_lHuiller = true
569 [ # # ][ # # ]: 0 : for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
[ # # ][ # # ]
[ # # ]
570 : : {
571 [ # # ]: 0 : EntityHandle poly = *it;
572 : : int srcIndex, tgtIndex;
573 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
574 : :
575 [ # # ]: 0 : EntityHandle src = rs1[srcIndex - 1]; // big assumption, it should work for meshes where global id is the same
576 : : // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes
577 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
578 : : // EntityHandle target = rs2[tgtIndex];
579 : : // big assumption here, target and source are "parallel" ;we should have an index from
580 : : // source to target (so a deformed source corresponds to an arrival tgt)
581 : : /// TODO: VSM: Its unclear whether we need the source or destination radius here.
582 : 0 : double radius = Rsrc;
583 [ # # ]: 0 : double areap = intxAreas.area_spherical_element( mb, poly, radius );
584 : 0 : check_intx_area += areap;
585 : : // so the departure cell at time t (srcIndex) covers a portion of a tgtCell
586 : : // that quantity will be transported to the tgtCell at time t+dt
587 : : // the source corresponds to a target arrival
588 : : EntityHandle tgtArr;
589 [ # # ]: 0 : rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
590 [ # # ][ # # ]: 0 : if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
591 : : {
592 : : #ifdef MOAB_HAVE_MPI
593 [ # # ][ # # ]: 0 : if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
[ # # ][ # # ]
[ # # ][ # # ]
594 : : // maybe the element is remote, from another processor
595 : : int global_id_src;
596 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &src, 1, &global_id_src );MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding source gid" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
597 : : // find the
598 [ # # ]: 0 : int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
599 [ # # ]: 0 : if( index_in_remote == -1 )
600 [ # # ][ # # ]: 0 : MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
[ # # ][ # # ]
[ # # ]
601 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
602 : : remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
603 [ # # ]: 0 : currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
604 : : #endif
605 : : }
606 [ # # ]: 0 : else if( MB_SUCCESS == rval )
607 : : {
608 [ # # ]: 0 : int arrTgtIndex = rs2.index( tgtArr );
609 [ # # ][ # # ]: 0 : if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
[ # # ][ # # ]
[ # # ][ # # ]
610 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
611 [ # # ][ # # ]: 0 : newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
612 : : }
613 : :
614 : : else
615 [ # # ][ # # ]: 0 : MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " );
[ # # ][ # # ]
[ # # ][ # # ]
616 : : }
617 : : // now, send back the remote_cells_with_tracers to the processors they came from, with the
618 : : // updated values for the tracer mass in a cell
619 : : #ifdef MOAB_HAVE_MPI
620 [ # # ]: 0 : if( remote_cells_with_tracers )
621 : : {
622 : : // so this means that some cells will be sent back with tracer info to the procs they were
623 : : // sent from
624 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
[ # # ]
625 : : // now, look at the global id, find the proper "tgt" cell with that index and update its
626 : : // mass
627 : : // remote_cells->print("remote cells after routing");
628 [ # # ]: 0 : int n = remote_cells_with_tracers->get_n();
629 [ # # ]: 0 : for( int j = 0; j < n; j++ )
630 : : {
631 : 0 : EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j]; // entity handle sent back
632 [ # # ]: 0 : int arrTgtIndex = rs2.index( tgtCell );
633 [ # # ][ # # ]: 0 : if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
[ # # ][ # # ]
[ # # ][ # # ]
634 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
635 [ # # ]: 0 : newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
636 : : }
637 : : }
638 : : #endif /* MOAB_HAVE_MPI */
639 : : // now divide by target area (current)
640 : 0 : int j = 0;
641 [ # # ]: 0 : Range::iterator iter = rs2.begin();
642 : 0 : void* data = NULL; // used for stored area
643 : 0 : int count = 0;
644 [ # # ]: 0 : std::vector< double > total_mass_local( numTracers, 0. );
645 [ # # ][ # # ]: 0 : while( iter != rs2.end() )
[ # # ]
646 : : {
647 [ # # ][ # # ]: 0 : rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
648 : 0 : double* ptrArea = (double*)data;
649 [ # # ][ # # ]: 0 : for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
650 : : {
651 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
652 : : {
653 [ # # ][ # # ]: 0 : total_mass_local[k] += newValues[j * numTracers + k];
654 [ # # ]: 0 : newValues[j * numTracers + k] /= ( *ptrArea );
655 : : }
656 : : }
657 : : }
658 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
659 : :
660 : : #ifdef MOAB_HAVE_MPI
661 [ # # ]: 0 : std::vector< double > total_mass( numTracers, 0. );
662 : 0 : double total_intx_area = 0;
663 : : int mpi_err =
664 [ # # ][ # # ]: 0 : MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
[ # # ]
665 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
666 : : // now reduce total area
667 [ # # ]: 0 : mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
668 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
669 [ # # ]: 0 : if( my_rank == 0 )
670 : : {
671 [ # # ]: 0 : for( int k = 0; k < numTracers; k++ )
672 [ # # ][ # # ]: 0 : std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
673 [ # # ][ # # ]: 0 : std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
[ # # ]
674 [ # # ][ # # ]: 0 : << total_intx_area << "\n";
675 : : }
676 : :
677 [ # # ]: 0 : if( remote_cells_with_tracers )
678 : : {
679 [ # # ]: 0 : delete remote_cells_with_tracers;
680 : 0 : remote_cells_with_tracers = NULL;
681 : : }
682 : : #else
683 : : for( int k = 0; k < numTracers; k++ )
684 : : std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n";
685 : : std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "
686 : : << check_intx_area << "\n";
687 : : #endif
688 : 0 : return MB_SUCCESS;
689 : : }
690 : :
691 : : #ifdef MOAB_HAVE_MPI
692 : 0 : ErrorCode Intx2MeshOnSphere::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
693 : : {
694 [ # # ]: 0 : localEnts.clear();
695 [ # # ][ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );MB_CHK_SET_ERR( rval, "can't get local ents" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
696 : :
697 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( localEnts, local_verts );MB_CHK_SET_ERR( rval, "can't get connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
698 [ # # ]: 0 : int num_local_verts = (int)local_verts.size();
699 : :
700 [ # # ]: 0 : assert( parcomm != NULL );
701 : :
702 [ # # ]: 0 : if( num_local_verts == 0 )
703 : : {
704 : : // it is probably point cloud, get the local vertices from set
705 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( euler_set, 0, local_verts );MB_CHK_SET_ERR( rval, "can't get local vertices from set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
706 [ # # ]: 0 : num_local_verts = (int)local_verts.size();
707 [ # # ]: 0 : localEnts = local_verts;
708 : : }
709 : : // will use 6 gnomonic planes to decide boxes for each gnomonic plane
710 : : // each gnomonic box will be 2d, min, max
711 : : double gnom_box[24];
712 [ # # ]: 0 : for( int i = 0; i < 6; i++ )
713 : : {
714 : 0 : gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
715 : 0 : gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
716 : : }
717 : :
718 : : // there are 6 gnomonic planes; some elements could be on the corners, and affect multiple
719 : : // planes decide what gnomonic planes will be affected by each cell some elements could appear
720 : : // in multiple gnomonic planes !
721 [ # # ]: 0 : std::vector< double > coords( 3 * num_local_verts );
722 [ # # ][ # # ]: 0 : rval = mb->get_coords( local_verts, &coords[0] );MB_CHK_SET_ERR( rval, "can't get vertex coords" );ERRORR( rval, "can't get coords of vertices " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
723 : : // decide each local vertex to what gnomonic plane it belongs
724 : :
725 [ # # ]: 0 : std::vector< int > gnplane;
726 [ # # ]: 0 : gnplane.resize( num_local_verts );
727 [ # # ]: 0 : for( int i = 0; i < num_local_verts; i++ )
728 : : {
729 [ # # ][ # # ]: 0 : CartVect pos( &coords[3 * i] );
730 : : int pl;
731 [ # # ]: 0 : IntxUtils::decide_gnomonic_plane( pos, pl );
732 [ # # ]: 0 : gnplane[i] = pl;
733 : : }
734 : :
735 [ # # ][ # # ]: 0 : for( Range::iterator it = localEnts.begin(); it != localEnts.end(); it++ )
[ # # ][ # # ]
[ # # ]
736 : : {
737 [ # # ]: 0 : EntityHandle cell = *it;
738 [ # # ]: 0 : EntityType typeCell = mb->type_from_handle( cell ); // could be vertex, for point cloud
739 : : // get coordinates, and decide gnomonic planes for it
740 : : int nnodes;
741 : 0 : const EntityHandle* conn = NULL;
742 : : EntityHandle c[1];
743 [ # # ]: 0 : if( typeCell != MBVERTEX )
744 : : {
745 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
746 : : }
747 : : else
748 : : {
749 : 0 : nnodes = 1;
750 : 0 : c[0] = cell; // actual node
751 : 0 : conn = &c[0];
752 : : }
753 : : // get coordinates of vertices involved with this
754 [ # # ]: 0 : std::vector< double > elco( 3 * nnodes );
755 [ # # ]: 0 : std::set< int > planes;
756 [ # # ]: 0 : for( int i = 0; i < nnodes; i++ )
757 : : {
758 [ # # ]: 0 : int ix = local_verts.index( conn[i] );
759 [ # # ][ # # ]: 0 : planes.insert( gnplane[ix] );
760 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
761 : : {
762 [ # # ][ # # ]: 0 : elco[3 * i + j] = coords[3 * ix + j];
763 : : }
764 : : }
765 : : // now, augment the boxes for all planes involved
766 [ # # ][ # # ]: 0 : for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
[ # # ]
767 : : {
768 [ # # ]: 0 : int pl = *st;
769 [ # # ]: 0 : for( int i = 0; i < nnodes; i++ )
770 : : {
771 [ # # ][ # # ]: 0 : CartVect pos( &elco[3 * i] );
772 : : double c2[2];
773 : : IntxUtils::gnomonic_projection( pos, Rdest, pl, c2[0],
774 [ # # ]: 0 : c2[1] ); // 2 coordinates
775 : : //
776 [ # # ]: 0 : for( int k = 0; k < 2; k++ )
777 : : {
778 : 0 : double val = c2[k];
779 [ # # ]: 0 : if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val; // min in k direction
780 [ # # ]: 0 : if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
781 : 0 : gnom_box[4 * ( pl - 1 ) + 2 + k] = val; // max in k direction
782 : : }
783 : : }
784 : : }
785 : 0 : }
786 : :
787 [ # # ][ # # ]: 0 : int numprocs = parcomm->proc_config().proc_size();
788 [ # # ]: 0 : allBoxes.resize( 24 * numprocs ); // 6 gnomonic planes , 4 doubles for each for 2d box
789 : :
790 [ # # ][ # # ]: 0 : my_rank = parcomm->proc_config().proc_rank();
791 [ # # ]: 0 : for( int k = 0; k < 24; k++ )
792 [ # # ]: 0 : allBoxes[24 * my_rank + k] = gnom_box[k];
793 : :
794 : : // now communicate to get all boxes
795 : : int mpi_err;
796 : : #if( MPI_VERSION >= 2 )
797 : : // use "in place" option
798 [ # # ]: 0 : mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 24, MPI_DOUBLE,
799 [ # # ][ # # ]: 0 : parcomm->proc_config().proc_comm() );
[ # # ]
800 : : #else
801 : : {
802 : : std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
803 : : mpi_err = MPI_Allgather( &allBoxes[24 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
804 : : parcomm->proc_config().proc_comm() );
805 : : allBoxes = allBoxes_tmp;
806 : : }
807 : : #endif
808 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
809 : :
810 : : #ifdef VERBOSE
811 : : if( my_rank == 0 )
812 : : {
813 : : std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
814 : : << " on second mesh \n";
815 : : for( int i = 0; i < numprocs; i++ )
816 : : {
817 : : std::cout << "task: " << i << " \n";
818 : : for( int pl = 1; pl <= 6; pl++ )
819 : : {
820 : : std::cout << " plane " << pl << " min: \t" << allBoxes[24 * i + 4 * ( pl - 1 )] << " \t"
821 : : << allBoxes[24 * i + 4 * ( pl - 1 ) + 1] << "\n";
822 : : std::cout << " \t max: \t" << allBoxes[24 * i + 4 * ( pl - 1 ) + 2] << " \t"
823 : : << allBoxes[24 * i + 4 * ( pl - 1 ) + 3] << "\n";
824 : : }
825 : : }
826 : : }
827 : : #endif
828 : :
829 : 0 : return MB_SUCCESS;
830 : : }
831 : : //#define VERBOSE
832 : : // this will use the bounding boxes for the (euler)/ fix mesh that are already established
833 : : // will distribute the mesh to other procs, so that on each task, the covering set covers the local
834 : : // bounding box this means it will cover the second (local) mesh set; So the covering set will cover
835 : : // completely the second local mesh set (in intersection)
836 : 0 : ErrorCode Intx2MeshOnSphere::construct_covering_set( EntityHandle& initial_distributed_set, EntityHandle& covering_set )
837 : : {
838 [ # # ]: 0 : assert( parcomm != NULL );
839 [ # # ][ # # ]: 0 : if( 1 == parcomm->proc_config().proc_size() )
[ # # ]
840 : : {
841 : 0 : covering_set = initial_distributed_set; // nothing to move around, it must be serial
842 : 0 : return MB_SUCCESS;
843 : : }
844 : :
845 : : // primary element came from, in the joint communicator ; this will be forwarded by coverage
846 : : // mesh needed for tag migrate later on
847 : 0 : int defaultInt = -1; // no processor, so it was not migrated from somewhere else
848 : : ErrorCode rval = mb->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag,
849 [ # # ][ # # ]: 0 : MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
850 : :
851 : : // mark on the coverage mesh where this element came from
852 : : Tag sendProcTag; /// for coverage mesh, will store the sender
853 : : rval = mb->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag, MB_TAG_DENSE | MB_TAG_CREAT,
854 [ # # ][ # # ]: 0 : &defaultInt );MB_CHK_SET_ERR( rval, "can't create sending processor tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
855 : :
856 : : // this information needs to be forwarded to coverage mesh, if this mesh was already migrated
857 : : // from somewhere else
858 [ # # ]: 0 : Range meshCells;
859 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( initial_distributed_set, 2, meshCells );MB_CHK_SET_ERR( rval, "can't get cells by dimension from mesh set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
860 : :
861 : : // look at the value of orgSendProcTag for one mesh cell; if -1, no need to forward that; if
862 : : // !=-1, we know that this mesh was migrated, we need to find out more about origin of cell
863 : 0 : int orig_sender = -1;
864 : 0 : EntityHandle oneCell = 0;
865 : : // decide if we need to transfer global DOFs info attached to each HOMME coarse cell; first we
866 : : // need to decide if the mesh has that tag; will affect the size of the tuple list involved in
867 : : // the crystal routing
868 : 0 : int size_gdofs_tag = 0;
869 [ # # ]: 0 : std::vector< int > valsDOFs;
870 : : Tag gdsTag;
871 [ # # ]: 0 : rval = mb->tag_get_handle( "GLOBAL_DOFS", gdsTag );
872 : :
873 [ # # ][ # # ]: 0 : if( meshCells.size() > 0 )
874 : : {
875 [ # # ]: 0 : oneCell = meshCells[0]; // it is possible we do not have any cells, even after migration
876 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( orgSendProcTag, &oneCell, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sending processor value" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
877 [ # # ]: 0 : if( gdsTag )
878 : : {
879 : : DataType dtype;
880 [ # # ]: 0 : rval = mb->tag_get_data_type( gdsTag, dtype );
881 [ # # ][ # # ]: 0 : if( MB_SUCCESS == rval && MB_TYPE_INTEGER == dtype )
882 : : {
883 : : // find the values on first cell
884 : 0 : int lenTag = 0;
885 [ # # ]: 0 : rval = mb->tag_get_length( gdsTag, lenTag );
886 [ # # ][ # # ]: 0 : if( MB_SUCCESS == rval && lenTag > 0 )
887 : : {
888 [ # # ]: 0 : valsDOFs.resize( lenTag );
889 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gdsTag, &oneCell, 1, &valsDOFs[0] );
890 [ # # ][ # # ]: 0 : if( MB_SUCCESS == rval && valsDOFs[0] > 0 )
[ # # ][ # # ]
891 : : {
892 : : // first value positive means we really need to transport this data during
893 : : // coverage
894 : 0 : size_gdofs_tag = lenTag;
895 : : }
896 : : }
897 : : }
898 : : }
899 : : }
900 : :
901 : : // another collective call, to see if the mesh is migrated and if the GLOBAL_DOFS tag need to be
902 : : // transferred over to the coverage mesh it is possible that there is no initial mesh source
903 : : // mesh on the task, so we do not know that info from the tag but TupleList needs to be sized
904 : : // uniformly for all tasks do a collective MPI_MAX to see if it is migrated and if we have
905 : : // (collectively) a GLOBAL_DOFS task
906 : :
907 : : int local_int_array[2], global_int_array[2];
908 : 0 : local_int_array[0] = orig_sender;
909 : 0 : local_int_array[1] = size_gdofs_tag;
910 : : // now reduce over all processors
911 : : int mpi_err =
912 [ # # ][ # # ]: 0 : MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
[ # # ]
913 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
914 : 0 : orig_sender = global_int_array[0];
915 : 0 : size_gdofs_tag = global_int_array[1];
916 : : #ifdef VERBOSE
917 : : std::cout << "proc: " << my_rank << " size_gdofs_tag:" << size_gdofs_tag << "\n";
918 : : #endif
919 [ # # ]: 0 : valsDOFs.resize( size_gdofs_tag );
920 : :
921 : : // finally, we have correct global info, we can decide if mesh was migrated and if we have
922 : : // global dofs tag that need to be sent with coverage info
923 : 0 : int migrated_mesh = 0;
924 [ # # ]: 0 : if( orig_sender != -1 ) migrated_mesh = 1; //
925 : : // if size_gdofs_tag>0, we are sure valsDOFs got resized to what we need
926 : :
927 : : // get all mesh verts1
928 [ # # ]: 0 : Range mesh_verts;
929 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( meshCells, mesh_verts );MB_CHK_SET_ERR( rval, "can't get mesh vertices" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
930 [ # # ]: 0 : int num_mesh_verts = (int)mesh_verts.size();
931 : :
932 : : // now see the mesh points positions; to what boxes should we send them?
933 [ # # ]: 0 : std::vector< double > coords_mesh( 3 * num_mesh_verts );
934 [ # # ][ # # ]: 0 : rval = mb->get_coords( mesh_verts, &coords_mesh[0] );MB_CHK_SET_ERR( rval, "can't get mesh points position" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
935 : :
936 : : // decide gnomonic plane for each vertex, as in the compute boxes
937 [ # # ]: 0 : std::vector< int > gnplane;
938 [ # # ]: 0 : gnplane.resize( num_mesh_verts );
939 [ # # ]: 0 : for( int i = 0; i < num_mesh_verts; i++ )
940 : : {
941 [ # # ][ # # ]: 0 : CartVect pos( &coords_mesh[3 * i] );
942 : : int pl;
943 [ # # ]: 0 : IntxUtils::decide_gnomonic_plane( pos, pl );
944 [ # # ]: 0 : gnplane[i] = pl;
945 : : }
946 : :
947 [ # # ]: 0 : std::vector< int > gids( num_mesh_verts );
948 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, mesh_verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get vertices gids" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
949 : :
950 : : // ranges to send to each processor; will hold vertices and elements (quads/ polygons)
951 : : // will look if the box of the mesh cell covers bounding box(es) (within tolerances)
952 [ # # ]: 0 : std::map< int, Range > Rto;
953 [ # # ][ # # ]: 0 : int numprocs = parcomm->proc_config().proc_size();
954 : :
955 [ # # ][ # # ]: 0 : for( Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
[ # # ][ # # ]
[ # # ]
956 : : {
957 [ # # ]: 0 : EntityHandle q = *eit;
958 : : const EntityHandle* conn;
959 : : int num_nodes;
960 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity on cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
961 : :
962 : : // first decide what planes need to consider
963 [ # # ]: 0 : std::set< int > planes; // if this list contains more than 3 planes, we have a very bad mesh!!!
964 [ # # ]: 0 : std::vector< double > elco( 3 * num_nodes );
965 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
966 : : {
967 : 0 : EntityHandle v = conn[i];
968 [ # # ]: 0 : int index = mesh_verts.index( v );
969 [ # # ][ # # ]: 0 : planes.insert( gnplane[index] );
970 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
971 : : {
972 [ # # ][ # # ]: 0 : elco[3 * i + j] = coords_mesh[3 * index + j]; // extract from coords
973 : : }
974 : : }
975 : : // now loop over all planes that need to be considered for this element
976 [ # # ][ # # ]: 0 : for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
[ # # ]
977 : : {
978 [ # # ]: 0 : int pl = *st; // gnomonic plane considered
979 : 0 : double qmin[2] = { DBL_MAX, DBL_MAX };
980 : 0 : double qmax[2] = { -DBL_MAX, -DBL_MAX };
981 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
982 : : {
983 [ # # ][ # # ]: 0 : CartVect dp( &elco[3 * i] ); // uses constructor for CartVect that takes a
984 : : // pointer to double
985 : : // gnomonic projection
986 : : double c2[2];
987 [ # # ]: 0 : IntxUtils::gnomonic_projection( dp, Rsrc, pl, c2[0], c2[1] ); // 2 coordinates
988 [ # # ]: 0 : for( int j = 0; j < 2; j++ )
989 : : {
990 [ # # ]: 0 : if( qmin[j] > c2[j] ) qmin[j] = c2[j];
991 [ # # ]: 0 : if( qmax[j] < c2[j] ) qmax[j] = c2[j];
992 : : }
993 : : }
994 : : // now decide if processor p should be interested in this cell, by looking at plane pl
995 : : // 2d box this is one of the few size n loops;
996 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ ) // each cell q can be sent to more than one processor
997 : : {
998 [ # # ]: 0 : double procMin1 = allBoxes[24 * p + 4 * ( pl - 1 )]; // these were determined before
999 : : //
1000 [ # # ]: 0 : if( procMin1 >= DBL_MAX ) // the processor has no targets on this plane
1001 : 0 : continue;
1002 [ # # ]: 0 : double procMin2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
1003 [ # # ]: 0 : double procMax1 = allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
1004 [ # # ]: 0 : double procMax2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
1005 : : // test overlap of 2d boxes
1006 [ # # ][ # # ]: 0 : if( procMin1 > qmax[0] + box_error || procMin2 > qmax[1] + box_error ) continue; //
1007 [ # # ][ # # ]: 0 : if( qmin[0] > procMax1 + box_error || qmin[1] > procMax2 + box_error ) continue;
1008 : : // good to be inserted
1009 [ # # ][ # # ]: 0 : Rto[p].insert( q );
1010 : : }
1011 : : }
1012 : 0 : }
1013 : :
1014 : : // here, we will use crystal router to send each cell to designated tasks (mesh migration)
1015 : :
1016 : : // a better implementation would be to use pcomm send / recv entities; a good test case
1017 : : // pcomm send / receives uses point to point communication, not global gather / scatter
1018 : :
1019 : : // now, build TLv and TLq (tuple list for vertices and cells, separately sent)
1020 : 0 : size_t numq = 0;
1021 : 0 : size_t numv = 0;
1022 : :
1023 : : // merge the list of vertices to be sent
1024 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ )
1025 : : {
1026 [ # # ]: 0 : if( p == (int)my_rank ) continue; // do not "send" it to current task, because it is already here
1027 [ # # ]: 0 : Range& range_to_P = Rto[p];
1028 : : // add the vertices to it
1029 [ # # ][ # # ]: 0 : if( range_to_P.empty() ) continue; // nothing to send to proc p
1030 : : #ifdef VERBOSE
1031 : : std::cout << " proc : " << my_rank << " to proc " << p << " send " << range_to_P.size() << " cells "
1032 : : << " psize: " << range_to_P.psize() << "\n";
1033 : : #endif
1034 [ # # ]: 0 : Range vertsToP;
1035 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( range_to_P, vertsToP );MB_CHK_SET_ERR( rval, "can't get connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1036 [ # # ]: 0 : numq = numq + range_to_P.size();
1037 [ # # ]: 0 : numv = numv + vertsToP.size();
1038 [ # # ][ # # ]: 0 : range_to_P.merge( vertsToP );
1039 : 0 : }
1040 : :
1041 [ # # ]: 0 : TupleList TLv; // send vertices with a different tuple list
1042 [ # # ]: 0 : TupleList TLq;
1043 [ # # ]: 0 : TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, 3 real coordinates
1044 [ # # ]: 0 : TLv.enableWriteAccess();
1045 : :
1046 : : // add also GLOBAL_DOFS info, if found on the mesh cell; it should be found only on HOMME cells!
1047 : : int sizeTuple =
1048 : 0 : 2 + max_edges_1 + migrated_mesh + size_gdofs_tag; // max edges could be up to MAXEDGES :) for polygons
1049 : : TLq.initialize( sizeTuple, 0, 0, 0,
1050 [ # # ]: 0 : numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v), plus
1051 : : // original sender if set (migrated mesh case)
1052 : : // we will not send the entity handle, global ID should be more than enough
1053 : : // we will not need more than 2B vertices
1054 : : // if we need more than 2B, we will need to use a different marker anyway (GLOBAL ID is not
1055 : : // enough then)
1056 : :
1057 [ # # ]: 0 : TLq.enableWriteAccess();
1058 : : #ifdef VERBOSE
1059 : : std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1060 : : #endif
1061 : :
1062 [ # # ]: 0 : for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1063 : : {
1064 [ # # ]: 0 : if( to_proc == (int)my_rank ) continue;
1065 [ # # ]: 0 : Range& range_to_P = Rto[to_proc];
1066 [ # # ]: 0 : Range V = range_to_P.subset_by_type( MBVERTEX );
1067 : :
1068 [ # # ][ # # ]: 0 : for( Range::iterator it = V.begin(); it != V.end(); ++it )
[ # # ][ # # ]
[ # # ]
1069 : : {
1070 [ # # ]: 0 : EntityHandle v = *it;
1071 [ # # ]: 0 : int index = mesh_verts.index( v ); //
1072 [ # # ]: 0 : assert( -1 != index );
1073 [ # # ]: 0 : int n = TLv.get_n(); // current size of tuple list
1074 : 0 : TLv.vi_wr[2 * n] = to_proc; // send to processor
1075 [ # # ]: 0 : TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the second_mesh_verts range
1076 [ # # ]: 0 : TLv.vr_wr[3 * n] = coords_mesh[3 * index]; // departure position, of the node local_verts[i]
1077 [ # # ]: 0 : TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
1078 [ # # ]: 0 : TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
1079 [ # # ]: 0 : TLv.inc_n(); // increment tuple list size
1080 : : }
1081 : : // also, prep the 2d cells for sending ...
1082 [ # # ][ # # ]: 0 : Range Q = range_to_P.subset_by_dimension( 2 );
1083 [ # # ][ # # ]: 0 : for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
1084 : : {
1085 [ # # ]: 0 : EntityHandle q = *it; // this is a second mesh cell (or src, lagrange set)
1086 : : int global_id;
1087 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_SET_ERR( rval, "can't get gid for polygon" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1088 [ # # ]: 0 : int n = TLq.get_n(); // current size
1089 : 0 : TLq.vi_wr[sizeTuple * n] = to_proc; //
1090 : 0 : TLq.vi_wr[sizeTuple * n + 1] =
1091 : 0 : global_id; // global id of element, used to identify it for debug purposes only
1092 : : const EntityHandle* conn4;
1093 : : int num_nodes; // could be up to MAXEDGES; max_edges?;
1094 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn4, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity for cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1095 [ # # ]: 0 : if( num_nodes > max_edges_1 )
1096 : : {
1097 [ # # ][ # # ]: 0 : mb->list_entities( &q, 1 );MB_CHK_SET_ERR( MB_FAILURE, "too many nodes in a cell (" << num_nodes << "," << max_edges_1 << ")" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1098 : : }
1099 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
1100 : : {
1101 : 0 : EntityHandle v = conn4[i];
1102 [ # # ]: 0 : int index = mesh_verts.index( v );
1103 [ # # ]: 0 : assert( -1 != index );
1104 [ # # ]: 0 : TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1105 : : }
1106 [ # # ]: 0 : for( int k = num_nodes; k < max_edges_1; k++ )
1107 : : {
1108 : 0 : TLq.vi_wr[sizeTuple * n + 2 + k] =
1109 : 0 : 0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1110 : : }
1111 : 0 : int currentIndexIntTuple = 2 + max_edges_1;
1112 : : // is the mesh migrated before or not?
1113 [ # # ]: 0 : if( migrated_mesh )
1114 : : {
1115 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( orgSendProcTag, &q, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sender for polygon, in migrate scenario" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1116 : 0 : TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender; // should be different than -1
1117 : 0 : currentIndexIntTuple++;
1118 : : }
1119 : : // GLOBAL_DOFS info, if available
1120 [ # # ]: 0 : if( size_gdofs_tag )
1121 : : {
1122 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gdsTag, &q, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't get gdofs data in HOMME" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1123 [ # # ]: 0 : for( int i = 0; i < size_gdofs_tag; i++ )
1124 : : {
1125 : 0 : TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
1126 [ # # ]: 0 : valsDOFs[i]; // should be different than 0 or -1
1127 : : }
1128 : : }
1129 : :
1130 [ # # ]: 0 : TLq.inc_n(); // increment tuple list size
1131 : : }
1132 : 0 : } // end for loop over total number of processors
1133 : :
1134 : : // now we are done populating the tuples; route them to the appropriate processors
1135 : : // this does the communication magic
1136 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
[ # # ]
1137 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
[ # # ]
1138 : :
1139 : : // the first mesh elements are in localEnts; we do not need them at all
1140 : :
1141 : : // maps from global ids to new vertex and cell handles, that are added
1142 : :
1143 [ # # ]: 0 : std::map< int, EntityHandle > globalID_to_vertex_handle;
1144 : : // we already have some vertices from second mesh set; they are already in the processor, even
1145 : : // before receiving other verts from neighbors this is an inverse map from gid to vertex handle,
1146 : : // which is local here, we do not want to duplicate vertices their identifier is the global ID!!
1147 : : // it must be unique per mesh ! (I mean, second mesh); gid gor first mesh is not needed here
1148 : 0 : int k = 0;
1149 [ # # ][ # # ]: 0 : for( Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
[ # # ][ # # ]
[ # # ]
1150 : : {
1151 [ # # ][ # # ]: 0 : globalID_to_vertex_handle[gids[k]] = *vit;
[ # # ]
1152 : : }
1153 : : /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one?
1154 : 0 : globalID_to_eh.clear(); // we do not really need it, but we keep it for debugging mostly
1155 : :
1156 : : // now, look at every TLv, and see if we have to create a vertex there or not
1157 [ # # ]: 0 : int n = TLv.get_n(); // the size of the points received
1158 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1159 : : {
1160 : 0 : int globalId = TLv.vi_rd[2 * i + 1];
1161 [ # # ][ # # ]: 0 : if( globalID_to_vertex_handle.find( globalId ) ==
[ # # ]
1162 : : globalID_to_vertex_handle.end() ) // we do not have locally this vertex (yet)
1163 : : // so we have to create it, and add to the inverse map
1164 : : {
1165 : : EntityHandle new_vert;
1166 : 0 : double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1167 [ # # ][ # # ]: 0 : rval = mb->create_vertex( dp_pos, new_vert );MB_CHK_SET_ERR( rval, "can't create new vertex " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1168 [ # # ]: 0 : globalID_to_vertex_handle[globalId] = new_vert; // now add it to the map
1169 : : // set the GLOBAL ID tag on the new vertex
1170 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_vert, 1, &globalId );MB_CHK_SET_ERR( rval, "can't set global ID tag on new vertex " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1171 : : }
1172 : : }
1173 : :
1174 : : // now, all necessary vertices should be created
1175 : : // look in the local list of 2d cells for this proc, and add all those cells to covering set
1176 : : // also
1177 : :
1178 [ # # ]: 0 : Range& local = Rto[my_rank];
1179 [ # # ]: 0 : Range local_q = local.subset_by_dimension( 2 );
1180 : :
1181 [ # # ][ # # ]: 0 : for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
[ # # ][ # # ]
[ # # ]
1182 : : {
1183 [ # # ]: 0 : EntityHandle q = *it; // these are from lagr cells, local
1184 : : int gid_el;
1185 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &gid_el );MB_CHK_SET_ERR( rval, "can't get global id of cell " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1186 [ # # ]: 0 : assert( gid_el >= 0 );
1187 [ # # ]: 0 : globalID_to_eh[gid_el] = q; // do we need this? yes, now we do; parent tags are now using it heavily
1188 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( sendProcTag, &q, 1, &my_rank );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1189 : : }
1190 : :
1191 : : // now look at all elements received through; we do not want to duplicate them
1192 [ # # ]: 0 : n = TLq.get_n(); // number of elements received by this processor
1193 : : // a cell should be received from one proc only; so why are we so worried about duplicated
1194 : : // elements? a vertex can be received from multiple sources, that is fine
1195 : :
1196 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1197 : : {
1198 : 0 : int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1199 : : // int from_proc=TLq.vi_rd[sizeTuple * i ]; // we do not need from_proc anymore
1200 : :
1201 : : // do we already have a quad with this global ID, represented? no way !
1202 : : // if (globalID_to_eh.find(globalIdEl) == globalID_to_eh.end())
1203 : : //{
1204 : : // construct the conn triangle , quad or polygon
1205 : : EntityHandle new_conn[MAXEDGES]; // we should use std::vector with max_edges_1
1206 : 0 : int nnodes = -1;
1207 [ # # ]: 0 : for( int j = 0; j < max_edges_1; j++ )
1208 : : {
1209 : 0 : int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1210 [ # # ]: 0 : if( vgid == 0 )
1211 : 0 : new_conn[j] = 0; // this can actually happen for polygon mesh (when we have less
1212 : : // number of vertices than max_edges)
1213 : : else
1214 : : {
1215 [ # # ][ # # ]: 0 : assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
[ # # ]
1216 [ # # ]: 0 : new_conn[j] = globalID_to_vertex_handle[vgid];
1217 : 0 : nnodes = j + 1; // nodes are at the beginning, and are variable number
1218 : : }
1219 : : }
1220 : : EntityHandle new_element;
1221 : : //
1222 : 0 : EntityType entType = MBQUAD;
1223 [ # # ]: 0 : if( nnodes > 4 ) entType = MBPOLYGON;
1224 [ # # ]: 0 : if( nnodes < 4 ) entType = MBTRI;
1225 [ # # ][ # # ]: 0 : rval = mb->create_element( entType, new_conn, nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element for second mesh " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1226 : :
1227 [ # # ]: 0 : globalID_to_eh[globalIdEl] = new_element;
1228 [ # # ]: 0 : local_q.insert( new_element );
1229 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set gid for cell " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1230 : 0 : int currentIndexIntTuple = 2 + max_edges_1;
1231 [ # # ]: 0 : if( migrated_mesh )
1232 : : {
1233 : 0 : orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1234 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( orgSendProcTag, &new_element, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't set original sender for cell, in migrate scenario" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1235 : 0 : currentIndexIntTuple++; // add one more
1236 : : }
1237 : : // check if we need to retrieve and set GLOBAL_DOFS data
1238 [ # # ]: 0 : if( size_gdofs_tag )
1239 : : {
1240 [ # # ]: 0 : for( int j = 0; j < size_gdofs_tag; j++ )
1241 : : {
1242 [ # # ]: 0 : valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
1243 : : }
1244 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gdsTag, &new_element, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't set GLOBAL_DOFS data on coverage mesh" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1245 : : }
1246 : : // store also the processor this coverage element came from
1247 : 0 : int from_proc = TLq.vi_rd[sizeTuple * i];
1248 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( sendProcTag, &new_element, 1, &from_proc );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1249 : : }
1250 : :
1251 : : // now, create a new set, covering_set
1252 [ # # ][ # # ]: 0 : rval = mb->add_entities( covering_set, local_q );MB_CHK_SET_ERR( rval, "can't add entities to new mesh set " );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1253 : : #ifdef VERBOSE
1254 : : std::cout << " proc " << my_rank << " add " << local_q.size() << " cells to covering set \n";
1255 : : #endif
1256 : 0 : return MB_SUCCESS;
1257 : : }
1258 : :
1259 : : #endif // MOAB_HAVE_MPI
1260 : : //#undef VERBOSE
1261 : : #undef CHECK_CONVEXITY
1262 : :
1263 [ + - ][ + - ]: 4 : } /* namespace moab */
|
