1/*
2 * Intx2MeshOnSphere.cpp
3 *
4 * Created on: Oct 3, 2012
5 */67#ifdef _MSC_VER /* windows */8#define _USE_MATH_DEFINES // For M_PI9#endif1011#include"moab/IntxMesh/Intx2MeshOnSphere.hpp"12#include"moab/IntxMesh/IntxUtils.hpp"13#include"moab/GeomUtil.hpp"14#include"moab/BoundBox.hpp"15#include"moab/MeshTopoUtil.hpp"16#ifdef MOAB_HAVE_MPI17#include"moab/ParallelComm.hpp"18#endif19#include"MBTagConventions.hpp"2021#include<cassert>2223// #define ENABLE_DEBUG24#define CHECK_CONVEXITY25namespace moab
26 {
2728 Intx2MeshOnSphere::Intx2MeshOnSphere( Interface* mbimpl, IntxAreaUtils::AreaMethod amethod )
29 : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
30 {
31 }
3233 Intx2MeshOnSphere::~Intx2MeshOnSphere() {}
3435/*
36 * return also the area for robustness verification
37 */38doubleIntx2MeshOnSphere::setup_tgt_cell( EntityHandle tgt, int& nsTgt )
39 {
4041// get coordinates of the target quad, to decide the gnomonic plane42double cellArea = 0;
4344int num_nodes;
45 ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
4647 nsTgt = num_nodes;
48// account for possible padded polygons49while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
50 nsTgt--;
5152// CartVect coords[4];53 rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
5455 CartVect middle = tgtCoords[0];
56for( int i = 1; i < nsTgt; i++ )
57 middle += tgtCoords[i];
58 middle = 1. / nsTgt * middle;
5960 IntxUtils::decide_gnomonic_plane( middle, plane ); // output the plane61for( int j = 0; j < nsTgt; j++ )
62 {
63// populate coords in the plane for intersection64// they should be oriented correctly, positively65 rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
66 }
6768for( int j = 1; j < nsTgt - 1; j++ )
69 cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
7071// take target coords in order and compute area in plane72return cellArea;
73 }
7475/* the elements are convex for sure, then do a gnomonic projection of both,
76 * compute intersection in the plane, then go back to the sphere for the points
77 * */78ErrorCode Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt,
79 EntityHandle src,
80double* P,
81int& nP,
82double& area,
83int markb[MAXEDGES],
84int markr[MAXEDGES],
85int& nsBlue,
86int& nsTgt,
87bool check_boxes_first )
88 {
89// the area will be used from now on, to see how well we fill the target cell with polygons90// the points will be at most 40; they will describe a convex patch, after the points will be91// ordered and collapsed (eliminate doubles)9293// CartVect srccoords[4];94int num_nodes = 0;
95 ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
96 nsBlue = num_nodes;
97// account for possible padded polygons98while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
99 nsBlue--;
100 rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
101102 area = 0.;
103 nP = 0; // number of intersection points we are marking the boundary of src!104if( check_boxes_first )
105 {
106// look at the boxes formed with vertices; if they are far away, return false early107// make sure the target is setup already108setup_tgt_cell( tgt, nsTgt ); // we do not need area here109// use here gnomonic plane (plane) to see where source is110bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
111int planeb;
112 CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
113 IntxUtils::decide_gnomonic_plane( mid3, planeb );
114if( !overlap3d && ( plane != planeb ) ) // plane was set at setup_tgt_cell115return MB_SUCCESS; // no error, but no intersection, decide early to get out116// if same plane, still check for gnomonic plane in 2d117// if no overlap in 2d, get out118if( !overlap3d && plane == planeb ) // CHECK 2D too119 {
120for( int j = 0; j < nsBlue; j++ )
121 {
122 rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j],
123 srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
124 }
125bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
126if( !overlap2d ) return MB_SUCCESS; // we are sure they are not overlapping in 2d , either127 }
128 }
129#ifdef ENABLE_DEBUG130if( dbg_1 )
131 {
132 std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
133for( int j = 0; j < nsTgt; j++ )
134 {
135 std::cout << tgtCoords[j] << "\n";
136 }
137 std::cout << "src " << mb->id_from_handle( src ) << "\n";
138for( int j = 0; j < nsBlue; j++ )
139 {
140 std::cout << srcCoords[j] << "\n";
141 }
142 mb->list_entities( &tgt, 1 );
143 mb->list_entities( &src, 1 );
144 }
145#endif146147for( int j = 0; j < nsBlue; j++ )
148 {
149 rval = IntxUtils::gnomonic_projection( srcCoords[j], Rsrc, plane, srcCoords2D[2 * j], srcCoords2D[2 * j + 1] );MB_CHK_ERR( rval );
150 }
151152#ifdef ENABLE_DEBUG153if( dbg_1 )
154 {
155 std::cout << "gnomonic plane: " << plane << "\n";
156 std::cout << " target src\n";
157for( int j = 0; j < nsTgt; j++ )
158 {
159 std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
160 }
161for( int j = 0; j < nsBlue; j++ )
162 {
163 std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
164 }
165 }
166#endif167168 rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
169170int side[MAXEDGES] = { 0 }; // this refers to what side? source or tgt?171int extraPoints =
172 IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
173if( extraPoints >= 1 )
174 {
175for( int k = 0; k < nsBlue; k++ )
176 {
177if( side[k] )
178 {
179// this means that vertex k of source is inside convex tgt; mark edges k-1 and k in180// src,181// as being "intersected" by tgt; (even though they might not be intersected by182// other edges, the fact that their apex is inside, is good enough)183 markb[k] = 1;
184 markb[( k + nsBlue - 1 ) % nsBlue] =
1851; // it is the previous edge, actually, but instead of doing -1, it is186// better to do modulo +3 (modulo 4)187// null side b for next call188 side[k] = 0;
189 }
190 }
191 }
192 nP += extraPoints;
193194 extraPoints =
195 IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
196if( extraPoints >= 1 )
197 {
198for( int k = 0; k < nsTgt; k++ )
199 {
200if( side[k] )
201 {
202// this is to mark that target edges k-1 and k are intersecting src203 markr[k] = 1;
204 markr[( k + nsTgt - 1 ) % nsTgt] =
2051; // it is the previous edge, actually, but instead of doing -1, it is206// better to do modulo +3 (modulo 4)207// null side b for next call208 }
209 }
210 }
211 nP += extraPoints;
212213// now sort and orient the points in P, such that they are forming a convex polygon214// this will be the foundation of our new mesh215// this works if the polygons are convex216 IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 ); // nP should be at most 8 in the end ?217// if there are more than 3 points, some area will be positive218219if( nP >= 3 )
220 {
221for( int k = 1; k < nP - 1; k++ )
222 area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
223#ifdef CHECK_CONVEXITY224// each edge should be large enough that we can compute angles between edges225for( int k = 0; k < nP; k++ )
226 {
227int k1 = ( k + 1 ) % nP;
228int k2 = ( k1 + 1 ) % nP;
229double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
230if( orientedArea < 0 )
231 {
232 std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
233 << " " << src << " \n";
234 }
235 }
236#endif237 }
238239return MB_SUCCESS; // no error240 }
241242// this method will also construct the triangles/quads/polygons in the new mesh243// if we accept planar polygons, we just save them244// also, we could just create new vertices every time, and merge only in the end;245// could be too expensive, and the tolerance for merging could be an246// interesting topic247ErrorCode Intx2MeshOnSphere::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsBlue, double* iP, int nP )
248 {
249#ifdef ENABLE_DEBUG250// first of all, check against target and source vertices251//252if( dbg_1 )
253 {
254 std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
255 << "\n";
256for( int n = 0; n < nP; n++ )
257 std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
258 }
259#endif260261// get the edges for the target triangle; the extra points will be on those edges, saved as262// lists (unordered)263264// first get the list of edges adjacent to the target cell265// use the neighTgtEdgeTag266 EntityHandle adjTgtEdges[MAXEDGES];
267 ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge target tag" );
268// we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small269// potatoes some of them will be handles to the initial vertices from source or target meshes270271 std::vector< EntityHandle > foundIds;
272 foundIds.resize( nP );
273#ifdef CHECK_CONVEXITY274int npBefore1 = nP;
275int oldNodes = 0;
276int otherIntx = 0;
277 moab::IntxAreaUtils areaAdaptor;
278#endif279for( int i = 0; i < nP; i++ )
280 {
281double* pp = &iP[2 * i]; // iP+2*i282// project the point back on the sphere283 CartVect pos;
284 IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
285int found = 0;
286// first, are they on vertices from target or src?287// priority is the target mesh (mb2?)288int j = 0;
289 EntityHandle outNode = (EntityHandle)0;
290for( j = 0; j < nsTgt && !found; j++ )
291 {
292// int node = tgtTri.v[j];293double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
294if( d2 < epsilon_1 / 1000 ) // two orders of magnitude smaller than it should, to avoid concave polygons295 {
296297 foundIds[i] = tgtConn[j]; // no new node298 found = 1;
299#ifdef CHECK_CONVEXITY300 oldNodes++;
301#endif302#ifdef ENABLE_DEBUG303if( dbg_1 )
304 std::cout << " target node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
305 << " 2d coords:" << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << " d2: " << d2
306 << " \n";
307#endif308 }
309 }
310311for( j = 0; j < nsBlue && !found; j++ )
312 {
313// int node = srcTri.v[j];314double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
315if( d2 < epsilon_1 / 1000 )
316 {
317// suspect is srcConn[j] corresponding in mbOut318319 foundIds[i] = srcConn[j]; // no new node320 found = 1;
321#ifdef CHECK_CONVEXITY322 oldNodes++;
323#endif324#ifdef ENABLE_DEBUG325if( dbg_1 )
326 std::cout << " source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
327 << " \n";
328#endif329 }
330 }
331332if( !found )
333 {
334// find the edge it belongs, first, on the red element335// look at the minimum area, not at the first below some tolerance336double minArea = 1.e+38;
337int index_min = -1;
338for( j = 0; j < nsTgt; j++ )
339 {
340int j1 = ( j + 1 ) % nsTgt;
341double area = fabs( IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp ) );
342// how to check if pp is between redCoords2D[j] and redCoords2D[j1] ?343// they should form a straight line; the sign should be -1344double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
345 IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
346 IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
347if( area < minArea && checkx < 2 * epsilon_1 ) // round off error or not?348 {
349 index_min = j;
350 minArea = area;
351 }
352 }
353// verify that index_min is valid354assert( index_min >= 0 );
355356if( minArea < epsilon_1 / 2 ) // we found the smallest area, so we think we found the357// target edge it belongs358 {
359// found the edge; now find if there is a point in the list here360// std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];361int indx = TgtEdges.index( adjTgtEdges[index_min] );
362if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)363 {
364 std::cerr << " error in adjacent target edge: " << mb->id_from_handle( adjTgtEdges[index_min] )
365 << "\n";
366return MB_FAILURE;
367 }
368 std::vector< EntityHandle >* expts = extraNodesVec[indx];
369// if the points pp is between extra points, then just give that id370// if not, create a new point, (check the id)371// get the coordinates of the extra points so far372int nbExtraNodesSoFar = expts->size();
373if( nbExtraNodesSoFar > 0 )
374 {
375 std::vector< CartVect > coords1;
376 coords1.resize( nbExtraNodesSoFar );
377 mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
378// std::list<int>::iterator it;379for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
380 {
381// int pnt = *it;382double d2 = ( pos - coords1[k] ).length();
383if( d2 < 2 * epsilon_1 ) // is this below machine precision?384 {
385 found = 1;
386 foundIds[i] = ( *expts )[k];
387#ifdef CHECK_CONVEXITY388 otherIntx++;
389#endif390 }
391 }
392 }
393if( !found )
394 {
395// create a new point in 2d (at the intersection)396// foundIds[i] = m_num2dPoints;397// expts.push_back(m_num2dPoints);398// need to create a new node in mbOut399// this will be on the edge, and it will be added to the local list400 rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
401 ( *expts ).push_back( outNode );
402// CID 181168; avoid leak storage error403 rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
404 foundIds[i] = outNode;
405 found = 1;
406 }
407 }
408 }
409if( !found )
410 {
411 std::cout << " target quad: ";
412for( int j1 = 0; j1 < nsTgt; j1++ )
413 {
414 std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
415 }
416 std::cout << " a point pp is not on a target quad " << *pp << " " << pp[1] << " target quad "417 << mb->id_from_handle( tgt ) << " \n";
418return MB_FAILURE;
419 }
420 }
421#ifdef ENABLE_DEBUG422if( dbg_1 )
423 {
424 std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
425for( int i1 = 0; i1 < nP; i1++ )
426 std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
427 }
428#endif429// first, find out if we have nodes collapsed; shrink them430// we may have to reduce nP431// it is possible that some nodes are collapsed after intersection only432// nodes will always be in order (convex intersection)433#ifdef CHECK_CONVEXITY434int npBefore2 = nP;
435#endif436correct_polygon( &foundIds[0], nP );
437// now we can build the triangles, from P array, with foundIds438// we will put them in the out set439if( nP >= 3 )
440 {
441 EntityHandle polyNew;
442 rval = mb->create_element( MBPOLYGON, &foundIds[0], nP, polyNew );MB_CHK_ERR( rval );
443 rval = mb->add_entities( outSet, &polyNew, 1 );MB_CHK_ERR( rval );
444445// tag it with the global ids from target and source elements446int globalID;
447 rval = mb->tag_get_data( gid, &src, 1, &globalID );MB_CHK_ERR( rval );
448 rval = mb->tag_set_data( srcParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
449// if(!parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<450// " : Blue = " << globalID << ", " << mb->id_from_handle(src) << "\t\n";451 rval = mb->tag_get_data( gid, &tgt, 1, &globalID );MB_CHK_ERR( rval );
452 rval = mb->tag_set_data( tgtParentTag, &polyNew, 1, &globalID );MB_CHK_ERR( rval );
453// if(parcomm->rank()) std::cout << "Setting parent for " << mb->id_from_handle(polyNew) <<454// " : target = " << globalID << ", " << mb->id_from_handle(tgt) << "\n";455456 counting++;
457 rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
458if( orgSendProcTag )
459 {
460int org_proc = -1;
461 rval = mb->tag_get_data( orgSendProcTag, &src, 1, &org_proc );MB_CHK_ERR( rval );
462 rval = mb->tag_set_data( orgSendProcTag, &polyNew, 1, &org_proc );MB_CHK_ERR( rval ); // yet another tag463 }
464#ifdef CHECK_CONVEXITY465// each edge should be large enough that we can compute angles between edges466 std::vector< double > coords;
467 coords.resize( 3 * nP );
468 rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
469std::vector< CartVect > posi( nP );
470 rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
471472for( int k = 0; k < nP; k++ )
473 {
474int k1 = ( k + 1 ) % nP;
475int k2 = ( k1 + 1 ) % nP;
476double orientedArea =
477 areaAdaptor.area_spherical_triangle( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
478if( orientedArea < 0 )
479 {
480 std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
481for( int i = 0; i < nP; i++ )
482 {
483int nexti = ( i + 1 ) % nP;
484double lengthEdge = ( posi[i] - posi[nexti] ).length();
485 std::cout << " " << foundIds[i] << " edge en:" << lengthEdge << "\n";
486 }
487 std::cout << " old verts: " << oldNodes << " other intx:" << otherIntx << "\n";
488489 std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
490 << " target, src:" << tgt << " " << src << " \n";
491 }
492 }
493#endif494495#ifdef ENABLE_DEBUG496if( dbg_1 )
497 {
498 std::cout << "Counting: " << counting << "\n";
499 std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
500for( int i1 = 0; i1 < nP; i1++ )
501 std::cout << " " << mb->id_from_handle( foundIds[i1] );
502 std::cout << " plane: " << plane << "\n";
503std::vector< CartVect > posi( nP );
504 mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
505for( int i1 = 0; i1 < nP; i1++ )
506 std::cout << foundIds[i1] << " " << posi[i1] << "\n";
507508 std::stringstream fff;
509 fff << "file0" << counting << ".vtk";
510 rval = mb->write_mesh( fff.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
511 }
512#endif513 }
514// else {515// std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";516// }517// disable_debug();518return MB_SUCCESS;
519 }
520521ErrorCode Intx2MeshOnSphere::update_tracer_data( EntityHandle out_set, Tag& tagElem, Tag& tagArea )
522 {
523 EntityHandle dum = 0;
524525 Tag corrTag;
526 ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE,
527 &dum ); // it should have been created528MB_CHK_SET_ERR( rval, "can't get correlation tag" );
529530// get all polygons out of out_set; then see where are they coming from531 Range polys;
532 rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
533534// rs2 is the target range, arrival; rs1 is src, departure;535// there is a connection between rs1 and rs2, through the corrTag536// corrTag is __correlation537// basically, mb->tag_get_data(corrTag, &(tgtPoly), 1, &srcPoly);538// also, mb->tag_get_data(corrTag, &(srcPoly), 1, &tgtPoly);539// we start from rs2 existing, then we have to update something540541// tagElem will have multiple tracers542int numTracers = 0;
543 rval = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
544if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
545546std::vector< double > currentVals( rs2.size() * numTracers );
547 rval = mb->tag_get_data( tagElem, rs2, ¤tVals[0] );MB_CHK_SET_ERR( rval, "can't get existing tracers values" );
548549// create new tuple list for tracers to other processors, from remote_cells550#ifdef MOAB_HAVE_MPI551if( remote_cells )
552 {
553int n = remote_cells->get_n();
554if( n > 0 )
555 {
556 remote_cells_with_tracers = newTupleList();
557 remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
558 n ); // tracers are in these tuples559 remote_cells_with_tracers->enableWriteAccess();
560for( int i = 0; i < n; i++ )
561 {
562 remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i];
563 remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
564// remote_cells->vr_wr[i] = 0.; will have a different tuple for communication565 remote_cells_with_tracers->vul_wr[i] =
566 remote_cells->vul_wr[i]; // this is the corresponding target cell (arrival)567for( int k = 0; k < numTracers; k++ )
568 remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0; // initialize tracers to be transported569 remote_cells_with_tracers->inc_n();
570 }
571 }
572delete remote_cells;
573 remote_cells = NULL;
574 }
575#endif576// for each polygon, we have 2 indices: target and source parents577// we need index source to update index tgt?578std::vector< double > newValues( rs2.size() * numTracers,
5790. ); // initialize with 0 all of them580// area of the polygon * conc on target (old) current quantity581// finally, divide by the area of the tgt582double check_intx_area = 0.;
583moab::IntxAreaUtils intxAreas( this->areaMethod ); // use_lHuiller = true584for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
585 {
586 EntityHandle poly = *it;
587int srcIndex, tgtIndex;
588 rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
589590 EntityHandle src = rs1[srcIndex - 1]; // big assumption, it should work for meshes where global id is the same591// as element handle (ordered from 1 to number of elements); should be OK for Homme meshes592 rval = mb->tag_get_data( tgtParentTag, &poly, 1, &tgtIndex );MB_CHK_SET_ERR( rval, "can't get target tag" );
593// EntityHandle target = rs2[tgtIndex];594// big assumption here, target and source are "parallel" ;we should have an index from595// source to target (so a deformed source corresponds to an arrival tgt)596/// TODO: VSM: Its unclear whether we need the source or destination radius here.597double radius = Rsrc;
598double areap = intxAreas.area_spherical_element( mb, poly, radius );
599 check_intx_area += areap;
600// so the departure cell at time t (srcIndex) covers a portion of a tgtCell601// that quantity will be transported to the tgtCell at time t+dt602// the source corresponds to a target arrival603 EntityHandle tgtArr;
604 rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
605if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
606 {
607#ifdef MOAB_HAVE_MPI608if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
609// maybe the element is remote, from another processor610int global_id_src;
611 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" );
612// find the613int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
614if( index_in_remote == -1 )
615MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
616for( int k = 0; k < numTracers; k++ )
617 remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
618 currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
619#endif620 }
621elseif( MB_SUCCESS == rval )
622 {
623int arrTgtIndex = rs2.index( tgtArr );
624if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
625for( int k = 0; k < numTracers; k++ )
626 newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
627 }
628629else630MB_CHK_SET_ERR( rval, "can't get arrival target for corresponding " );
631 }
632// now, send back the remote_cells_with_tracers to the processors they came from, with the633// updated values for the tracer mass in a cell634#ifdef MOAB_HAVE_MPI635if( remote_cells_with_tracers )
636 {
637// so this means that some cells will be sent back with tracer info to the procs they were638// sent from639 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
640// now, look at the global id, find the proper "tgt" cell with that index and update its641// mass642// remote_cells->print("remote cells after routing");643int n = remote_cells_with_tracers->get_n();
644for( int j = 0; j < n; j++ )
645 {
646 EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j]; // entity handle sent back647int arrTgtIndex = rs2.index( tgtCell );
648if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
649for( int k = 0; k < numTracers; k++ )
650 newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
651 }
652 }
653#endif/* MOAB_HAVE_MPI */654// now divide by target area (current)655int j = 0;
656 Range::iterator iter = rs2.begin();
657void* data = NULL; // used for stored area658int count = 0;
659std::vector< double > total_mass_local( numTracers, 0. );
660while( iter != rs2.end() )
661 {
662 rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
663double* ptrArea = (double*)data;
664for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
665 {
666for( int k = 0; k < numTracers; k++ )
667 {
668 total_mass_local[k] += newValues[j * numTracers + k];
669 newValues[j * numTracers + k] /= ( *ptrArea );
670 }
671 }
672 }
673 rval = mb->tag_set_data( tagElem, rs2, &newValues[0] );MB_CHK_SET_ERR( rval, "can't set new values tag" );
674675#ifdef MOAB_HAVE_MPI676std::vector< double > total_mass( numTracers, 0. );
677double total_intx_area = 0;
678int mpi_err =
679MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
680if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
681// now reduce total area682 mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
683if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
684if( my_rank == 0 )
685 {
686for( int k = 0; k < numTracers; k++ )
687 std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass[k] << "\n";
688 std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "689 << total_intx_area << "\n";
690 }
691692if( remote_cells_with_tracers )
693 {
694delete remote_cells_with_tracers;
695 remote_cells_with_tracers = NULL;
696 }
697#else698for( int k = 0; k < numTracers; k++ )
699 std::cout << "total mass now tracer k=" << k + 1 << " " << total_mass_local[k] << "\n";
700 std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * Rsrc * Rsrc << " "701 << check_intx_area << "\n";
702#endif703return MB_SUCCESS;
704 }
705706#ifdef MOAB_HAVE_MPI707ErrorCode Intx2MeshOnSphere::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts, bool gnomonic )
708 {
709if( !gnomonic )
710 {
711return Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts, gnomonic );
712 }
713// so here, we know that the logic is for gnomonic == true714 localEnts.clear();
715 ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );MB_CHK_SET_ERR( rval, "can't get local ents" );
716717 rval = mb->get_connectivity( localEnts, local_verts );MB_CHK_SET_ERR( rval, "can't get connectivity" );
718int num_local_verts = (int)local_verts.size();
719720assert( parcomm != NULL );
721722if( num_local_verts == 0 )
723 {
724// it is probably point cloud, get the local vertices from set725 rval = mb->get_entities_by_dimension( euler_set, 0, local_verts );MB_CHK_SET_ERR( rval, "can't get local vertices from set" );
726 num_local_verts = (int)local_verts.size();
727 localEnts = local_verts;
728 }
729// will use 6 gnomonic planes to decide boxes for each gnomonic plane730// each gnomonic box will be 2d, min, max731double gnom_box[24];
732for( int i = 0; i < 6; i++ )
733 {
734 gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
735 gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
736 }
737738// there are 6 gnomonic planes; some elements could be on the corners, and affect multiple739// planes decide what gnomonic planes will be affected by each cell some elements could appear740// in multiple gnomonic planes !741std::vector< double > coords( 3 * num_local_verts );
742 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 " );
743// decide each local vertex to what gnomonic plane it belongs744745 std::vector< int > gnplane;
746 gnplane.resize( num_local_verts );
747for( int i = 0; i < num_local_verts; i++ )
748 {
749CartVect pos( &coords[3 * i] );
750int pl;
751 IntxUtils::decide_gnomonic_plane( pos, pl );
752 gnplane[i] = pl;
753 }
754755for( Range::iterator it = localEnts.begin(); it != localEnts.end(); it++ )
756 {
757 EntityHandle cell = *it;
758 EntityType typeCell = mb->type_from_handle( cell ); // could be vertex, for point cloud759// get coordinates, and decide gnomonic planes for it760int nnodes;
761const EntityHandle* conn = NULL;
762 EntityHandle c[1];
763if( typeCell != MBVERTEX )
764 {
765 rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
766 }
767else768 {
769 nnodes = 1;
770 c[0] = cell; // actual node771 conn = &c[0];
772 }
773// get coordinates of vertices involved with this774std::vector< double > elco( 3 * nnodes );
775 std::set< int > planes;
776for( int i = 0; i < nnodes; i++ )
777 {
778int ix = local_verts.index( conn[i] );
779 planes.insert( gnplane[ix] );
780for( int j = 0; j < 3; j++ )
781 {
782 elco[3 * i + j] = coords[3 * ix + j];
783 }
784 }
785// now, augment the boxes for all planes involved786for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
787 {
788int pl = *st;
789for( int i = 0; i < nnodes; i++ )
790 {
791CartVect pos( &elco[3 * i] );
792double c2[2];
793 IntxUtils::gnomonic_projection( pos, Rdest, pl, c2[0],
794 c2[1] ); // 2 coordinates795//796for( int k = 0; k < 2; k++ )
797 {
798double val = c2[k];
799if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val; // min in k direction800if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
801 gnom_box[4 * ( pl - 1 ) + 2 + k] = val; // max in k direction802 }
803 }
804 }
805 }
806807int numprocs = parcomm->proc_config().proc_size();
808 allBoxes.resize( 24 * numprocs ); // 6 gnomonic planes , 4 doubles for each for 2d box809810 my_rank = parcomm->proc_config().proc_rank();
811for( int k = 0; k < 24; k++ )
812 allBoxes[24 * my_rank + k] = gnom_box[k];
813814// now communicate to get all boxes815int mpi_err;
816#if( MPI_VERSION >= 2 )817// use "in place" option818 mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 24, MPI_DOUBLE,
819 parcomm->proc_config().proc_comm() );
820#else821 {
822std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
823 mpi_err = MPI_Allgather( &allBoxes[24 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
824 parcomm->proc_config().proc_comm() );
825 allBoxes = allBoxes_tmp;
826 }
827#endif828if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
829830#ifdef VERBOSE831if( my_rank == 0 )
832 {
833 std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
834 << " on second mesh \n";
835for( int i = 0; i < numprocs; i++ )
836 {
837 std::cout << "task: " << i << " \n";
838for( int pl = 1; pl <= 6; pl++ )
839 {
840 std::cout << " plane " << pl << " min: \t" << allBoxes[24 * i + 4 * ( pl - 1 )] << " \t"841 << allBoxes[24 * i + 4 * ( pl - 1 ) + 1] << "\n";
842 std::cout << " \t max: \t" << allBoxes[24 * i + 4 * ( pl - 1 ) + 2] << " \t"843 << allBoxes[24 * i + 4 * ( pl - 1 ) + 3] << "\n";
844 }
845 }
846 }
847#endif848849return MB_SUCCESS;
850 }
851//#define VERBOSE852// this will use the bounding boxes for the (euler)/ fix mesh that are already established853// will distribute the mesh to other procs, so that on each task, the covering set covers the local854// bounding box this means it will cover the second (local) mesh set; So the covering set will cover855// completely the second local mesh set (in intersection)856// now, when covering set needs to have extra layers, we will increase dramatically the box_eps, from something close to 0,857// to something larger than the "source mesh size" * sqrt(3) for each layer needed858// so the first step is finding the global diagonal mesh size in the source mesh859ErrorCode Intx2MeshOnSphere::construct_covering_set( EntityHandle& initial_distributed_set,
860 EntityHandle& covering_set,
861bool gnomonic,
862int nb_ghost_layers )
863 {
864// primary element came from, in the joint communicator ; this will be forwarded by coverage865// mesh needed for tag migrate later on866int defaultInt = -1; // no processor, so it was not migrated from somewhere else867 ErrorCode rval = mb->tag_get_handle( "orig_sending_processor", 1, MB_TYPE_INTEGER, orgSendProcTag,
868 MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create original sending processor tag" );
869870assert( parcomm != NULL );
871 Range meshCells;
872 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" );
873874if( 1 == parcomm->proc_config().proc_size() )
875 {
876// move all initial cells to coverage set877 rval = mb->add_entities( covering_set, meshCells );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
878// if point cloud source, add vertices879if( 0 == meshCells.size() || max_edges_1 == 0 )
880 {
881// add vertices from the source set882 Range verts;
883 rval = mb->get_entities_by_dimension( initial_distributed_set, 0, verts );MB_CHK_SET_ERR( rval, "can't get vertices from mesh set" );
884 rval = mb->add_entities( covering_set, verts );MB_CHK_SET_ERR( rval, "can't add primary ents to covering set" );
885 }
886return MB_SUCCESS;
887 }
888889// mark on the coverage mesh where this element came from890 Tag sendProcTag; /// for coverage mesh, will store the sender891 rval = mb->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag, MB_TAG_DENSE | MB_TAG_CREAT,
892 &defaultInt );MB_CHK_SET_ERR( rval, "can't create sending processor tag" );
893894// this information needs to be forwarded to coverage mesh, if this mesh was already migrated895// from somewhere else896// look at the value of orgSendProcTag for one mesh cell; if -1, no need to forward that; if897// !=-1, we know that this mesh was migrated, we need to find out more about origin of cell898int orig_sender = -1;
899 EntityHandle oneCell = 0;
900// decide if we need to transfer global DOFs info attached to each HOMME coarse cell; first we901// need to decide if the mesh has that tag; will affect the size of the tuple list involved in902// the crystal routing903int size_gdofs_tag = 0;
904 std::vector< int > valsDOFs;
905 Tag gdsTag;
906 rval = mb->tag_get_handle( "GLOBAL_DOFS", gdsTag );
907908if( meshCells.size() > 0 )
909 {
910 oneCell = meshCells[0]; // it is possible we do not have any cells, even after migration911 rval = mb->tag_get_data( orgSendProcTag, &oneCell, 1, &orig_sender );MB_CHK_SET_ERR( rval, "can't get original sending processor value" );
912if( gdsTag )
913 {
914 DataType dtype;
915 rval = mb->tag_get_data_type( gdsTag, dtype );
916if( MB_SUCCESS == rval && MB_TYPE_INTEGER == dtype )
917 {
918// find the values on first cell919int lenTag = 0;
920 rval = mb->tag_get_length( gdsTag, lenTag );
921if( MB_SUCCESS == rval && lenTag > 0 )
922 {
923 valsDOFs.resize( lenTag );
924 rval = mb->tag_get_data( gdsTag, &oneCell, 1, &valsDOFs[0] );
925if( MB_SUCCESS == rval && valsDOFs[0] > 0 )
926 {
927// first value positive means we really need to transport this data during928// coverage929 size_gdofs_tag = lenTag;
930 }
931 }
932 }
933 }
934 }
935936// another collective call, to see if the mesh is migrated and if the GLOBAL_DOFS tag need to be937// transferred over to the coverage mesh. It is possible that there is no initial mesh source938// mesh on the task, so we do not know that info from the tag but TupleList needs to be sized939// uniformly for all tasks. Do a collective MPI_MAX to see if it is migrated and if we have940// (collectively) a GLOBAL_DOFS task941942int local_int_array[2], global_int_array[2];
943 local_int_array[0] = orig_sender;
944 local_int_array[1] = size_gdofs_tag;
945// now reduce over all processors946int mpi_err =
947MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
948if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
949 orig_sender = global_int_array[0];
950 size_gdofs_tag = global_int_array[1];
951#ifdef VERBOSE952 std::cout << "proc: " << my_rank << " size_gdofs_tag:" << size_gdofs_tag << "\n";
953#endif954 valsDOFs.resize( size_gdofs_tag );
955956// finally, we have correct global info, we can decide if mesh was migrated and if we have957// global dofs tag that need to be sent with coverage info958int migrated_mesh = 0;
959if( orig_sender != -1 ) migrated_mesh = 1; //960961// if size_gdofs_tag>0, we are sure valsDOFs got resized to what we need962963// get all mesh verts1964 Range mesh_verts;
965 rval = mb->get_connectivity( meshCells, mesh_verts );MB_CHK_SET_ERR( rval, "can't get mesh vertices" );
966size_t num_mesh_verts = mesh_verts.size();
967968// now see the mesh points positions; to what boxes should we send them?969std::vector< double > coords_mesh( 3 * num_mesh_verts );
970 rval = mb->get_coords( mesh_verts, &coords_mesh[0] );MB_CHK_SET_ERR( rval, "can't get mesh points position" );
971972// decide gnomonic plane for each vertex, as in the compute boxes973 std::vector< int > gnplane;
974if( gnomonic )
975 {
976 gnplane.resize( num_mesh_verts );
977for( size_t i = 0; i < num_mesh_verts; i++ )
978 {
979CartVect pos( &coords_mesh[3 * i] );
980int pl;
981 IntxUtils::decide_gnomonic_plane( pos, pl );
982 gnplane[i] = pl;
983 }
984 }
985986std::vector< int > gids( num_mesh_verts );
987 rval = mb->tag_get_data( gid, mesh_verts, &gids[0] );MB_CHK_SET_ERR( rval, "can't get vertices gids" );
988989// ranges to send to each processor; will hold vertices and elements (quads/ polygons)990// will look if the box of the mesh cell covers bounding box(es) (within tolerances)991 std::map< int, Range > Rto;
992int numprocs = parcomm->proc_config().proc_size();
993994// now, box error is pretty small, in general995// for bilinear mesh, we need an extra layer, which we will get by increasing the epsilon to catch the extra layer996// it will depend on the size of the source mesh997// so we will compute the max diagonal length for each cell, on the sphere, so we will modify box_error998if( nb_ghost_layers > 0 )
999 {
1000double diagonal;
1001 rval = IntxUtils::max_diagonal( mb, meshCells, max_edges_1, diagonal );MB_CHK_SET_ERR( rval, "can't get max diagonal" );
1002//1003double global_diag = 0;
1004 mpi_err = MPI_Allreduce( &diagonal, &global_diag, 1, MPI_DOUBLE, MPI_MAX, parcomm->proc_config().proc_comm() );
1005if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
1006double extra_thickness = global_diag * nb_ghost_layers;
1007if( gnomonic ) extra_thickness *= sqrt( 3. );
1008 box_error += extra_thickness; //1009if( !my_rank )
1010 std::cout << "ghost_layers:" << nb_ghost_layers << " max diagonal:" << global_diag
1011 << " extra thickness:" << extra_thickness << " box_error:" << box_error << "\n";
1012 }
1013for( Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
1014 {
1015 EntityHandle q = *eit;
1016const EntityHandle* conn;
1017int num_nodes;
1018 rval = mb->get_connectivity( q, conn, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity on cell" );
10191020// first decide what planes need to consider1021 std::set< int > planes; // if this list contains more than 3 planes, we have a very bad mesh!!!1022std::vector< double > elco( 3 * num_nodes );
1023for( int i = 0; i < num_nodes; i++ )
1024 {
1025 EntityHandle v = conn[i];
1026int index = mesh_verts.index( v );
1027if( gnomonic ) planes.insert( gnplane[index] );
1028for( int j = 0; j < 3; j++ )
1029 {
1030 elco[3 * i + j] = coords_mesh[3 * index + j]; // extract from coords1031 }
1032 }
1033if( gnomonic )
1034 {
1035// now loop over all planes that need to be considered for this element1036for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
1037 {
1038int pl = *st; // gnomonic plane considered1039double qmin[2] = { DBL_MAX, DBL_MAX };
1040double qmax[2] = { -DBL_MAX, -DBL_MAX };
1041for( int i = 0; i < num_nodes; i++ )
1042 {
1043CartVect dp( &elco[3 * i] ); // uses constructor for CartVect that takes a1044// pointer to double1045// gnomonic projection1046double c2[2];
1047 IntxUtils::gnomonic_projection( dp, Rsrc, pl, c2[0], c2[1] ); // 2 coordinates1048for( int j = 0; j < 2; j++ )
1049 {
1050if( qmin[j] > c2[j] ) qmin[j] = c2[j];
1051if( qmax[j] < c2[j] ) qmax[j] = c2[j];
1052 }
1053 }
1054// now decide if processor p should be interested in this cell, by looking at plane pl1055// 2d box this is one of the few size n loops;1056for( int p = 0; p < numprocs; p++ ) // each cell q can be sent to more than one processor1057 {
1058double procMin1 = allBoxes[24 * p + 4 * ( pl - 1 )]; // these were determined before1059//1060if( procMin1 >= DBL_MAX ) // the processor has no targets on this plane1061continue;
1062double procMin2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
1063double procMax1 = allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
1064double procMax2 = allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
1065// test overlap of 2d boxes1066if( procMin1 > qmax[0] + box_error || procMin2 > qmax[1] + box_error ) continue; //1067if( qmin[0] > procMax1 + box_error || qmin[1] > procMax2 + box_error ) continue;
1068// good to be inserted1069 Rto[p].insert( q );
1070 }
1071 }
1072 }
1073else// regular 3d box; one box per processor1074 {
1075for( int p = 0; p < numprocs; p++ )
1076 {
1077 BoundBox box( &allBoxes[6 * p] );
1078bool insert = false;
1079for( int i = 0; i < num_nodes; i++ )
1080 {
1081if( box.contains_point( &elco[3 * i], box_error ) )
1082 {
1083 insert = true;
1084break;
1085 }
1086 }
1087if( insert ) Rto[p].insert( q );
1088 }
1089 }
1090 }
10911092// here, we will use crystal router to send each cell to designated tasks (mesh migration)10931094// a better implementation would be to use pcomm send / recv entities; a good test case1095// pcomm send / receives uses point to point communication, not global gather / scatter10961097// now, build TLv and TLq (tuple list for vertices and cells, separately sent)1098size_t numq = 0;
1099size_t numv = 0;
11001101// merge the list of vertices to be sent1102for( int p = 0; p < numprocs; p++ )
1103 {
1104if( p == (int)my_rank ) continue; // do not "send" it to current task, because it is already here1105 Range& range_to_P = Rto[p];
1106// add the vertices to it1107if( range_to_P.empty() ) continue; // nothing to send to proc p1108#ifdef VERBOSE1109 std::cout << " proc : " << my_rank << " to proc " << p << " send " << range_to_P.size() << " cells "1110 << " psize: " << range_to_P.psize() << "\n";
1111#endif1112 Range vertsToP;
1113 rval = mb->get_connectivity( range_to_P, vertsToP );MB_CHK_SET_ERR( rval, "can't get connectivity" );
1114 numq = numq + range_to_P.size();
1115 numv = numv + vertsToP.size();
1116 range_to_P.merge( vertsToP );
1117 }
11181119 TupleList TLv; // send vertices with a different tuple list1120 TupleList TLq;
1121 TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, 3 real coordinates1122 TLv.enableWriteAccess();
11231124// add also GLOBAL_DOFS info, if found on the mesh cell; it should be found only on HOMME cells!1125int sizeTuple =
11262 + max_edges_1 + migrated_mesh + size_gdofs_tag; // max edges could be up to MAXEDGES :) for polygons1127 TLq.initialize( sizeTuple, 0, 0, 0,
1128 numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v), plus1129// original sender if set (migrated mesh case)1130// we will not send the entity handle, global ID should be more than enough1131// we will not need more than 2B vertices TODO 2B vertices or cells1132// if we need more than 2B, we will need to use a different marker anyway (GLOBAL ID is not1133// enough then)11341135 TLq.enableWriteAccess();
1136#ifdef VERBOSE1137 std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1138#endif11391140for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1141 {
1142if( to_proc == (int)my_rank ) continue;
1143 Range& range_to_P = Rto[to_proc];
1144 Range V = range_to_P.subset_by_type( MBVERTEX );
11451146for( Range::iterator it = V.begin(); it != V.end(); ++it )
1147 {
1148 EntityHandle v = *it;
1149int index = mesh_verts.index( v ); //1150assert( -1 != index );
1151int n = TLv.get_n(); // current size of tuple list1152 TLv.vi_wr[2 * n] = to_proc; // send to processor1153 TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the second_mesh_verts range1154 TLv.vr_wr[3 * n] = coords_mesh[3 * index]; // departure position, of the node local_verts[i]1155 TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
1156 TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
1157 TLv.inc_n(); // increment tuple list size1158 }
1159// also, prep the 2d cells for sending ...1160 Range Q = range_to_P.subset_by_dimension( 2 );
1161for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
1162 {
1163 EntityHandle q = *it; // this is a second mesh cell (or src, lagrange set)1164int global_id;
1165 rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_SET_ERR( rval, "can't get gid for polygon" );
1166int n = TLq.get_n(); // current size1167 TLq.vi_wr[sizeTuple * n] = to_proc; //1168 TLq.vi_wr[sizeTuple * n + 1] =
1169 global_id; // global id of element, used to identify it for debug purposes only1170const EntityHandle* conn4;
1171int num_nodes; // could be up to MAXEDGES; max_edges?;1172 rval = mb->get_connectivity( q, conn4, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity for cell" );
1173if( num_nodes > max_edges_1 )
1174 {
1175 mb->list_entities( &q, 1 );
1176MB_CHK_SET_ERR( MB_FAILURE, "too many nodes in a cell (" << num_nodes << "," << max_edges_1 << ")" );
1177 }
1178for( int i = 0; i < num_nodes; i++ )
1179 {
1180 EntityHandle v = conn4[i];
1181int index = mesh_verts.index( v );
1182assert( -1 != index );
1183 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1184 }
1185for( int k = num_nodes; k < max_edges_1; k++ )
1186 {
1187 TLq.vi_wr[sizeTuple * n + 2 + k] =
11880; // fill the rest of node ids with 0; we know that the node ids start from 1!1189 }
1190int currentIndexIntTuple = 2 + max_edges_1;
1191// is the mesh migrated before or not?1192if( migrated_mesh )
1193 {
1194// case of extra work, maybe need to check if it is ghost ? yes, next loop !1195 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" );
1196 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender; // should be different than -11197 currentIndexIntTuple++;
1198 }
1199// GLOBAL_DOFS info, if available1200if( size_gdofs_tag )
1201 {
1202 rval = mb->tag_get_data( gdsTag, &q, 1, &valsDOFs[0] );MB_CHK_SET_ERR( rval, "can't get gdofs data in HOMME" );
1203for( int i = 0; i < size_gdofs_tag; i++ )
1204 {
1205 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
1206 valsDOFs[i]; // should be different than 0 or -11207 }
1208 }
12091210 TLq.inc_n(); // increment tuple list size1211 }
1212 } // end for loop over total number of processors12131214// now we are done populating the tuples; route them to the appropriate processors1215// this does the communication magic1216 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1217 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
12181219// the first mesh elements are in localEnts; we do not need them at all12201221// maps from global ids to new vertex and cell handles, that are added12221223 std::map< int, EntityHandle > globalID_to_vertex_handle;
1224// we already have some vertices from second mesh set; they are already in the processor, even1225// before receiving other verts from neighbors this is an inverse map from gid to vertex handle,1226// which is local here, we do not want to duplicate vertices their identifier is the global ID!!1227// it must be unique per mesh ! (I mean, first mesh, source); gid for second mesh is not needed here1228int k = 0;
1229for( Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
1230 {
1231 globalID_to_vertex_handle[gids[k]] = *vit;
1232 }
1233/*std::map<int, EntityHandle> globalID_to_eh;*/// do we need this one?1234 globalID_to_eh.clear(); // we need it now in case of extra work, to not duplicate cells12351236// now, look at every TLv, and see if we have to create a vertex there or not1237int n = TLv.get_n(); // the size of the points received1238for( int i = 0; i < n; i++ )
1239 {
1240int globalId = TLv.vi_rd[2 * i + 1];
1241if( globalID_to_vertex_handle.find( globalId ) ==
1242 globalID_to_vertex_handle.end() ) // we do not have locally this vertex (yet)1243// so we have to create it, and add to the inverse map1244 {
1245 EntityHandle new_vert;
1246double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1247 rval = mb->create_vertex( dp_pos, new_vert );MB_CHK_SET_ERR( rval, "can't create new vertex " );
1248 globalID_to_vertex_handle[globalId] = new_vert; // now add it to the map1249// set the GLOBAL ID tag on the new vertex1250 rval = mb->tag_set_data( gid, &new_vert, 1, &globalId );MB_CHK_SET_ERR( rval, "can't set global ID tag on new vertex " );
1251 }
1252 }
12531254// now, all necessary vertices should be created1255// look in the local list of 2d cells for this proc, and add all those cells to covering set1256// also12571258 Range& local = Rto[my_rank];
1259 Range local_q = local.subset_by_dimension( 2 );
12601261for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1262 {
1263 EntityHandle q = *it; // these are from source cells, local1264int gid_el;
1265 rval = mb->tag_get_data( gid, &q, 1, &gid_el );MB_CHK_SET_ERR( rval, "can't get global id of cell " );
1266assert( gid_el >= 0 );
1267 globalID_to_eh[gid_el] = q; // do we need this? yes, now we do; parent tags are now using it heavily1268 rval = mb->tag_set_data( sendProcTag, &q, 1, &my_rank );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
1269 }
12701271// now look at all elements received through; we do not want to duplicate them1272 n = TLq.get_n(); // number of elements received by this processor1273// a cell should be received from one proc only; so why are we so worried about duplicated1274// elements? a vertex can be received from multiple sources, that is fine12751276for( int i = 0; i < n; i++ )
1277 {
1278int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1279// int from_proc=TLq.vi_rd[sizeTuple * i ]; // we do not need from_proc anymore12801281// do we already have a cell with this global ID, represented?1282// yes, it could happen for extraWork !1283if( globalID_to_eh.find( globalIdEl ) != globalID_to_eh.end() ) continue;
1284// construct the conn triangle , quad or polygon1285 EntityHandle new_conn[MAXEDGES]; // we should use std::vector with max_edges_11286int nnodes = -1;
1287for( int j = 0; j < max_edges_1; j++ )
1288 {
1289int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID1290if( vgid == 0 )
1291 new_conn[j] = 0; // this can actually happen for polygon mesh (when we have less1292// number of vertices than max_edges)1293else1294 {
1295assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
1296 new_conn[j] = globalID_to_vertex_handle[vgid];
1297 nnodes = j + 1; // nodes are at the beginning, and are variable number1298 }
1299 }
1300 EntityHandle new_element;
1301//1302 EntityType entType = MBQUAD;
1303if( nnodes > 4 ) entType = MBPOLYGON;
1304if( nnodes < 4 ) entType = MBTRI;
1305 rval = mb->create_element( entType, new_conn, nnodes, new_element );MB_CHK_SET_ERR( rval, "can't create new element for second mesh " );
13061307 globalID_to_eh[globalIdEl] = new_element;
1308 local_q.insert( new_element );
1309 rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );MB_CHK_SET_ERR( rval, "can't set gid for cell " );
1310int currentIndexIntTuple = 2 + max_edges_1;
1311if( migrated_mesh )
1312 {
1313 orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1314 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" );
1315 currentIndexIntTuple++; // add one more1316 }
1317// store also the processor this coverage element came from1318int from_proc = TLq.vi_rd[sizeTuple * i];
1319 rval = mb->tag_set_data( sendProcTag, &new_element, 1, &from_proc );MB_CHK_SET_ERR( rval, "can't set sender for cell" );
13201321// check if we need to retrieve and set GLOBAL_DOFS data1322if( size_gdofs_tag )
1323 {
1324for( int j = 0; j < size_gdofs_tag; j++ )
1325 {
1326 valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
1327 }
1328 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" );
1329 }
1330 }
13311332// now, add to the covering_set the elements created in the local_q range1333 rval = mb->add_entities( covering_set, local_q );MB_CHK_SET_ERR( rval, "can't add entities to new mesh set " );
1334#ifdef VERBOSE1335 std::cout << " proc " << my_rank << " add " << local_q.size() << " cells to covering set \n";
1336#endif1337return MB_SUCCESS;
1338 }
13391340#endif// MOAB_HAVE_MPI1341//#undef VERBOSE1342#undef CHECK_CONVEXITY13431344 } /* namespace moab */