Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::Intx2Mesh Class Referenceabstract

#include <Intx2Mesh.hpp>

+ Inheritance diagram for moab::Intx2Mesh:
+ Collaboration diagram for moab::Intx2Mesh:

Public Member Functions

 Intx2Mesh (Interface *mbimpl)
 
virtual ~Intx2Mesh ()
 
ErrorCode intersect_meshes (EntityHandle mbs1, EntityHandle mbs2, EntityHandle &outputSet)
 
ErrorCode intersect_meshes_kdtree (EntityHandle mbset1, EntityHandle mbset2, EntityHandle &outputSet)
 
virtual ErrorCode computeIntersectionBetweenTgtAndSrc (EntityHandle tgt, EntityHandle src, double *P, int &nP, double &area, int markb[MAXEDGES], int markr[MAXEDGES], int &nsidesSrc, int &nsidesTgt, bool check_boxes_first=false)=0
 
virtual ErrorCode findNodes (EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double *iP, int nP)=0
 
virtual double setup_tgt_cell (EntityHandle tgt, int &nsTgt)=0
 
virtual ErrorCode FindMaxEdgesInSet (EntityHandle eset, int &max_edges)
 
virtual ErrorCode FindMaxEdges (EntityHandle set1, EntityHandle set2)
 
virtual ErrorCode createTags ()
 
virtual ErrorCode filterByMask (Range &cells)
 
ErrorCode DetermineOrderedNeighbors (EntityHandle inputSet, int max_edges, Tag &neighTag)
 
void set_error_tolerance (double eps)
 
void clean ()
 
void set_box_error (double berror)
 
ErrorCode create_departure_mesh_2nd_alg (EntityHandle &euler_set, EntityHandle &covering_lagr_set)
 
ErrorCode create_departure_mesh_3rd_alg (EntityHandle &lagr_set, EntityHandle &covering_set)
 
void correct_polygon (EntityHandle *foundIds, int &nP)
 

Protected Attributes

Interfacemb
 
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 EntityHandletgtConn
 
const EntityHandlesrcConn
 
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
 

Detailed Description

Definition at line 55 of file Intx2Mesh.hpp.

Constructor & Destructor Documentation

◆ Intx2Mesh()

moab::Intx2Mesh::Intx2Mesh ( Interface mbimpl)

Definition at line 28 of file Intx2Mesh.cpp.

29  : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
30  countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), imaskTag( 0 ),
31  tgtConn( NULL ), srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ),
32  my_rank( 0 )
33 #ifdef MOAB_HAVE_MPI
34  ,
35  parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
36 #endif
37  ,
38  max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
39 {
40  gid = mbimpl->globalId_tag();
41 }

References gid, and moab::Interface::globalId_tag().

◆ ~Intx2Mesh()

moab::Intx2Mesh::~Intx2Mesh ( )
virtual

Definition at line 43 of file Intx2Mesh.cpp.

44 {
45  // TODO Auto-generated destructor stub
46 #ifdef MOAB_HAVE_MPI
47  if( remote_cells )
48  {
49  delete remote_cells;
50  remote_cells = NULL;
51  }
52 #endif
53 }

Member Function Documentation

◆ clean()

void moab::Intx2Mesh::clean ( )

Definition at line 836 of file Intx2Mesh.cpp.

837 {
838  //
839  int indx = 0;
840  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
841  {
842  delete extraNodesVec[indx];
843  }
844  // extraNodesMap.clear();
845  extraNodesVec.clear();
846  // also, delete some bit tags, used to mark processed tgts and srcs
848  counting = 0; // reset counting to original value
849 }

References moab::Range::begin(), counting, moab::Range::end(), extraNodesVec, mb, moab::Interface::tag_delete(), TgtEdges, and TgtFlagTag.

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

◆ computeIntersectionBetweenTgtAndSrc()

virtual ErrorCode moab::Intx2Mesh::computeIntersectionBetweenTgtAndSrc ( EntityHandle  tgt,
EntityHandle  src,
double *  P,
int &  nP,
double &  area,
int  markb[MAXEDGES],
int  markr[MAXEDGES],
int &  nsidesSrc,
int &  nsidesTgt,
bool  check_boxes_first = false 
)
pure virtual

◆ correct_polygon()

void moab::Intx2Mesh::correct_polygon ( EntityHandle foundIds,
int &  nP 
)

Definition at line 854 of file Intx2Mesh.cpp.

