 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#include <Intx2MeshOnSphere.hpp>
Inheritance diagram for moab::Intx2MeshOnSphere: Collaboration diagram for moab::Intx2MeshOnSphere:Public Attributes | |
| const IntxAreaUtils::AreaMethod | areaMethod |
Private Attributes | |
| int | plane |
| double | Rsrc |
| double | Rdest |
Additional Inherited Members | |
| Protected Attributes inherited from moab::Intx2Mesh | |
| Interface * | mb |
| EntityHandle | mbs1 |
| EntityHandle | mbs2 |
| Range | rs1 |
| Range | rs2 |
| EntityHandle | outSet |
| Tag | gid |
| Tag | TgtFlagTag |
| Range | TgtEdges |
| Tag | tgtParentTag |
| Tag | srcParentTag |
| Tag | countTag |
| Tag | srcNeighTag |
| Tag | tgtNeighTag |
| Tag | neighTgtEdgeTag |
| Tag | orgSendProcTag |
| Tag | imaskTag |
| for coverage mesh, will store the original sender More... | |
| const EntityHandle * | tgtConn |
| const EntityHandle * | srcConn |
| CartVect | tgtCoords [MAXEDGES] |
| CartVect | srcCoords [MAXEDGES] |
| double | tgtCoords2D [MAXEDGES2] |
| double | srcCoords2D [MAXEDGES2] |
| std::vector< std::vector< EntityHandle > * > | extraNodesVec |
| double | epsilon_1 |
| double | epsilon_area |
| std::vector< double > | allBoxes |
| double | box_error |
| EntityHandle | localRoot |
| Range | localEnts |
| unsigned int | my_rank |
| int | max_edges_1 |
| int | max_edges_2 |
| int | counting |
Definition at line 16 of file Intx2MeshOnSphere.hpp.
| moab::Intx2MeshOnSphere::Intx2MeshOnSphere | ( | Interface * | mbimpl, |
| IntxAreaUtils::AreaMethod | amethod = IntxAreaUtils::lHuiller |
||
| ) |
Definition at line 28 of file Intx2MeshOnSphere.cpp.
29 : Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 ) 30 { 31 }
|
virtual |
Definition at line 33 of file Intx2MeshOnSphere.cpp.
33 {}
|
virtual |
Implements moab::Intx2Mesh.
Definition at line 78 of file Intx2MeshOnSphere.cpp.
88 {
89 // the area will be used from now on, to see how well we fill the target cell with polygons
90 // the points will be at most 40; they will describe a convex patch, after the points will be
91 // ordered and collapsed (eliminate doubles)
92
93 // CartVect srccoords[4];
94 int 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 polygons
98 while( srcConn[nsBlue - 2] == srcConn[nsBlue - 1] && nsBlue > 3 )
99 nsBlue--;
100 rval = mb->get_coords( srcConn, nsBlue, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
101
102 area = 0.;
103 nP = 0; // number of intersection points we are marking the boundary of src!
104 if( check_boxes_first )
105 {
106 // look at the boxes formed with vertices; if they are far away, return false early
107 // make sure the target is setup already
108 setup_tgt_cell( tgt, nsTgt ); // we do not need area here
109 // use here gnomonic plane (plane) to see where source is
110 bool overlap3d = GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsBlue, box_error );
111 int planeb;
112 CartVect mid3 = ( srcCoords[0] + srcCoords[1] + srcCoords[2] ) / 3;
113 IntxUtils::decide_gnomonic_plane( mid3, planeb );
114 if( !overlap3d && ( plane != planeb ) ) // plane was set at setup_tgt_cell
115 return MB_SUCCESS; // no error, but no intersection, decide early to get out
116 // if same plane, still check for gnomonic plane in 2d
117 // if no overlap in 2d, get out
118 if( !overlap3d && plane == planeb ) // CHECK 2D too
119 {
120 for( 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 }
125 bool overlap2d = GeomUtil::bounding_boxes_overlap_2d( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, box_error );
126 if( !overlap2d ) return MB_SUCCESS; // we are sure they are not overlapping in 2d , either
127 }
128 }
129 #ifdef ENABLE_DEBUG
130 if( dbg_1 )
131 {
132 std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
133 for( int j = 0; j < nsTgt; j++ )
134 {
135 std::cout << tgtCoords[j] << "\n";
136 }
137 std::cout << "src " << mb->id_from_handle( src ) << "\n";
138 for( 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 #endif
146
147 for( 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 }
151
152 #ifdef ENABLE_DEBUG
153 if( dbg_1 )
154 {
155 std::cout << "gnomonic plane: " << plane << "\n";
156 std::cout << " target src\n";
157 for( int j = 0; j < nsTgt; j++ )
158 {
159 std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n";
160 }
161 for( int j = 0; j < nsBlue; j++ )
162 {
163 std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
164 }
165 }
166 #endif
167
168 rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
169
170 int side[MAXEDGES] = { 0 }; // this refers to what side? source or tgt?
171 int extraPoints =
172 IntxUtils::borderPointsOfXinY2( srcCoords2D, nsBlue, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
173 if( extraPoints >= 1 )
174 {
175 for( int k = 0; k < nsBlue; k++ )
176 {
177 if( side[k] )
178 {
179 // this means that vertex k of source is inside convex tgt; mark edges k-1 and k in
180 // src,
181 // as being "intersected" by tgt; (even though they might not be intersected by
182 // other edges, the fact that their apex is inside, is good enough)
183 markb[k] = 1;
184 markb[( k + nsBlue - 1 ) % nsBlue] =
185 1; // it is the previous edge, actually, but instead of doing -1, it is
186 // better to do modulo +3 (modulo 4)
187 // null side b for next call
188 side[k] = 0;
189 }
190 }
191 }
192 nP += extraPoints;
193
194 extraPoints =
195 IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsBlue, &( P[2 * nP] ), side, epsilon_area );
196 if( extraPoints >= 1 )
197 {
198 for( int k = 0; k < nsTgt; k++ )
199 {
200 if( side[k] )
201 {
202 // this is to mark that target edges k-1 and k are intersecting src
203 markr[k] = 1;
204 markr[( k + nsTgt - 1 ) % nsTgt] =
205 1; // it is the previous edge, actually, but instead of doing -1, it is
206 // better to do modulo +3 (modulo 4)
207 // null side b for next call
208 }
209 }
210 }
211 nP += extraPoints;
212
213 // now sort and orient the points in P, such that they are forming a convex polygon
214 // this will be the foundation of our new mesh
215 // this works if the polygons are convex
216 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 positive
218
219 if( nP >= 3 )
220 {
221 for( int k = 1; k < nP - 1; k++ )
222 area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
223 #ifdef CHECK_CONVEXITY
224 // each edge should be large enough that we can compute angles between edges
225 for( int k = 0; k < nP; k++ )
226 {
227 int k1 = ( k + 1 ) % nP;
228 int k2 = ( k1 + 1 ) % nP;
229 double orientedArea = IntxUtils::area2D( &P[2 * k], &P[2 * k1], &P[2 * k2] );
230 if( orientedArea < 0 )
231 {
232 std::cout << " oriented area is negative: " << orientedArea << " k:" << k << " target, src:" << tgt
233 << " " << src << " \n";
234 }
235 }
236 #endif
237 }
238
239 return MB_SUCCESS; // no error
240 }
References moab::IntxUtils::area2D(), moab::IntxUtils::borderPointsOfXinY2(), moab::GeomUtil::bounding_boxes_overlap(), moab::GeomUtil::bounding_boxes_overlap_2d(), moab::Intx2Mesh::box_error, moab::IntxUtils::decide_gnomonic_plane(), moab::IntxUtils::EdgeIntersections2(), moab::Intx2Mesh::epsilon_1, moab::Intx2Mesh::epsilon_area, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::IntxUtils::gnomonic_projection(), moab::Interface::id_from_handle(), moab::Interface::list_entities(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_ERR, MB_SUCCESS, plane, Rsrc, setup_tgt_cell(), moab::IntxUtils::SortAndRemoveDoubles2(), moab::Intx2Mesh::srcConn, moab::Intx2Mesh::srcCoords, moab::Intx2Mesh::srcCoords2D, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.
|
virtual |
Implements moab::Intx2Mesh.
Definition at line 247 of file Intx2MeshOnSphere.cpp.
248 {
249 #ifdef ENABLE_DEBUG
250 // first of all, check against target and source vertices
251 //
252 if( dbg_1 )
253 {
254 std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
255 << "\n";
256 for( int n = 0; n < nP; n++ )
257 std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
258 }
259 #endif
260
261 // get the edges for the target triangle; the extra points will be on those edges, saved as
262 // lists (unordered)
263
264 // first get the list of edges adjacent to the target cell
265 // use the neighTgtEdgeTag
266 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 small
269 // potatoes some of them will be handles to the initial vertices from source or target meshes
270
271 std::vector< EntityHandle > foundIds;
272 foundIds.resize( nP );
273 #ifdef CHECK_CONVEXITY
274 int npBefore1 = nP;
275 int oldNodes = 0;
276 int otherIntx = 0;
277 moab::IntxAreaUtils areaAdaptor;
278 #endif
279 for( int i = 0; i < nP; i++ )
280 {
281 double* pp = &iP[2 * i]; // iP+2*i
282 // project the point back on the sphere
283 CartVect pos;
284 IntxUtils::reverse_gnomonic_projection( pp[0], pp[1], Rdest, plane, pos );
285 int found = 0;
286 // first, are they on vertices from target or src?
287 // priority is the target mesh (mb2?)
288 int j = 0;
289 EntityHandle outNode = (EntityHandle)0;
290 for( j = 0; j < nsTgt && !found; j++ )
291 {
292 // int node = tgtTri.v[j];
293 double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
294 if( d2 < epsilon_1 / 1000 ) // two orders of magnitude smaller than it should, to avoid concave polygons
295 {
296
297 foundIds[i] = tgtConn[j]; // no new node
298 found = 1;
299 #ifdef CHECK_CONVEXITY
300 oldNodes++;
301 #endif
302 #ifdef ENABLE_DEBUG
303 if( 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 #endif
308 }
309 }
310
311 for( j = 0; j < nsBlue && !found; j++ )
312 {
313 // int node = srcTri.v[j];
314 double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
315 if( d2 < epsilon_1 / 1000 )
316 {
317 // suspect is srcConn[j] corresponding in mbOut
318
319 foundIds[i] = srcConn[j]; // no new node
320 found = 1;
321 #ifdef CHECK_CONVEXITY
322 oldNodes++;
323 #endif
324 #ifdef ENABLE_DEBUG
325 if( dbg_1 )
326 std::cout << " source node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2
327 << " \n";
328 #endif
329 }
330 }
331
332 if( !found )
333 {
334 // find the edge it belongs, first, on the red element
335 // look at the minimum area, not at the first below some tolerance
336 double minArea = 1.e+38;
337 int index_min = -1;
338 for( j = 0; j < nsTgt; j++ )
339 {
340 int j1 = ( j + 1 ) % nsTgt;
341 double 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 -1
344 double checkx = IntxUtils::dist2( &tgtCoords2D[2 * j], pp ) +
345 IntxUtils::dist2( &tgtCoords2D[2 * j1], pp ) -
346 IntxUtils::dist2( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1] );
347 if( 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 valid
354 assert( index_min >= 0 );
355
356 if( minArea < epsilon_1 / 2 ) // we found the smallest area, so we think we found the
357 // target edge it belongs
358 {
359 // found the edge; now find if there is a point in the list here
360 // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
361 int indx = TgtEdges.index( adjTgtEdges[index_min] );
362 if( 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";
366 return MB_FAILURE;
367 }
368 std::vector< EntityHandle >* expts = extraNodesVec[indx];
369 // if the points pp is between extra points, then just give that id
370 // if not, create a new point, (check the id)
371 // get the coordinates of the extra points so far
372 int nbExtraNodesSoFar = expts->size();
373 if( 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;
379 for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
380 {
381 // int pnt = *it;
382 double d2 = ( pos - coords1[k] ).length();
383 if( d2 < 2 * epsilon_1 ) // is this below machine precision?
384 {
385 found = 1;
386 foundIds[i] = ( *expts )[k];
387 #ifdef CHECK_CONVEXITY
388 otherIntx++;
389 #endif
390 }
391 }
392 }
393 if( !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 mbOut
399 // this will be on the edge, and it will be added to the local list
400 rval = mb->create_vertex( pos.array(), outNode );MB_CHK_ERR( rval );
401 ( *expts ).push_back( outNode );
402 // CID 181168; avoid leak storage error
403 rval = mb->add_entities( outSet, &outNode, 1 );MB_CHK_ERR( rval );
404 foundIds[i] = outNode;
405 found = 1;
406 }
407 }
408 }
409 if( !found )
410 {
411 std::cout << " target quad: ";
412 for( 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";
418 return MB_FAILURE;
419 }
420 }
421 #ifdef ENABLE_DEBUG
422 if( dbg_1 )
423 {
424 std::cout << " candidate polygon: nP" << nP << " plane: " << plane << "\n";
425 for( int i1 = 0; i1 < nP; i1++ )
426 std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
427 }
428 #endif
429 // first, find out if we have nodes collapsed; shrink them
430 // we may have to reduce nP
431 // it is possible that some nodes are collapsed after intersection only
432 // nodes will always be in order (convex intersection)
433 #ifdef CHECK_CONVEXITY
434 int npBefore2 = nP;
435 #endif
436 correct_polygon( &foundIds[0], nP );
437 // now we can build the triangles, from P array, with foundIds
438 // we will put them in the out set
439 if( 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 );
444
445 // tag it with the global ids from target and source elements
446 int 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";
455
456 counting++;
457 rval = mb->tag_set_data( countTag, &polyNew, 1, &counting );MB_CHK_ERR( rval );
458 if( orgSendProcTag )
459 {
460 int 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 tag
463 }
464 #ifdef CHECK_CONVEXITY
465 // each edge should be large enough that we can compute angles between edges
466 std::vector< double > coords;
467 coords.resize( 3 * nP );
468 rval = mb->get_coords( &foundIds[0], nP, &coords[0] );MB_CHK_ERR( rval );
469 std::vector< CartVect > posi( nP );
470 rval = mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );MB_CHK_ERR( rval );
471
472 for( int k = 0; k < nP; k++ )
473 {
474 int k1 = ( k + 1 ) % nP;
475 int k2 = ( k1 + 1 ) % nP;
476 double orientedArea =
477 areaAdaptor.area_spherical_triangle( &coords[3 * k], &coords[3 * k1], &coords[3 * k2], Rdest );
478 if( orientedArea < 0 )
479 {
480 std::cout << " np before 1 , 2, current " << npBefore1 << " " << npBefore2 << " " << nP << "\n";
481 for( int i = 0; i < nP; i++ )
482 {
483 int nexti = ( i + 1 ) % nP;
484 double 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";
488
489 std::cout << "rank:" << my_rank << " oriented area in 3d is negative: " << orientedArea << " k:" << k
490 << " target, src:" << tgt << " " << src << " \n";
491 }
492 }
493 #endif
494
495 #ifdef ENABLE_DEBUG
496 if( dbg_1 )
497 {
498 std::cout << "Counting: " << counting << "\n";
499 std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
500 for( int i1 = 0; i1 < nP; i1++ )
501 std::cout << " " << mb->id_from_handle( foundIds[i1] );
502 std::cout << " plane: " << plane << "\n";
503 std::vector< CartVect > posi( nP );
504 mb->get_coords( &foundIds[0], nP, &( posi[0][0] ) );
505 for( int i1 = 0; i1 < nP; i1++ )
506 std::cout << foundIds[i1] << " " << posi[i1] << "\n";
507
508 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 #endif
513 }
514 // else {
515 // std::cout << "[[FAILURE]] Number of vertices in polygon is less than 3\n";
516 // }
517 // disable_debug();
518 return MB_SUCCESS;
519 }
References moab::Interface::add_entities(), moab::IntxUtils::area2D(), moab::IntxAreaUtils::area_spherical_triangle(), moab::CartVect::array(), moab::Intx2Mesh::correct_polygon(), moab::Intx2Mesh::counting, moab::Intx2Mesh::countTag, moab::Interface::create_element(), moab::Interface::create_vertex(), moab::IntxUtils::dist2(), moab::Intx2Mesh::epsilon_1, ErrorCode, moab::Intx2Mesh::extraNodesVec, moab::Interface::get_coords(), moab::Intx2Mesh::gid, moab::Interface::id_from_handle(), moab::Range::index(), length(), MAXEDGES, moab::Intx2Mesh::mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, moab::Intx2Mesh::my_rank, moab::Intx2Mesh::neighTgtEdgeTag, moab::Intx2Mesh::orgSendProcTag, moab::Intx2Mesh::outSet, plane, Rdest, moab::IntxUtils::reverse_gnomonic_projection(), moab::Intx2Mesh::srcConn, moab::Intx2Mesh::srcCoords2D, moab::Intx2Mesh::srcParentTag, moab::Interface::tag_get_data(), moab::Interface::tag_set_data(), moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords2D, moab::Intx2Mesh::TgtEdges, moab::Intx2Mesh::tgtParentTag, and moab::Interface::write_mesh().
|
inline |
Definition at line 27 of file Intx2MeshOnSphere.hpp.
28 { 29 Rdest = radius; 30 }
References Rdest.
Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().
|
inline |
Definition at line 23 of file Intx2MeshOnSphere.hpp.
24 { 25 Rsrc = radius; 26 }
References Rsrc.
Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().
|
virtual |
Implements moab::Intx2Mesh.
Definition at line 38 of file Intx2MeshOnSphere.cpp.
39 {
40
41 // get coordinates of the target quad, to decide the gnomonic plane
42 double cellArea = 0;
43
44 int num_nodes;
45 ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );MB_CHK_ERR_RET_VAL( rval, cellArea );
46
47 nsTgt = num_nodes;
48 // account for possible padded polygons
49 while( tgtConn[nsTgt - 2] == tgtConn[nsTgt - 1] && nsTgt > 3 )
50 nsTgt--;
51
52 // CartVect coords[4];
53 rval = mb->get_coords( tgtConn, nsTgt, &( tgtCoords[0][0] ) );MB_CHK_ERR_RET_VAL( rval, cellArea );
54
55 CartVect middle = tgtCoords[0];
56 for( int i = 1; i < nsTgt; i++ )
57 middle += tgtCoords[i];
58 middle = 1. / nsTgt * middle;
59
60 IntxUtils::decide_gnomonic_plane( middle, plane ); // output the plane
61 for( int j = 0; j < nsTgt; j++ )
62 {
63 // populate coords in the plane for intersection
64 // they should be oriented correctly, positively
65 rval = IntxUtils::gnomonic_projection( tgtCoords[j], Rdest, plane, tgtCoords2D[2 * j], tgtCoords2D[2 * j + 1] );MB_CHK_ERR_RET_VAL( rval, cellArea );
66 }
67
68 for( int j = 1; j < nsTgt - 1; j++ )
69 cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
70
71 // take target coords in order and compute area in plane
72 return cellArea;
73 }
References moab::IntxUtils::area2D(), moab::IntxUtils::decide_gnomonic_plane(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::IntxUtils::gnomonic_projection(), moab::Intx2Mesh::mb, MB_CHK_ERR_RET_VAL, plane, Rdest, moab::Intx2Mesh::tgtConn, moab::Intx2Mesh::tgtCoords, and moab::Intx2Mesh::tgtCoords2D.
Referenced by computeIntersectionBetweenTgtAndSrc().
| ErrorCode moab::Intx2MeshOnSphere::update_tracer_data | ( | EntityHandle | out_set, |
| Tag & | tagElem, | ||
| Tag & | tagArea | ||
| ) |
TODO: VSM: Its unclear whether we need the source or destination radius here.
Definition at line 521 of file Intx2MeshOnSphere.cpp.
522 {
523 EntityHandle dum = 0;
524
525 Tag corrTag;
526 ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE,
527 &dum ); // it should have been created
528 MB_CHK_SET_ERR( rval, "can't get correlation tag" );
529
530 // get all polygons out of out_set; then see where are they coming from
531 Range polys;
532 rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_SET_ERR( rval, "can't get polygons out" );
533
534 // rs2 is the target range, arrival; rs1 is src, departure;
535 // there is a connection between rs1 and rs2, through the corrTag
536 // corrTag is __correlation
537 // 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 something
540
541 // tagElem will have multiple tracers
542 int numTracers = 0;
543 rval = mb->tag_get_length( tagElem, numTracers );MB_CHK_SET_ERR( rval, "can't get number of tracers in simulation" );
544 if( numTracers < 1 ) MB_CHK_SET_ERR( MB_FAILURE, "no tracers data" );
545
546 std::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" );
548
549 // create new tuple list for tracers to other processors, from remote_cells
550 #ifdef MOAB_HAVE_MPI
551 if( remote_cells )
552 {
553 int n = remote_cells->get_n();
554 if( n > 0 )
555 {
556 remote_cells_with_tracers = new TupleList();
557 remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
558 n ); // tracers are in these tuples
559 remote_cells_with_tracers->enableWriteAccess();
560 for( 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 communication
565 remote_cells_with_tracers->vul_wr[i] =
566 remote_cells->vul_wr[i]; // this is the corresponding target cell (arrival)
567 for( int k = 0; k < numTracers; k++ )
568 remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0; // initialize tracers to be transported
569 remote_cells_with_tracers->inc_n();
570 }
571 }
572 delete remote_cells;
573 remote_cells = NULL;
574 }
575 #endif
576 // for each polygon, we have 2 indices: target and source parents
577 // we need index source to update index tgt?
578 std::vector< double > newValues( rs2.size() * numTracers,
579 0. ); // initialize with 0 all of them
580 // area of the polygon * conc on target (old) current quantity
581 // finally, divide by the area of the tgt
582 double check_intx_area = 0.;
583 moab::IntxAreaUtils intxAreas( this->areaMethod ); // use_lHuiller = true
584 for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
585 {
586 EntityHandle poly = *it;
587 int srcIndex, tgtIndex;
588 rval = mb->tag_get_data( srcParentTag, &poly, 1, &srcIndex );MB_CHK_SET_ERR( rval, "can't get source tag" );
589
590 EntityHandle src = rs1[srcIndex - 1]; // big assumption, it should work for meshes where global id is the same
591 // as element handle (ordered from 1 to number of elements); should be OK for Homme meshes
592 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 from
595 // 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.
597 double radius = Rsrc;
598 double 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 tgtCell
601 // that quantity will be transported to the tgtCell at time t+dt
602 // the source corresponds to a target arrival
603 EntityHandle tgtArr;
604 rval = mb->tag_get_data( corrTag, &src, 1, &tgtArr );
605 if( 0 == tgtArr || MB_TAG_NOT_FOUND == rval )
606 {
607 #ifdef MOAB_HAVE_MPI
608 if( !remote_cells_with_tracers ) MB_CHK_SET_ERR( MB_FAILURE, "no remote cells, failure\n" );
609 // maybe the element is remote, from another processor
610 int 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 the
613 int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
614 if( index_in_remote == -1 )
615 MB_CHK_SET_ERR( MB_FAILURE, "can't find the global id element in remote cells\n" );
616 for( 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 #endif
620 }
621 else if( MB_SUCCESS == rval )
622 {
623 int arrTgtIndex = rs2.index( tgtArr );
624 if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
625 for( int k = 0; k < numTracers; k++ )
626 newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
627 }
628
629 else
630 MB_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 the
633 // updated values for the tracer mass in a cell
634 #ifdef MOAB_HAVE_MPI
635 if( remote_cells_with_tracers )
636 {
637 // so this means that some cells will be sent back with tracer info to the procs they were
638 // sent from
639 ( 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 its
641 // mass
642 // remote_cells->print("remote cells after routing");
643 int n = remote_cells_with_tracers->get_n();
644 for( int j = 0; j < n; j++ )
645 {
646 EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j]; // entity handle sent back
647 int arrTgtIndex = rs2.index( tgtCell );
648 if( -1 == arrTgtIndex ) MB_CHK_SET_ERR( MB_FAILURE, "can't find the target arrival index" );
649 for( 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)
655 int j = 0;
656 Range::iterator iter = rs2.begin();
657 void* data = NULL; // used for stored area
658 int count = 0;
659 std::vector< double > total_mass_local( numTracers, 0. );
660 while( iter != rs2.end() )
661 {
662 rval = mb->tag_iterate( tagArea, iter, rs2.end(), count, data );MB_CHK_SET_ERR( rval, "can't tag iterate" );
663 double* ptrArea = (double*)data;
664 for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
665 {
666 for( 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" );
674
675 #ifdef MOAB_HAVE_MPI
676 std::vector< double > total_mass( numTracers, 0. );
677 double total_intx_area = 0;
678 int mpi_err =
679 MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
680 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
681 // now reduce total area
682 mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
683 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
684 if( my_rank == 0 )
685 {
686 for( 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 }
691
692 if( remote_cells_with_tracers )
693 {
694 delete remote_cells_with_tracers;
695 remote_cells_with_tracers = NULL;
696 }
697 #else
698 for( 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 #endif
703 return MB_SUCCESS;
704 }
References moab::IntxAreaUtils::area_spherical_element(), areaMethod, moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), ErrorCode, moab::Interface::get_entities_by_dimension(), moab::Intx2Mesh::gid, moab::Range::index(), moab::Intx2Mesh::mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TAG_NOT_FOUND, MB_TYPE_HANDLE, moab::Intx2Mesh::my_rank, moab::Intx2Mesh::rs1, moab::Intx2Mesh::rs2, Rsrc, moab::Range::size(), moab::Intx2Mesh::srcParentTag, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_get_length(), moab::Interface::tag_iterate(), moab::Interface::tag_set_data(), and moab::Intx2Mesh::tgtParentTag.
| const IntxAreaUtils::AreaMethod moab::Intx2MeshOnSphere::areaMethod |
Definition at line 59 of file Intx2MeshOnSphere.hpp.
Referenced by update_tracer_data().
|
private |
Definition at line 62 of file Intx2MeshOnSphere.hpp.
Referenced by computeIntersectionBetweenTgtAndSrc(), findNodes(), and setup_tgt_cell().
|
private |
Definition at line 63 of file Intx2MeshOnSphere.hpp.
Referenced by findNodes(), set_radius_destination_mesh(), and setup_tgt_cell().
|
private |
Definition at line 63 of file Intx2MeshOnSphere.hpp.
Referenced by computeIntersectionBetweenTgtAndSrc(), set_radius_source_mesh(), and update_tracer_data().
1.9.1.