855 {
856  int i = 0;
857  while( i < nP )
858  {
859  int nextIndex = ( i + 1 ) % nP;
860  if( nodes[i] == nodes[nextIndex] )
861  {
862 #ifdef ENABLE_DEBUG
863  // we need to reduce nP, and collapse nodes
864  if( dbg_1 )
865  {
866  std::cout << " nodes duplicated in list: ";
867  for( int j = 0; j < nP; j++ )
868  std::cout << nodes[j] << " ";
869  std::cout << "\n";
870  std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
871  << "\n";
872  }
873 #endif
874  // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
875  // then next thing does nothing
876  // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
877  for( int k = i; k < nP - 1; k++ )
878  nodes[k] = nodes[k + 1];
879  nP--; // decrease the number of nodes; also, decrease i, just if we may need to check
880  // again
881  i--;
882  }
883  i++;
884  }
885  return;
886 }

Referenced by moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), and moab::IntxRllCssphere::findNodes().

◆ create_departure_mesh_2nd_alg()

ErrorCode moab::Intx2Mesh::create_departure_mesh_2nd_alg ( EntityHandle euler_set,
EntityHandle covering_lagr_set 
)

◆ create_departure_mesh_3rd_alg()

ErrorCode moab::Intx2Mesh::create_departure_mesh_3rd_alg ( EntityHandle lagr_set,
EntityHandle covering_set 
)

◆ createTags()

ErrorCode moab::Intx2Mesh::createTags ( )
virtual

Definition at line 93 of file Intx2Mesh.cpp.

94 {
97  if( countTag ) mb->tag_delete( countTag );
98 
99  unsigned char def_data_bit = 0; // unused by default
100  // maybe the tgt tag is better to be deleted every time, and recreated;
101  // or is it easy to set all values to something again? like 0?
102  ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" );
103  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
104  // this is the only call that is potentially NlogN, in the whole method
105  rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
106 
107  // now, create a map from each edge to a list of potential new nodes on a tgt edge
108  // this memory has to be cleaned up
109  // change it to a vector, and use the index in range of tgt edges
110  int indx = 0;
111  extraNodesVec.resize( TgtEdges.size() );
112  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
113  {
114  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
115  extraNodesVec[indx] = nv;
116  }
117 
118  int defaultInt = -1;
119 
120  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
121  &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
122 
123  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
124  &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
125 
126  rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
127 
128  // for each cell in set 1, determine its neigh in set 1 (could be null too)
129  // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
130  rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
131  rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" );
132 
133  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
134  // them each time edges were for sure created before (tgtEdges)
135  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
136  // if we have a tag with this name, it could be of a different size, so delete it
137  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
139  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
140  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
141  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
142  {
143  EntityHandle tgtCell = *rit;
144  int num_nodes = 0;
145  rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
146  // account for padded polygons
147  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
148  num_nodes--;
149 
150  int i = 0;
151  for( i = 0; i < num_nodes; i++ )
152  {
153  EntityHandle v[2] = { tgtConn[i],
154  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
155  std::vector< EntityHandle > adj_entities;
156  rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
157  if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
158  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
159  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
160  // the tag is dense
161  }
162  // now set the value of the tag
163  rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
164  }
165  return MB_SUCCESS;
166 }

References moab::Range::begin(), countTag, DetermineOrderedNeighbors(), moab::Range::end(), ErrorCode, extraNodesVec, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::INTERSECT, max_edges_1, max_edges_2, mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_BIT, MB_TYPE_HANDLE, MB_TYPE_INTEGER, mbs1, mbs2, neighTgtEdgeTag, rs2, moab::Range::size(), srcNeighTag, srcParentTag, moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), tgtConn, TgtEdges, TgtFlagTag, tgtNeighTag, tgtParentTag, and moab::Interface::UNION.

Referenced by intersect_meshes().

◆ DetermineOrderedNeighbors()

ErrorCode moab::Intx2Mesh::DetermineOrderedNeighbors ( EntityHandle  inputSet,
int  max_edges,
Tag neighTag 
)

Definition at line 168 of file Intx2Mesh.cpp.

169 {
170  Range cells;
171  ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" );
172 
173  std::vector< EntityHandle > neighbors( max_edges );
174  std::vector< EntityHandle > zeroh( max_edges, 0 );
175  // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
176  // one on src mesh
177  rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" );
178 
179  for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
180  {
181  EntityHandle cell = *cit;
182  int nnodes = 3;
183  // will get the nnodes ordered neighbors;
184  // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
185  const EntityHandle* conn4;
186  rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
187  int nsides = nnodes;
188  // account for possible padded polygons
189  while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
190  nsides--;
191 
192  for( int i = 0; i < nsides; i++ )
193  {
194  EntityHandle v[2];
195  v[0] = conn4[i];
196  v[1] = conn4[( i + 1 ) % nsides];
197  // get all cells adjacent to these 2 vertices on the edge
198  std::vector< EntityHandle > adjcells;
199  std::vector< EntityHandle > cellsInSet;
200  rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" );
201  // now look for the cells contained in the input set;
202  // the input set should be a correct mesh, not overlapping cells, and manifold
203  size_t siz = adjcells.size();
204  for( size_t j = 0; j < siz; j++ )
205  if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
206  siz = cellsInSet.size();
207 
208  if( siz > 2 )
209  {
210  std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
211  << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" ); // non-manifold
212  }
213  if( siz == 1 )
214  {
215  // it must be the border of the input mesh;
216  neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
217  // borders do not appear for a sphere in serial, but they do appear for
218  // parallel processing anyway
219  continue;
220  }
221  // here siz ==2, it is either the first or second
222  if( cell == cellsInSet[0] )
223  neighbors[i] = cellsInSet[1];
224  else
225  neighbors[i] = cellsInSet[0];
226  }
227  // fill the rest with 0
228  for( int i = nsides; i < max_edges; i++ )
229  neighbors[i] = 0;
230  // now simply set the neighbors tag; the last few positions will not be used, but for
231  // simplicity will keep them all (MAXEDGES)
232  rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
233  }
234  return MB_SUCCESS;
235 }

References moab::Range::begin(), moab::Interface::contains_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), moab::Interface::INTERSECT, moab::Interface::list_entities(), mb, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by createTags().

◆ filterByMask()

ErrorCode moab::Intx2Mesh::filterByMask ( Range cells)
virtual

Definition at line 817 of file Intx2Mesh.cpp.

818 {
819  if( !imaskTag ) return MB_SUCCESS; // nothing to do
820  size_t sz = cells.size();
821  std::vector< int > masks( sz );
822 
823  ErrorCode rval = mb->tag_get_data( imaskTag, cells, &masks[0] );MB_CHK_ERR( rval );
824  Range cellsToRemove;
825  size_t indx = 0;
826  for( Range::iterator eit = cells.begin(); eit != cells.end(); ++eit, ++indx )
827  {
828  if( masks[indx] ) continue;
829  cellsToRemove.insert( *eit );
830  }
831  cells = subtract( cells, cellsToRemove );
832  return MB_SUCCESS;
833 }

References moab::Range::begin(), moab::Range::end(), ErrorCode, imaskTag, moab::Range::insert(), mb, MB_CHK_ERR, MB_SUCCESS, moab::Range::size(), moab::subtract(), and moab::Interface::tag_get_data().

Referenced by intersect_meshes(), and intersect_meshes_kdtree().

◆ FindMaxEdges()

ErrorCode moab::Intx2Mesh::FindMaxEdges ( EntityHandle  set1,
EntityHandle  set2 
)
virtual

Definition at line 85 of file Intx2Mesh.cpp.

86 {
87  ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" );
88  rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" );
89 
90  return MB_SUCCESS;
91 }

References ErrorCode, FindMaxEdgesInSet(), max_edges_1, max_edges_2, MB_CHK_SET_ERR, and MB_SUCCESS.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

◆ FindMaxEdgesInSet()

ErrorCode moab::Intx2Mesh::FindMaxEdgesInSet ( EntityHandle  eset,
int &  max_edges 
)
virtual

Definition at line 55 of file Intx2Mesh.cpp.

56 {
57  Range cells;
58  ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );
59 
60  max_edges = 0; // can be 0 for point clouds
61  for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
62  {
63  EntityHandle cell = *cit;
64  const EntityHandle* conn4;
65  int nnodes = 3;
66  rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
67  if( nnodes > max_edges ) max_edges = nnodes;
68  }
69  // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
70  // global) but it is better to be consistent
71 #ifdef MOAB_HAVE_MPI
72  if( parcomm )
73  {
74  int local_max_edges = max_edges;
75  // now reduce max_edges over all processors
76  int mpi_err =
77  MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
78  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
79  }
80 #endif
81 
82  return MB_SUCCESS;
83 }

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, and MB_SUCCESS.

Referenced by FindMaxEdges().

◆ findNodes()

virtual ErrorCode moab::Intx2Mesh::findNodes ( EntityHandle  tgt,
int  nsTgt,
EntityHandle  src,
int  nsSrc,
double *  iP,
int  nP 
)
pure virtual

◆ intersect_meshes()

ErrorCode moab::Intx2Mesh::intersect_meshes ( EntityHandle  mbs1,
EntityHandle  mbs2,
EntityHandle outputSet 
)

Definition at line 455 of file Intx2Mesh.cpp.

456 {
457  ErrorCode rval;
458  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
459  mbs2 = mbset2;
460  outSet = outputSet;
461 #ifdef VERBOSE
462  std::stringstream ffs, fft;
463  ffs << "source_rank0" << my_rank << ".vtk";
464  rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
465  fft << "target_rank0" << my_rank << ".vtk";
466  rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
467 
468 #endif
469  // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
470 
471  EntityHandle startSrc = 0, startTgt = 0;
472 
473  rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
474  rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
475 
476  // filter rs1 and rs2 by mask; remove everything with 0 mask
477  // get the mask tag if it exists; if not, leave it uninitialized (NULL)
478  mb->tag_get_handle( "GRID_IMASK", imaskTag );
479  if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
480 
481  rval = filterByMask( rs1 );MB_CHK_ERR( rval );
482  rval = filterByMask( rs2 );MB_CHK_ERR( rval );
483 
484  createTags(); // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
485 
486  Range rs22 = rs2; // a copy of the initial range; we will remove from it elements as we
487  // advance ; rs2 is needed for marking the polygon to the tgt parent
488 
489  // create the local kdd tree with source elements; will use it to search
490  // more efficiently for the seeds in advancing front;
491  // some of the target cells will not be covered by source cells, and they need to be eliminated
492  // early from contention
493 
494  // build a kd tree with the rs1 (source) cells
495  FileOptions kdOpts("PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;");
496  AdaptiveKDTree kd( mb );
497  kd.parse_options( kdOpts );
498  EntityHandle tree_root = 0;
499  rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
500 
501  while( !rs22.empty() )
502  {
503 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
504  if( rs22.size() < rs2.size() )
505  {
506  std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
507  << "\n";
508  std::stringstream ffo;
509  ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
510  rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
511  }
512 #endif
513  bool seedFound = false;
514  for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
515  {
516  startTgt = *it;
517  int found = 0;
518  // find vertex positions
519  const EntityHandle* conn = NULL;
520  int nnodes = 0;
521  rval = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
522  // find leaves close to those positions
523  std::vector< double > positions;
524  positions.resize( nnodes * 3 );
525  rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
526  // find leaves within a distance from each vertex of target
527  // in those leaves, collect all cells; we will try for an intx in there, instead of
528  // looping over all rs1 cells, as before
529  Range close_source_cells;
530  std::vector< EntityHandle > leaves;
531  for( int i = 0; i < nnodes; i++ )
532  {
533  leaves.clear();
534  rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
535 
536  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
537  {
538  Range tmp;
539  rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
540 
541  close_source_cells.merge( tmp.begin(), tmp.end() );
542  }
543  }
544 
545  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
546  {
547  startSrc = *it2;
548  double area = 0;
549  // if area is > 0 , we have intersections
550  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
551  //
552  int nP = 0;
553  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
554  int nsTgt, nsSrc;
555  rval =
556  computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
557  if( area > 0 )
558  {
559  found = 1;
560  seedFound = true;
561  break; // found 2 elements that intersect; these will be the seeds
562  }
563  }
564  if( found )
565  break;
566  else
567  {
568 #if defined( VERBOSE )
569  std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
570  << " not intx with any source\n";
571 #endif
572  rs22.erase( startTgt );
573  }
574  }
575  if( !seedFound ) continue; // continue while(!rs22.empty())
576 
577  std::queue< EntityHandle > srcQueue; // these are corresponding to Ta,
578  srcQueue.push( startSrc );
579  std::queue< EntityHandle > tgtQueue;
580  tgtQueue.push( startTgt );
581 
582  Range toResetSrcs; // will be used to reset src flags for every tgt element processed
583 
584  /*if (my_rank==0)
585  dbg_1 = 1;*/
586  unsigned char used = 1;
587  // mark the start tgt quad as used, so it will not come back again
588  rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
589  while( !tgtQueue.empty() )
590  {
591  // flags for the side : 0 means a src cell not found on side
592  // a paired src not found yet for the neighbors of tgt
593  Range nextSrc[MAXEDGES]; // there are new ranges of possible next src cells for
594  // seeding the side j of tgt cell
595 
596  EntityHandle currentTgt = tgtQueue.front();
597  tgtQueue.pop();
598  int nsidesTgt; // will be initialized now
599  double areaTgtCell = setup_tgt_cell( currentTgt, nsidesTgt ); // this is the area in the gnomonic plane
600  double recoveredArea = 0;
601  // get the neighbors of tgt, and if they are solved already, do not bother with that
602  // side of tgt
603  EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
604  rval = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
605 #ifdef ENABLE_DEBUG
606  if( dbg_1 )
607  {
608  std::cout << "Next: neighbors for current tgt ";
609  for( int kk = 0; kk < nsidesTgt; kk++ )
610  {
611  if( tgtNeighbors[kk] > 0 )
612  std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
613  else
614  std::cout << 0 << " ";
615  }
616  std::cout << std::endl;
617  }
618 #endif
619  // now get the status of neighbors; if already solved, make them 0, so not to bother
620  // anymore on that side of tgt
621  for( int j = 0; j < nsidesTgt; j++ )
622  {
623  EntityHandle tgtNeigh = tgtNeighbors[j];
624  unsigned char status = 1;
625  if( tgtNeigh == 0 ) continue;
626  rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval ); // status 0 is unused
627  if( 1 == status ) tgtNeighbors[j] = 0; // so will not look anymore on this side of tgt
628  }
629 
630 #ifdef ENABLE_DEBUG
631  if( dbg_1 )
632  {
633  std::cout << "reset sources: ";
634  for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
635  std::cout << mb->id_from_handle( *itr ) << " ";
636  std::cout << std::endl;
637  }
638 #endif
639  EntityHandle currentSrc = srcQueue.front();
640  // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
641  // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
642  // pairs;
643  // at every moment, the queue contains pairs of cells that intersect, and they form the
644  // "advancing front"
645  srcQueue.pop();
646  toResetSrcs.clear(); // empty the range of used srcs, will have to be set unused again,
647  // at the end of tgt element processing
648  toResetSrcs.insert( currentSrc );
649  // mb2->set_tag_data
650  std::queue< EntityHandle > localSrc;
651  localSrc.push( currentSrc );
652 #ifdef VERBOSE
653  int countingStart = counting;
654 #endif
655  // will advance-front search in the neighborhood of tgt cell, until we finish processing
656  // all
657  // possible src cells; localSrc queue will contain all possible src cells that cover
658  // the current tgt cell
659  while( !localSrc.empty() )
660  {
661  //
662  EntityHandle srcT = localSrc.front();
663  localSrc.pop();
664  double P[10 * MAXEDGES], area; //
665  int nP = 0;
666  int nb[MAXEDGES] = { 0 };
667  int nr[MAXEDGES] = { 0 };
668 
669  int nsidesSrc; ///
670  // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
671  // intersection points could include the vertices of initial elements
672  // nb [j] = 0 means no intersection on the side j for element src (markers)
673  // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
674  // tgt poly. A potential next poly in the tgt queue is the tgt poly that is
675  // adjacent to this side
676  rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
677  nsidesTgt );MB_CHK_ERR( rval );
678  if( nP > 0 )
679  {
680 #ifdef ENABLE_DEBUG
681  if( dbg_1 )
682  {
683  for( int k = 0; k < 3; k++ )
684  {
685  std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
686  }
687  }
688 #endif
689 
690  // intersection found: output P and original triangles if nP > 2
691  EntityHandle neighbors[MAXEDGES] = { 0 };
692  rval = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
693  if( rval != MB_SUCCESS )
694  {
695  std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
696  return MB_FAILURE;
697  }
698 
699  // add neighbors to the localSrc queue, if they are not marked
700  for( int nn = 0; nn < nsidesSrc; nn++ )
701  {
702  EntityHandle neighbor = neighbors[nn];
703  if( neighbor > 0 && nb[nn] > 0 ) // advance across src boundary nn
704  {
705  if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
706  {
707  localSrc.push( neighbor );
708 #ifdef ENABLE_DEBUG
709  if( dbg_1 )
710  {
711  std::cout << " local src elem " << mb->id_from_handle( neighbor )
712  << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
713  mb->list_entities( &neighbor, 1 );
714  }
715 #endif
716  toResetSrcs.insert( neighbor );
717  }
718  }
719  }
720  // n(find(nc>0))=ac; % ac is starting candidate for neighbor
721  for( int nn = 0; nn < nsidesTgt; nn++ )
722  {
723  if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
724  nextSrc[nn].insert( srcT ); // potential src cell that can intersect
725  // the tgt neighbor nn
726  }
727  if( nP > 1 )
728  { // this will also construct triangles/polygons in the new mesh, if needed
729  rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
730  }
731 
732  recoveredArea += area;
733  }
734 #ifdef ENABLE_DEBUG
735  else if( dbg_1 )
736  {
737  std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
738  << mb->id_from_handle( srcT ) << "\n";
739  }
740 #endif
741  } // end while (!localSrc.empty())
742  recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fraction
743 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
744  if( fabs( recoveredArea ) > epsilon_1 )
745  {
746 #ifdef VERBOSE
747  std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
748  << " fraction error recovery:" << recoveredArea
749  << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
750  << "\n";
751 #endif
752  }
753 #endif
754  // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
755  // mesh)
756  rs22.erase( currentTgt );
757  // also, look at its neighbors, and add to the seeds a next one
758 
759  for( int j = 0; j < nsidesTgt; j++ )
760  {
761  EntityHandle tgtNeigh = tgtNeighbors[j];
762  if( tgtNeigh == 0 || nextSrc[j].size() == 0 ) // if tgt is bigger than src, there could be no src
763  // to advance on that side
764  continue;
765  int nsidesTgt2 = 0;
766  setup_tgt_cell( tgtNeigh,
767  nsidesTgt2 ); // find possible intersection with src cell from nextSrc
768  for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
769  {
770  EntityHandle nextB = *nit;
771  // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
772  // src quad
773  double P[10 * MAXEDGES], area; //
774  int nP = 0;
775  int nb[MAXEDGES] = { 0 };
776  int nr[MAXEDGES] = { 0 };
777 
778  int nsidesSrc; ///
780  /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
781  if( area > 0 )
782  {
783  tgtQueue.push( tgtNeigh );
784  srcQueue.push( nextB );
785 #ifdef ENABLE_DEBUG
786  if( dbg_1 )
787  std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
788  << mb->id_from_handle( nextB ) << std::endl;
789 #endif
790  rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
791  break; // so we are done with this side of tgt, we have found a proper next
792  // seed
793  }
794  }
795  }
796 
797  } // end while (!tgtQueue.empty())
798  }
799 #ifdef ENABLE_DEBUG
800  if( dbg_1 )
801  {
802  for( int k = 0; k < 6; k++ )
803  mout_1[k].close();
804  }
805 #endif
806  // before cleaning up , we need to settle the position of the intersection points
807  // on the boundary edges
808  // this needs to be collective, so we should maybe wait something
809 #ifdef MOAB_HAVE_MPI
810  rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
811 #endif
812 
813  this->clean();
814  return MB_SUCCESS;
815 }

References moab::Range::begin(), moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), counting, createTags(), moab::AdaptiveKDTree::distance_search(), moab::Range::empty(), moab::Range::end(), epsilon_1, moab::Range::erase(), ErrorCode, filterByMask(), moab::Range::find(), findNodes(), moab::Range::front(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::id_from_handle(), moab::ID_FROM_HANDLE(), imaskTag, moab::Range::insert(), moab::Interface::list_entities(), MAXEDGES, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, mbs1, mbs2, moab::Range::merge(), my_rank, nr, outSet, moab::AdaptiveKDTree::parse_options(), rs1, rs2, setup_tgt_cell(), size, moab::Range::size(), srcNeighTag, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), TgtFlagTag, tgtNeighTag, and moab::Interface::write_mesh().

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().

◆ intersect_meshes_kdtree()

ErrorCode moab::Intx2Mesh::intersect_meshes_kdtree ( EntityHandle  mbset1,
EntityHandle  mbset2,
EntityHandle outputSet 
)

Definition at line 239 of file Intx2Mesh.cpp.

240 {
241  ErrorCode rval;
242  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
243  mbs2 = mbset2;
244  outSet = outputSet;
245  rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
246  rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
247  // from create tags, copy relevant ones
250  if( countTag ) mb->tag_delete( countTag );
251 
252  // filter rs1 and rs2 by mask; remove everything with 0 mask
253  // get the mask tag if it exists; if not, leave it uninitialized (NULL)
254  rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
255  if( imaskTag != NULL && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
256  rval = filterByMask( rs1 );MB_CHK_ERR( rval );
257  rval = filterByMask( rs2 );MB_CHK_ERR( rval );
258  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
259  // this is the only call that is potentially NlogN, in the whole method
260  rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
261 
262  int index = 0;
263  extraNodesVec.resize( TgtEdges.size() );
264  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, index++ )
265  {
266  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
267  extraNodesVec[index] = nv;
268  }
269 
270  int defaultInt = -1;
271  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
272  &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
273 
274  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
275  &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
276 
277  rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
278 
279  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
280  // them each time edges were for sure created before (tgtEdges)
281  // if we have a tag with this name, it could be of a different size, so delete it
282  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
284  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
285  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
286  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
287 
288  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
289  {
290  EntityHandle tgtCell = *rit;
291  int num_nodes = 0;
292  rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
293  // account for padded polygons
294  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
295  num_nodes--;
296 
297  for( int i = 0; i < num_nodes; i++ )
298  {
299  EntityHandle v[2] = { tgtConn[i],
300  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
301  std::vector< EntityHandle > adj_entities;
302  rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
303  if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
304  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
305  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
306  // the tag is dense
307  }
308  // now set the value of the tag
309  rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
310  }
311 
312  // find out max edge on source mesh;
313  double max_length = 0;
314  {
315  std::vector< double > coords( 3 * max_edges_1, 0.0 );
316  for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
317  {
318  const EntityHandle* conn = NULL;
319  int nnodes;
320  rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
321  while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
322  nnodes--;
323  rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" );
324  for( int j = 0; j < nnodes; j++ )
325  {
326  int next = ( j + 1 ) % nnodes;
327  double edge_length =
328  ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
329  ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
330  ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
331  if( edge_length > max_length ) max_length = edge_length;
332  }
333  }
334  max_length = std::sqrt( max_length );
335  }
336 
337  // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
338  double tolerance = 1.e-15;
339  if( max_length < 1. )
340  {
341  // basically, the sag for an arc of length max_length on a circle of radius 1
342  tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
344  tolerance = 3 * tolerance; // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
345  // be more generous, use 1.5 ~= sqrt(2.)
346 
347  if( !my_rank )
348  {
349  std::cout << " max edge length: " << max_length << " tolerance for kd tree: " << tolerance << "\n";
350  std::cout << " box overlap tolerance: " << box_error << "\n";
351  }
352  }
353 #ifdef MOAB_HAVE_MPI
354  // reduce box tolerance on every task, if needed
355  double min_box_eps;
356  MPI_Allreduce( &box_error, &min_box_eps, 1, MPI_DOUBLE, MPI_MIN, parcomm->comm() );
357  box_error = min_box_eps;
358 #endif
359 
360  // create the kd tree on source cells, and intersect all targets in an expensive loop
361  // build a kd tree with the rs1 (source) cells
362  FileOptions kdOpts("PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;");
363  AdaptiveKDTree kd( mb );
364  kd.parse_options( kdOpts );
365  EntityHandle tree_root = 0;
366  rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
367 
368  for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
369  {
370  EntityHandle tcell = *it;
371  // find vertex positions
372  const EntityHandle* conn = NULL;
373  int nnodes = 0;
374  rval = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
375  // find leaves close to those positions
376  double areaTgtCell = setup_tgt_cell( tcell, nnodes ); // this is the area in the gnomonic plane
377  double recoveredArea = 0;
378  std::vector< double > positions;
379  positions.resize( nnodes * 3 );
380  rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
381 
382  // distance to search will be based on average edge length
383  double av_len = 0;
384  for( int k = 0; k < nnodes; k++ )
385  {
386  int ik = ( k + 1 ) % nnodes;
387  double len1 = 0;
388  for( int j = 0; j < 3; j++ )
389  {
390  double len2 = positions[3 * k + j] - positions[3 * ik + j];
391  len1 += len2 * len2;
392  }
393  av_len += sqrt( len1 );
394  }
395  if( nnodes > 0 ) av_len /= nnodes;
396  // find leaves within a distance from each vertex of target
397  // in those leaves, collect all cells; we will try for an intx in there
398  Range close_source_cells;
399  std::vector< EntityHandle > leaves;
400  for( int i = 0; i < nnodes; i++ )
401  {
402  leaves.clear();
403  rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
404 
405  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
406  {
407  Range tmp;
408  rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
409 
410  close_source_cells.merge( tmp.begin(), tmp.end() );
411  }
412  }
413 #ifdef VERBOSE
414  if( close_source_cells.empty() )
415  {
416  std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
417  << mb->id_from_handle( tcell ) << "\n";
418  }
419 #endif
420  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
421  {
422  EntityHandle startSrc = *it2;
423  double area = 0;
424  // if area is > 0 , we have intersections
425  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
426  //
427  int nP = 0;
428  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
429  int nsTgt, nsSrc;
430  rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
431  if( area > 0 )
432  {
433  if( nP > 1 )
434  { // this will also construct triangles/polygons in the new mesh, if needed
435  rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
436  }
437  recoveredArea += area;
438  }
439  }
440  recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fract
441  }
442  // before cleaning up , we need to settle the position of the intersection points
443  // on the boundary edges
444  // this needs to be collective, so we should maybe wait something
445 #ifdef MOAB_HAVE_MPI
446  rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
447 #endif
448 
449  this->clean();
450  return MB_SUCCESS;
451 }

References moab::Range::begin(), box_error, moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), countTag, moab::AdaptiveKDTree::distance_search(), edge_length(), moab::Range::empty(), moab::Range::end(), epsilon_1, ErrorCode, extraNodesVec, filterByMask(), findNodes(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::id_from_handle(), imaskTag, moab::Interface::INTERSECT, max_edges_1, max_edges_2, MAXEDGES, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, MB_TYPE_INTEGER, mbs1, mbs2, moab::Range::merge(), my_rank, neighTgtEdgeTag, nr, outSet, moab::AdaptiveKDTree::parse_options(), rs1, rs2, setup_tgt_cell(), moab::Range::size(), srcParentTag, moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), tgtConn, TgtEdges, tgtParentTag, moab::tolerance, and moab::Interface::UNION.

Referenced by moab::TempestRemapper::ComputeOverlapMesh().

◆ set_box_error()

void moab::Intx2Mesh::set_box_error ( double  berror)
inline

Definition at line 138 of file Intx2Mesh.hpp.

139  {
140  box_error = berror;
141  }

References box_error.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

◆ set_error_tolerance()

void moab::Intx2Mesh::set_error_tolerance ( double  eps)
inline

Definition at line 121 of file Intx2Mesh.hpp.

122  {
123  epsilon_1 = eps;
124  epsilon_area = eps * sqrt( eps );
125  }

References epsilon_1, and epsilon_area.

Referenced by moab::TempestRemapper::ConstructCoveringSet(), and main().

◆ setup_tgt_cell()

virtual double moab::Intx2Mesh::setup_tgt_cell ( EntityHandle  tgt,
int &  nsTgt 
)
pure virtual

Member Data Documentation

◆ allBoxes

std::vector< double > moab::Intx2Mesh::allBoxes
protected

Definition at line 253 of file Intx2Mesh.hpp.

◆ box_error

◆ counting

int moab::Intx2Mesh::counting
protected

◆ countTag

◆ epsilon_1

◆ epsilon_area

◆ extraNodesVec

std::vector< std::vector< EntityHandle >* > moab::Intx2Mesh::extraNodesVec
protected

◆ gid

Tag moab::Intx2Mesh::gid
protected

◆ imaskTag

Tag moab::Intx2Mesh::imaskTag
protected

for coverage mesh, will store the original sender

Definition at line 226 of file Intx2Mesh.hpp.

Referenced by filterByMask(), intersect_meshes(), and intersect_meshes_kdtree().

◆ localEnts

Range moab::Intx2Mesh::localEnts
protected

Definition at line 257 of file Intx2Mesh.hpp.

◆ localRoot

EntityHandle moab::Intx2Mesh::localRoot
protected

Definition at line 256 of file Intx2Mesh.hpp.

◆ max_edges_1

int moab::Intx2Mesh::max_edges_1
protected

◆ max_edges_2

int moab::Intx2Mesh::max_edges_2
protected

◆ mb

◆ mbs1

EntityHandle moab::Intx2Mesh::mbs1
protected

Definition at line 199 of file Intx2Mesh.hpp.

Referenced by createTags(), intersect_meshes(), and intersect_meshes_kdtree().

◆ mbs2

EntityHandle moab::Intx2Mesh::mbs2
protected

Definition at line 200 of file Intx2Mesh.hpp.

Referenced by createTags(), intersect_meshes(), and intersect_meshes_kdtree().

◆ my_rank

unsigned int moab::Intx2Mesh::my_rank
protected

◆ neighTgtEdgeTag

◆ orgSendProcTag

Tag moab::Intx2Mesh::orgSendProcTag
protected

Definition at line 225 of file Intx2Mesh.hpp.

Referenced by moab::Intx2MeshOnSphere::findNodes().

◆ outSet

◆ rs1

◆ rs2

◆ srcConn

◆ srcCoords

◆ srcCoords2D

◆ srcNeighTag

Tag moab::Intx2Mesh::srcNeighTag
protected

Definition at line 218 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

◆ srcParentTag

◆ tgtConn

◆ tgtCoords

◆ tgtCoords2D

◆ TgtEdges

◆ TgtFlagTag

Tag moab::Intx2Mesh::TgtFlagTag
protected

Definition at line 208 of file Intx2Mesh.hpp.

Referenced by clean(), createTags(), and intersect_meshes().

◆ tgtNeighTag

Tag moab::Intx2Mesh::tgtNeighTag
protected

Definition at line 220 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

◆ tgtParentTag


The documentation for this class was generated from the following files: