Mesh Oriented datABase  (version 5.6.0)
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

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
 

Detailed Description

Definition at line 50 of file Intx2Mesh.hpp.

Constructor & Destructor Documentation

◆ Intx2Mesh()

moab::Intx2Mesh::Intx2Mesh ( Interface *  mbimpl)

Definition at line 43 of file Intx2Mesh.cpp.

44  : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
45  countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), imaskTag( 0 ),
46  tgtConn( nullptr ), srcConn( nullptr ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ),
47  my_rank( 0 )
48 #ifdef MOAB_HAVE_MPI
49  ,
50  parcomm( nullptr ), remote_cells( nullptr ), remote_cells_with_tracers( nullptr )
51 #endif
52  ,
53  max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
54 {
55  gid = mbimpl->globalId_tag();
56 }

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

◆ ~Intx2Mesh()

moab::Intx2Mesh::~Intx2Mesh ( )
virtual

Definition at line 58 of file Intx2Mesh.cpp.

59 {
60  // TODO Auto-generated destructor stub
61 #ifdef MOAB_HAVE_MPI
62  if( remote_cells )
63  {
64  delete remote_cells;
65  remote_cells = nullptr;
66  }
67 #endif
68 }

Member Function Documentation

◆ clean()

void moab::Intx2Mesh::clean ( )

Definition at line 1025 of file Intx2Mesh.cpp.

1026 {
1027  //
1028  int indx = 0;
1029  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
1030  {
1031  delete extraNodesVec[indx];
1032  }
1033  // extraNodesMap.clear();
1034  extraNodesVec.clear();
1035  // also, delete some bit tags, used to mark processed tgts and srcs
1036  mb->tag_delete( TgtFlagTag );
1037  counting = 0; // reset counting to original value
1038 }

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 1043 of file Intx2Mesh.cpp.

1044 {
1045  int i = 0;
1046  while( i < nP )
1047  {
1048  int nextIndex = ( i + 1 ) % nP;
1049  if( nodes[i] == nodes[nextIndex] )
1050  {
1051 #ifdef ENABLE_DEBUG
1052  // we need to reduce nP, and collapse nodes
1053  if( dbg_1 )
1054  {
1055  std::cout << " nodes duplicated in list: ";
1056  for( int j = 0; j < nP; j++ )
1057  std::cout << nodes[j] << " ";
1058  std::cout << "\n";
1059  std::cout << " node " << nodes[i] << " at index " << i << " is duplicated" << "\n";
1060  }
1061 #endif
1062  // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
1063  // then next thing does nothing
1064  // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
1065  for( int k = i; k < nP - 1; k++ )
1066  nodes[k] = nodes[k + 1];
1067  nP--; // decrease the number of nodes; also, decrease i, just if we may need to check
1068  // again
1069  i--;
1070  }
1071  i++;
1072  }
1073  return;
1074 }

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 128 of file Intx2Mesh.cpp.

129 {
132  if( countTag ) mb->tag_delete( countTag );
133 
134  unsigned char def_data_bit = 0; // unused by default
135  // maybe the tgt tag is better to be deleted every time, and recreated;
136  // or is it easy to set all values to something again? like 0?
137  MB_CHK_SET_ERR( mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit ),
138  "can't get tgt flag tag" );
139  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
140  // this is the only call that is potentially NlogN, in the whole method
141  MB_CHK_SET_ERR( mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION ), "can't get adjacent tgt edges" );
142 
143  // now, create a map from each edge to a list of potential new nodes on a tgt edge
144  // this memory has to be cleaned up
145  // change it to a vector, and use the index in range of tgt edges
146  int indx = 0;
147  extraNodesVec.resize( TgtEdges.size() );
148  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
149  {
150  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
151  extraNodesVec[indx] = nv;
152  }
153 
154  int defaultInt = -1;
155 
157  &defaultInt ),
158  "can't create TargetParent tag" );
159 
161  &defaultInt ),
162  "can't create SourceParent tag" );
163 
165  &defaultInt ),
166  "can't create Counting tag" );
167 
168  // for each cell in set 1, determine its neigh in set 1 (could be nullptr too)
169  // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
171  "can't determine neighbors for set 1" );
173  "can't determine neighbors for set 2" );
174 
175  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
176  // them each time edges were for sure created before (tgtEdges)
177  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
178  // if we have a tag with this name, it could be of a different size, so delete it if it exists
179  if( mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag ) == MB_SUCCESS ) mb->tag_delete( neighTgtEdgeTag );
181  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] ),
182  "can't create target edge neighbors tag" );
183  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); ++rit )
184  {
185  EntityHandle tgtCell = *rit;
186  int num_nodes = 0;
187  MB_CHK_SET_ERR( mb->get_connectivity( tgtCell, tgtConn, num_nodes ), "can't get target connectivity" );
188  // account for padded polygons
189  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
190  num_nodes--;
191 
192  int i = 0;
193  for( i = 0; i < num_nodes; i++ )
194  {
195  EntityHandle v[2] = { tgtConn[i],
196  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
197  std::vector< EntityHandle > adj_entities;
198  MB_CHK_SET_ERR( mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT ),
199  "can't get adjacencies" );
200  if( !adj_entities.size() ) MB_CHK_SET_ERR( MB_FAILURE, "no adjacencies found" ); // get out , big error
201  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
202  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
203  // the tag is dense
204  }
205  // zero out the rest
206  for( i = num_nodes; i < max_edges_2; i++ )
207  zeroh[i] = 0;
208  // now set the value of the tag
209  MB_CHK_SET_ERR( mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) ),
210  "can't set edge target edge neighbors tag" );
211  }
212  return MB_SUCCESS;
213 }

References moab::Range::begin(), countTag, DetermineOrderedNeighbors(), moab::Range::end(), 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 215 of file Intx2Mesh.cpp.

216 {
217  Range cells;
218  MB_CHK_SET_ERR( mb->get_entities_by_dimension( inputSet, 2, cells ), "can't get cells in set" );
219 
220  std::vector< EntityHandle > neighbors( max_edges );
221  std::vector< EntityHandle > zeroh( max_edges, 0 );
222  // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
223  // one on src mesh
225  &zeroh[0] ),
226  "can't create neighbors tag" );
227 
228  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
229  {
230  EntityHandle cell = *cit;
231  int nnodes = 3;
232  // will get the nnodes ordered neighbors;
233  // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
234  const EntityHandle* conn4;
235  MB_CHK_SET_ERR( mb->get_connectivity( cell, conn4, nnodes ), "can't get connectivity of a cell" );
236  int nsides = nnodes;
237  // account for possible padded polygons
238  while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
239  nsides--;
240 
241  for( int i = 0; i < nsides; i++ )
242  {
243  EntityHandle v[2];
244  v[0] = conn4[i];
245  v[1] = conn4[( i + 1 ) % nsides];
246  // get all cells adjacent to these 2 vertices on the edge
247  std::vector< EntityHandle > adjcells;
248  std::vector< EntityHandle > cellsInSet;
249  MB_CHK_SET_ERR( mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT ),
250  "can't get adjacency to 2 verts" );
251  // now look for the cells contained in the input set;
252  // the input set should be a correct mesh, not overlapping cells, and manifold
253  size_t siz = adjcells.size();
254  for( size_t j = 0; j < siz; j++ )
255  if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
256  siz = cellsInSet.size();
257 
258  if( siz > 2 )
259  {
260  std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
261  << std::endl;
262  MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" ); // non-manifold
263  }
264  if( siz == 1 )
265  {
266  // it must be the border of the input mesh;
267  neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
268  // borders do not appear for a sphere in serial, but they do appear for
269  // parallel processing anyway
270  continue;
271  }
272  // here siz ==2, it is either the first or second
273  if( cell == cellsInSet[0] )
274  neighbors[i] = cellsInSet[1];
275  else
276  neighbors[i] = cellsInSet[0];
277  }
278  // fill the rest with 0
279  for( int i = nsides; i < max_edges; i++ )
280  neighbors[i] = 0;
281  // now simply set the neighbors tag; the last few positions will not be used, but for
282  // simplicity will keep them all (MAXEDGES)
283  MB_CHK_SET_ERR( mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] ), "can't set neighbors tag" );
284  }
285  return MB_SUCCESS;
286 }

References moab::Range::begin(), moab::Interface::contains_entities(), moab::Range::end(), 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 1006 of file Intx2Mesh.cpp.

1007 {
1008  if( !imaskTag ) return MB_SUCCESS; // nothing to do
1009  size_t sz = cells.size();
1010  std::vector< int > masks( sz );
1011 
1012  MB_CHK_SET_ERR( mb->tag_get_data( imaskTag, cells, &masks[0] ), "can't get mask tag" );
1013  Range cellsToRemove;
1014  size_t indx = 0;
1015  for( Range::iterator eit = cells.begin(); eit != cells.end(); ++eit, ++indx )
1016  {
1017  if( masks[indx] ) continue;
1018  cellsToRemove.insert( *eit );
1019  }
1020  cells = subtract( cells, cellsToRemove );
1021  return MB_SUCCESS;
1022 }

References moab::Range::begin(), moab::Range::end(), imaskTag, moab::Range::insert(), mb, MB_CHK_SET_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 100 of file Intx2Mesh.cpp.

101 {
102  MB_CHK_SET_ERR( FindMaxEdgesInSet( set1, max_edges_1 ), "can't determine max_edges in set 1" );
103  MB_CHK_SET_ERR( FindMaxEdgesInSet( set2, max_edges_2 ), "can't determine max_edges in set 2" );
104 
105  // Numerous routines below use fixed-size stack buffers dimensioned from
106  // MAXEDGES (coordinates, edge-neighbor tags, connectivity scratch, ...),
107  // but max_edges_1/2 are whatever the meshes actually contain. A cell with
108  // more vertices than MAXEDGES silently overruns those buffers, which shows
109  // up much later as corrupted geometry or a crash far from the cause.
110  // Fail here instead, where the diagnosis is obvious.
111  //
112  // Dual meshes of high-order spectral-element grids are the usual source of
113  // high-valence cells; if this triggers, either raise MAXEDGES in
114  // IntxUtils.hpp or regenerate the mesh with lower-valence cells.
115  const int max_edges = std::max( max_edges_1, max_edges_2 );
116  if( max_edges > MAXEDGES )
117  {
118  MB_SET_ERR( MB_FAILURE, "mesh contains a cell with "
119  << max_edges << " vertices, which exceeds MAXEDGES (" << MAXEDGES
120  << "); source mesh max = " << max_edges_1
121  << ", target mesh max = " << max_edges_2
122  << ". Increase MAXEDGES in moab/IntxMesh/IntxUtils.hpp and rebuild." );
123  }
124 
125  return MB_SUCCESS;
126 }

References FindMaxEdgesInSet(), max_edges_1, max_edges_2, MAXEDGES, MB_CHK_SET_ERR, MB_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 70 of file Intx2Mesh.cpp.

71 {
72  Range cells;
73  MB_CHK_SET_ERR( mb->get_entities_by_dimension( eset, 2, cells ), "can't get entities by dimension" );
74 
75  max_edges = 0; // can be 0 for point clouds
76  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
77  {
78  EntityHandle cell = *cit;
79  const EntityHandle* conn4;
80  int nnodes = 3;
81  MB_CHK_SET_ERR( mb->get_connectivity( cell, conn4, nnodes ), "can't get connectivity of a cell" );
82  if( nnodes > max_edges ) max_edges = nnodes;
83  }
84  // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
85  // global) but it is better to be consistent
86 #ifdef MOAB_HAVE_MPI
87  if( parcomm )
88  {
89  int local_max_edges = max_edges;
90  // now reduce max_edges over all processors
91  int mpi_err =
92  MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
93  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
94  }
95 #endif
96 
97  return MB_SUCCESS;
98 }

References moab::Range::begin(), moab::Range::end(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), mb, 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 580 of file Intx2Mesh.cpp.

581 {
582  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
583  mbs2 = mbset2;
584  outSet = outputSet;
585 #ifdef VERBOSE
586  std::stringstream ffs, fft;
587  ffs << "source_rank0" << my_rank << ".vtk";
588  MB_CHK_SET_ERR( mb->write_mesh( ffs.str().c_str(), &mbset1, 1 ), "can't write source mesh" );
589  fft << "target_rank0" << my_rank << ".vtk";
590  MB_CHK_SET_ERR( mb->write_mesh( fft.str().c_str(), &mbset2, 1 ), "can't write target mesh" );
591 
592 #endif
593  // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
594 
595  EntityHandle startSrc = 0, startTgt = 0;
596 
597  MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs1, 2, rs1 ), "can't get source entities by dimension" );
598  MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs2, 2, rs2 ), "can't get target entities by dimension" );
599 
600  // filter rs1 and rs2 by mask; remove everything with 0 mask
601  // get the mask tag if it exists; if not, leave it uninitialized (nullptr)
602  ErrorCode rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
603  if( imaskTag != nullptr && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
604 
605  MB_CHK_SET_ERR( filterByMask( rs1 ), "can't filter source by mask" );
606  MB_CHK_SET_ERR( filterByMask( rs2 ), "can't filter target by mask" );
607 
608  createTags(); // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
609 
610  Range rs22 = rs2; // a copy of the initial range; we will remove from it elements as we
611  // advance ; rs2 is needed for marking the polygon to the tgt parent
612 
613  // create the local kdd tree with source elements; will use it to search
614  // more efficiently for the seeds in advancing front;
615  // some of the target cells will not be covered by source cells, and they need to be eliminated
616  // early from contention
617  FileOptions kdOpts( "PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;" );
618  AdaptiveKDTree kd( mb );
619  kd.parse_options( kdOpts );
620  EntityHandle tree_root = 0;
621 
622  // build a kd tree with the rs1 (source) cells
623  MB_CHK_SET_ERR( kd.build_tree( rs1, &tree_root ), "can't build kd tree on source cells" );
624 #if defined( ENABLE_DEBUG )
625  for( auto it = rs22.begin(); it != rs22.end(); ++it )
626  {
627  EntityHandle cell = *it;
628  const EntityHandle* conn = nullptr;
629  int nnodes = 0;
630  rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_ERR( rval );
631  std::cout << " cell: \t" << " ht:" << mb->id_from_handle( cell ) << " nodes: " << nnodes
632  << " gt:" << global_id_ent( mb, cell, gid ) << " stat: " << char_stat_ent( mb, cell, TgtFlagTag )
633  << "\n";
634  }
635 #endif
636 
637  while( !rs22.empty() )
638  {
639 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
640  if( rs22.size() < rs2.size() )
641  {
642  std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
643  << " rs22.size():" << rs22.size() << "\n";
644  std::stringstream ffo;
645  ffo << "file0" << counting << "rank0" << my_rank << ".h5m";
646  MB_CHK_SET_ERR( mb->write_mesh( ffo.str().c_str(), &outSet, 1 ), "can't write output mesh" );
647  for( auto it = rs22.begin(); it != rs22.end(); ++it )
648  {
649  EntityHandle cell = *it;
650  const EntityHandle* conn = nullptr;
651  int nnodes = 0;
652  rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_ERR( rval );
653  std::cout << " cell: \t" << " ht:" << mb->id_from_handle( cell ) << " nodes: " << nnodes
654  << " g:" << global_id_ent( mb, cell, gid )
655  << " stat: " << char_stat_ent( mb, cell, TgtFlagTag ) << "\n";
656  }
657  }
658 #endif
659  bool seedFound = false;
660  Range verified_seeds;
661  for( Range::reverse_iterator it = rs22.rbegin(); it != rs22.rend(); ++it )
662  {
663  startTgt = *it;
664  unsigned char status = 0;
665  rval = mb->tag_get_data( TgtFlagTag, &startTgt, 1, &status );MB_CHK_ERR( rval );
666  if( 1 == status )
667  {
668  verified_seeds.insert( startTgt );
669  continue;
670  }
671  int found = 0;
672  // find vertex positions
673  const EntityHandle* conn = nullptr;
674  int nnodes = 0;
675  MB_CHK_SET_ERR( mb->get_connectivity( startTgt, conn, nnodes ), "can't get target connectivity" );
676  // find leaves close to those positions
677  std::vector< double > positions;
678  positions.resize( nnodes * 3 );
679  MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &positions[0] ), "can't get target coordinates" );
680  // find leaves within a distance from each vertex of target
681  // in those leaves, collect all cells; we will try for an intx in there, instead of
682  // looping over all rs1 cells, as before
683  Range close_source_cells;
684  std::vector< EntityHandle > leaves;
685  for( int i = 0; i < nnodes; i++ )
686  {
687  leaves.clear();
688  MB_CHK_SET_ERR( kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 ),
689  "can't search for leaves" );
690 
691  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
692  {
693  Range tmp;
694  MB_CHK_SET_ERR( mb->get_entities_by_dimension( *j, 2, tmp ), "can't get entities by dimension" );
695 
696  close_source_cells.merge( tmp.begin(), tmp.end() );
697  }
698  }
699 
700  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
701  {
702  startSrc = *it2;
703  double area = 0;
704  // if area is > 0 , we have intersections
705  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
706  //
707  int nP = 0;
708  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
709  int nsTgt, nsSrc;
710  MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc,
711  nsTgt, true ),
712  "can't compute intersection between target and source" );
713  if( area > 0 )
714  {
715  found = 1;
716  seedFound = true;
717  break; // found 2 elements that intersect; these will be the seeds
718  }
719  }
720  if( found )
721  break;
722  else
723  {
724 #ifdef ENABLE_DEBUG
725  std::cout << " on rank " << my_rank << " target cell " << " ht:" << mb->id_from_handle( startTgt )
726  << " g:" << global_id_ent( mb, startTgt, gid ) << " not intx with any source\n";
727 #endif
728  verified_seeds.insert( startTgt );
729  }
730  }
731  rs22 = subtract( rs22, verified_seeds );
732  if( !seedFound ) continue; // continue while(!rs22.empty())
733 
734  std::queue< EntityHandle > srcQueue; // these are corresponding to Ta,
735  srcQueue.push( startSrc );
736  std::queue< EntityHandle > tgtQueue;
737  tgtQueue.push( startTgt );
738 
739  unsigned char used = 1;
740  // mark the start tgt quad as used, so it will not come back again
741  MB_CHK_SET_ERR( mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used ), "can't set target flag" );
742  while( !tgtQueue.empty() )
743  {
744  // flags for the side : 0 means a src cell not found on side
745  // a paired src not found yet for the neighbors of tgt
746  Range nextSrc[MAXEDGES]; // there are new ranges of possible next src cells for
747  // seeding the side j of tgt cell
748 
749  EntityHandle currentTgt = tgtQueue.front();
750  tgtQueue.pop();
751  int nsidesTgt; // will be initialized now
752  double areaTgtCell = setup_tgt_cell( currentTgt, nsidesTgt ); // this is the area in the gnomonic plane
753  double recoveredArea = 0;
754  // get the neighbors of tgt, and if they are solved already, do not bother with that
755  // side of tgt
756  EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
757  MB_CHK_SET_ERR( mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors ),
758  "can't get target neighbors" );
759 #ifdef ENABLE_DEBUG
760  if( dbg_1 )
761  {
762  std::cout << "Next: neighbors for current tgt nsidesTgt: " << nsidesTgt << " ";
763  for( int kk = 0; kk < nsidesTgt; kk++ )
764  {
765  if( tgtNeighbors[kk] > 0 )
766  std::cout << " ht:" << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
767  else
768  std::cout << 0 << " ";
769  }
770  std::cout << std::endl;
771  }
772 #endif
773  // now get the status of neighbors; if already solved, make them 0, so not to bother
774  // anymore on that side of tgt
775  for( int j = 0; j < nsidesTgt; j++ )
776  {
777  EntityHandle tgtNeigh = tgtNeighbors[j];
778  unsigned char status = 1;
779  if( tgtNeigh == 0 ) continue;
780  MB_CHK_SET_ERR( mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status ),
781  "can't get target flag" ); // status 0 is unused
782  if( 1 == status ) tgtNeighbors[j] = 0; // so will not look anymore on this side of tgt
783  }
784 
785  EntityHandle currentSrc = srcQueue.front();
786  // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
787  // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
788  // pairs;
789  // at every moment, the queue contains pairs of cells that intersect, and they form the
790  // "advancing front"
791  srcQueue.pop();
792 
793  Range localSrc;
794  Range localSrcAlreadyTested;
795  localSrc.insert( currentSrc );
796 #ifdef VERBOSE
797  int countingStart = counting;
798 #endif
799  // will advance-front search in the neighborhood of tgt cell, until we finish processing
800  // all
801  // possible src cells; localSrc set will contain all possible src cells that cover
802  // the current tgt cell
803  while( !localSrc.empty() )
804  {
805  //
806  EntityHandle srcT = localSrc.pop_front(); // also remove from local range
807  double P[10 * MAXEDGES], area; //
808  int nP = 0;
809  int nb[MAXEDGES] = { 0 };
810  int nr[MAXEDGES] = { 0 };
811 
812  int nsidesSrc;
813  // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
814  // intersection points could include the vertices of initial elements
815  // nb [j] = 0 means no intersection on the side j for element src (markers)
816  // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
817  // tgt poly. A potential next poly in the tgt queue is the tgt poly that is
818  // adjacent to this side
819  MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr,
820  nsidesSrc, nsidesTgt ),
821  "can't compute intersection between target and source" );
822  localSrcAlreadyTested.insert( srcT );
823 
824  if( nP > 0 )
825  {
826 #ifdef ENABLE_DEBUG
827  std::cout << " srcT: " << " hs:" << mb->id_from_handle( srcT )
828  << " g:" << global_id_ent( mb, srcT, gid ) << " nsidesSrc:" << nsidesSrc << "\n";
829  std::cout << " currentTgt: " << " ht:" << mb->id_from_handle( currentTgt )
830  << " g:" << global_id_ent( mb, currentTgt, gid )
831  << " stat:" << char_stat_ent( mb, currentTgt, TgtFlagTag ) << " nsidesTgt:" << nsidesTgt
832  << "\n";
833  unsigned char status = 1;
834  rval = mb->tag_set_data( TgtFlagTag, &currentTgt, 1, &status );MB_CHK_ERR( rval );
835  if( dbg_1 )
836  {
837  for( int k = 0; k < nsidesSrc; k++ )
838  std::cout << " nb[" << k << "]=" << nb[k];
839  std::cout << "\n";
840  for( int k = 0; k < nsidesTgt; k++ )
841  std::cout << " nr[" << k << "]=" << nr[k];
842  std::cout << "\n";
843  }
844 #endif
845 
846  // intersection found: output P and original triangles if nP > 2
847  EntityHandle neighbors[MAXEDGES] = { 0 };
848 
849  MB_CHK_SET_ERR( mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors ),
850  "failed to get the neighbors for source element " << mb->id_from_handle( srcT ) );
851 
852  Range newPotentialSrc;
853  // add neighbors to the localSrc queue, if they are not marked
854  for( int nn = 0; nn < nsidesSrc; nn++ )
855  {
856  EntityHandle neighbor = neighbors[nn];
857  if( nb[nn] > 0 ) // advance across src boundary nn
858  {
859  if( neighbor > 0 )
860  {
861  if( localSrcAlreadyTested.index( neighbor ) < 0 ) // -1
862  {
863  localSrc.insert( neighbor );
864 #ifdef ENABLE_DEBUG
865  std::cout << " local src elem " << " hs:" << mb->id_from_handle( neighbor )
866  << " for tgt:" << "ht:" << mb->id_from_handle( currentTgt ) << "\n";
867 #endif
868  }
869  }
870  else // if it is on the boundary it is a special case, maybe we need to advance more on that side,
871  // because the boundary is non-convex
872  {
873  // find the ends of edge nn, and add adjacent sources to those vertices (if they are in rs1)
874  int NumNodesSrc = 0;
875  const EntityHandle* connS;
876  rval = mb->get_connectivity( srcT, connS, NumNodesSrc );MB_CHK_SET_ERR( rval, "can't get connectivity" );
877  EntityHandle v1 = connS[nn];
878  EntityHandle v2 = connS[( nn + 1 ) % NumNodesSrc];
879  Range adjacentSourceCells1, adjacentSourceCells2;
880  rval = mb->get_adjacencies( &v1, 1, 2, false, adjacentSourceCells1 );MB_CHK_SET_ERR( rval, "can't get adjacent cells" );
881  adjacentSourceCells1 = intersect( adjacentSourceCells1, rs1 );
882  rval = mb->get_adjacencies( &v2, 1, 2, false, adjacentSourceCells2 );MB_CHK_SET_ERR( rval, "can't get adjacent cells" );
883  adjacentSourceCells2 = intersect( adjacentSourceCells2, rs1 );
884  adjacentSourceCells1.merge( adjacentSourceCells2 );
885  Range potentialSrc = subtract( adjacentSourceCells1, localSrcAlreadyTested );
886  newPotentialSrc.merge( potentialSrc );
887  }
888  }
889  }
890  // these might come from non-convex boundary
891  if( !newPotentialSrc.empty() )
892  {
893  localSrc.merge( newPotentialSrc );
894  }
895  // n(find(nc>0))=ac; % ac is starting candidate for neighbor
896  for( int nn = 0; nn < nsidesTgt; nn++ )
897  {
898  if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
899  nextSrc[nn].insert( srcT ); // potential src cell that can intersect
900  // the tgt neighbor nn
901  }
902  if( nP > 1 )
903  { // this will also construct triangles/polygons in the new mesh, if needed
904  MB_CHK_SET_ERR( findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP ),
905  "can't find nodes" );
906 #ifdef ENABLE_DEBUG
907  std::cout << " intersect: " << " ht:" << mb->id_from_handle( currentTgt ) << " "
908  << " hs:" << mb->id_from_handle( srcT )
909  << " g:" << global_id_ent( mb, currentTgt, gid ) << " "
910  << " g:" << global_id_ent( mb, srcT, gid ) << " counting: " << counting << "\n";
911 #endif
912  }
913 
914  recoveredArea += area;
915  }
916 #ifdef ENABLE_DEBUG
917  else if( dbg_1 )
918  {
919  std::cout << " tgt, src, do not intersect: " << "ht:" << mb->id_from_handle( currentTgt ) << " "
920  << "hs:" << mb->id_from_handle( srcT ) << "\n";
921  }
922 #endif
923  } // end while (!localSrc.empty())
924  recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fraction
925 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
926  if( fabs( recoveredArea ) > epsilon_1 )
927  {
928 #ifdef VERBOSE
929  std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
930  << " fraction error recovery:" << recoveredArea
931  << " tgtID: " << "ht:" << mb->id_from_handle( currentTgt )
932  << " countingStart:" << countingStart << "\n";
933 #endif
934  }
935 #endif
936  // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
937  // mesh)
938  rs22.erase( currentTgt );
939 #ifdef ENABLE_DEBUG
940  std::cout << " remove cell: " << "ht:" << mb->id_from_handle( currentTgt )
941  << " g:" << global_id_ent( mb, currentTgt, gid ) << " rs22.size():" << rs22.size() << "\n";
942 #endif
943  // also, look at its neighbors, and add to the seeds a next one
944 
945  //int savedGnoPlane = plane;
946  for( int j = 0; j < nsidesTgt; j++ )
947  {
948  EntityHandle tgtNeigh = tgtNeighbors[j];
949  if( tgtNeigh == 0 || nextSrc[j].size() == 0 ) // if tgt is bigger than src, there could be no src
950  // to advance on that side
951  continue;
952  int nsidesTgt2 = 0;
953  // this also changes gnomonic plane //plane//
954  setup_tgt_cell( tgtNeigh, nsidesTgt2 ); // find possible intersection with src cell from nextSrc
955  for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
956  {
957  EntityHandle nextB = *nit;
958  // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
959  // src quad
960  double P[10 * MAXEDGES], area; //
961  int nP = 0;
962  int nb[MAXEDGES] = { 0 };
963  int nr[MAXEDGES] = { 0 };
964 
965  int nsidesSrc; ///
967  /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 ),
968  "can't compute intersection between target and source" );
969  if( area > 0 )
970  {
971  unsigned char is_used = 0;
972  rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &is_used );MB_CHK_ERR( rval );
973  if( 0 == is_used )
974  {
975  tgtQueue.push( tgtNeigh );
976  srcQueue.push( nextB );
977 #ifdef ENABLE_DEBUG
978  if( dbg_1 )
979  std::cout << "new polys pushed: src, tgt:" << " ht:" << mb->id_from_handle( tgtNeigh )
980  << " hs:" << mb->id_from_handle( nextB ) << " counting: " << counting
981  << std::endl;
982 #endif
983  MB_CHK_SET_ERR( mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used ),
984  "can't set target flag" );
985  }
986  break; // so we are done with this side of tgt, we have found a proper next
987  // seed
988  }
989  }
990  }
991 
992  } // end while (!tgtQueue.empty())
993  }
994 
995  // before cleaning up , we need to settle the position of the intersection points
996  // on the boundary edges
997  // this needs to be collective, so we should maybe wait something
998 #ifdef MOAB_HAVE_MPI
999  MB_CHK_SET_ERR( resolve_intersection_sharing(), "can't resolve intersection sharing" );
1000 #endif
1001 
1002  this->clean();
1003  return MB_SUCCESS;
1004 }

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(), findNodes(), moab::Range::front(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), gid, moab::Interface::id_from_handle(), imaskTag, moab::Range::index(), moab::Range::insert(), moab::intersect(), MAXEDGES, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, mbs1, mbs2, moab::Range::merge(), my_rank, outSet, moab::AdaptiveKDTree::parse_options(), moab::Range::pop_front(), moab::Range::rbegin(), moab::Range::rend(), rs1, rs2, setup_tgt_cell(), moab::Range::size(), srcNeighTag, moab::subtract(), 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 
)

Slow KD-tree-based mesh intersection routine (no advancing-front).

Overview:

  • Builds a KD-tree over source (mbs1) faces and, for each target (mbs2) face, queries nearby source leaves using a distance-based search around target vertices.
  • For the candidate source faces gathered from nearby KD-tree leaves, computes exact polygonal intersections in a gnomonic plane, accumulates overlap area, and creates intersection polygons/nodes in outSet via findNodes.
  • This path is intentionally simpler and potentially more expensive than the advancing-front algorithm used by intersect_meshes.

Inputs/Assumptions:

  • mbset1 (source) fully covers mbset2 (target) on the sphere.
  • Both sets contain 2D elements (triangles, quads, or generic convex polygons).
  • On-sphere intersection math is performed using gnomonic projection; tolerances are derived from maximum edge lengths on the source mesh.

High-level Steps: 1) Cache 2D entities of source (rs1) and target (rs2). Optionally filter by GRID_IMASK tag to exclude masked-out elements. 2) Precompute and tag target-edge adjacency (__tgtEdgeNeighbors) for quick access when locating/creating intersection points on target boundaries. 3) Estimate tolerances: compute maximum source edge length to derive KD-tree search tolerance and box overlap epsilon; reduce across ranks under MPI. 4) Build an AdaptiveKDTree on the source faces with spherical options (PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;). 5) For each target face:

  • Gather its vertex coordinates; compute an average edge length av_len.
  • For each target vertex, perform kd.distance_search within radius av_len to collect nearby KD-tree leaves; accumulate their contained 2D source faces into close_source_cells.
  • For each candidate source face in that range, call computeIntersectionBetweenTgtAndSrc to compute polygon intersection points and area; if area > 0, call findNodes to create nodes/polygons in outSet.
  • Track recovered area vs. the target cell area (diagnostic). 6) Under MPI, reconcile shared intersection points across process boundaries via resolve_intersection_sharing. 7) Cleanup transient state and return.

Complexity Notes:

  • Building the KD-tree is roughly O(N log N). For each target face, the search radius heuristic (av_len) aims to limit candidates; worst-case behavior can still approach quadratic if meshes overlap densely.

Key Data/Tags:

  • tgtParentTag, srcParentTag, countTag maintain provenance and counters for created intersection entities; they are (re)created in this routine.
  • __tgtEdgeNeighbors stores per-target-face edge handles to speed boundary ops.

Error handling:

  • Uses MB_CHK_ERR/MB_CHK_SET_ERR macros for MOAB ErrorCode propagation.
  • Cleans up tags and temporary state before returning on success.

Definition at line 343 of file Intx2Mesh.cpp.

344 {
345  ErrorCode rval;
346  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
347  mbs2 = mbset2;
348  outSet = outputSet;
349  MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs1, 2, rs1 ), "can't get source entities by dimension" );
350  MB_CHK_SET_ERR( mb->get_entities_by_dimension( mbs2, 2, rs2 ), "can't get target entities by dimension" );
351  // from create tags, copy relevant ones
354  if( countTag ) mb->tag_delete( countTag );
355 
356  // filter rs1 and rs2 by mask; remove everything with 0 mask
357  // get the mask tag if it exists; if not, leave it uninitialized (nullptr)
358  rval = mb->tag_get_handle( "GRID_IMASK", imaskTag );
359  if( imaskTag != nullptr && rval != MB_SUCCESS ) MB_CHK_SET_ERR( rval, "can't get GRID_IMASK tag" );
360  MB_CHK_SET_ERR( filterByMask( rs1 ), "can't filter source by mask" );
361  MB_CHK_SET_ERR( filterByMask( rs2 ), "can't filter target by mask" );
362  // create tgt edges if they do not exist yet; so when they are looked upon, they are found
363  // this is the only call that is potentially NlogN, in the whole method
365  "can't get adjacent target edges" );
366 
367  int index = 0;
368  extraNodesVec.resize( TgtEdges.size() );
369  for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, index++ )
370  {
371  std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
372  extraNodesVec[index] = nv;
373  }
374 
375  int defaultInt = -1;
376  // Now let us create the association tags to source and target parent, along with internal counters
378  &defaultInt ),
379  "can't create target parent tag" );
381  &defaultInt ),
382  "can't create source parent tag" );
384  &defaultInt ),
385  "can't create Counting tag" );
386 
387  // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
388  // them each time edges were for sure created before (tgtEdges)
389  // if we have a tag with this name, it could be of a different size, so delete it
390  rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
392  std::vector< EntityHandle > zeroh( max_edges_2, 0 );
394  MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] ),
395  "can't create tgt edge neighbors tag" );
396 
397  for( Range::iterator rit = rs2.begin(); rit != rs2.end(); ++rit )
398  {
399  EntityHandle tgtCell = *rit;
400  int num_nodes = 0;
401  MB_CHK_SET_ERR( mb->get_connectivity( tgtCell, tgtConn, num_nodes ), "can't get target connectivity" );
402  // account for padded polygons
403  while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
404  num_nodes--;
405 
406  for( int i = 0; i < num_nodes; i++ )
407  {
408  EntityHandle v[2] = { tgtConn[i],
409  tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
410  std::vector< EntityHandle > adj_entities;
411  rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
412  if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
413  zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
414  // also, even if number of edges is less than max_edges_2, they will be ignored, even if
415  // the tag is dense
416  }
417  // now set the value of the tag
418  MB_CHK_SET_ERR( mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) ),
419  "can't set edge target edge neighbors tag" );
420  }
421 
422  // find out max edge on source mesh;
423  double max_length = 0;
424  {
425  std::vector< double > coords( 3 * max_edges_1, 0.0 );
426  for( Range::iterator it = rs1.begin(); it != rs1.end(); ++it )
427  {
428  const EntityHandle* conn = nullptr;
429  int nnodes;
430  MB_CHK_SET_ERR( mb->get_connectivity( *it, conn, nnodes ), "can't get source connectivity" );
431  while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
432  nnodes--;
433  MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &coords[0] ), "can't get source coordinates" );
434  for( int j = 0; j < nnodes; j++ )
435  {
436  int next = ( j + 1 ) % nnodes;
437  double edge_length =
438  ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
439  ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
440  ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
441  if( edge_length > max_length ) max_length = edge_length;
442  }
443  }
444  max_length = std::sqrt( max_length );
445  }
446 
447  // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
448  double tolerance = 1.e-15;
449  if( max_length < 1. )
450  {
451  // basically, the sag for an arc of length max_length on a circle of radius 1
452  tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
454  tolerance = 3 * tolerance; // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
455  // be more generous, use 1.5 ~= sqrt(2.)
456 
457  if( !my_rank )
458  {
459  std::cout << " max edge length: " << max_length << " tolerance for kd tree: " << tolerance << "\n";
460  std::cout << " box overlap tolerance: " << box_error << "\n";
461  }
462  }
463 #ifdef MOAB_HAVE_MPI
464  // reduce box tolerance on every task, if needed
465  double min_box_eps = box_error;
466  if( nullptr != parcomm ) MPI_Allreduce( &box_error, &min_box_eps, 1, MPI_DOUBLE, MPI_MIN, parcomm->comm() );
467  box_error = min_box_eps;
468 #endif
469 
470  // create the kd tree on source cells, and intersect all targets in an expensive loop
471  // build a kd tree with the rs1 (source) cells
472  FileOptions kdOpts( "PLANE_SET=1;SPLITS_PER_DIR=2;SPHERICAL;RADIUS=1.0;" );
473  AdaptiveKDTree kd( mb );
474  kd.parse_options( kdOpts );
475  EntityHandle tree_root = 0;
476  MB_CHK_SET_ERR( kd.build_tree( rs1, &tree_root ), "can't build Kd-tree on source cells" );
477 
478  for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
479  {
480  EntityHandle tcell = *it;
481  // find vertex positions
482  const EntityHandle* conn = nullptr;
483  int nnodes = 0;
484  MB_CHK_SET_ERR( mb->get_connectivity( tcell, conn, nnodes ), "can't get target connectivity" );
485  // find leaves close to those positions
486  // setup_tgt_cell is called for its side effects (populates tgtConn /
487  // redCoords[] used by computeIntersectionBetweenTgtAndSrc below);
488  // the returned gnomonic-plane area is unused on this serial path.
489  (void)setup_tgt_cell( tcell, nnodes );
490  std::vector< double > positions;
491  positions.resize( nnodes * 3 );
492  MB_CHK_SET_ERR( mb->get_coords( conn, nnodes, &positions[0] ), "can't get target coordinates" );
493 
494  // distance to search will be based on average edge length
495  double av_len = 0;
496  for( int k = 0; k < nnodes; k++ )
497  {
498  int ik = ( k + 1 ) % nnodes;
499  double len1 = 0;
500  for( int j = 0; j < 3; j++ )
501  {
502  double len2 = positions[3 * k + j] - positions[3 * ik + j];
503  len1 += len2 * len2;
504  }
505  av_len += sqrt( len1 );
506  }
507  if( nnodes > 0 ) av_len /= nnodes;
508  // find leaves within a distance from each vertex of target
509  // in those leaves, collect all cells; we will try for an intx in there
510  Range close_source_cells;
511  std::vector< EntityHandle > leaves;
512  for( int i = 0; i < nnodes; i++ )
513  {
514  leaves.clear();
515  MB_CHK_SET_ERR( kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 ),
516  "can't search for leaves" );
517 
518  for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
519  {
520  Range tmp;
521  MB_CHK_SET_ERR( mb->get_entities_by_dimension( *j, 2, tmp ), "can't get entities by dimension" );
522 
523  close_source_cells.merge( tmp.begin(), tmp.end() );
524  }
525  }
526 #ifdef VERBOSE
527  if( close_source_cells.empty() )
528  {
529  std::cout << " there are no close source cells to target cell " << tcell
530  << " ht:" << mb->id_from_handle( tcell ) << "\n";
531  }
532 #endif
533  for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
534  {
535  EntityHandle startSrc = *it2;
536  double area = 0;
537  // if area is > 0 , we have intersections
538  double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
539  //
540  int nP = 0;
541  int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
542  int nsTgt, nsSrc;
543  MB_CHK_SET_ERR( computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt,
544  true ),
545  "can't compute intersection between target and source" );
546  if( area > 0 )
547  {
548  if( nP > 1 )
549  { // this will also construct triangles/polygons in the new mesh, if needed
550  MB_CHK_SET_ERR( findNodes( tcell, nnodes, startSrc, nsSrc, P, nP ), "can't find nodes" );
551 #ifdef ENABLE_DEBUG
552  std::cout << " intersect: " << " ht:" << mb->id_from_handle( tcell ) << " "
553  << " hs:" << mb->id_from_handle( startSrc ) << " g:" << global_id_ent( mb, tcell, gid )
554  << " g:" << global_id_ent( mb, startSrc, gid ) << " counting: " << counting << "\n";
555 #endif
556  }
557  // (Historical: serial path used to accumulate a per-cell
558  // recoveredArea / areaTgtCell ratio here, but neither value
559  // was ever read again. The verbose-print branch for that
560  // ratio lives only in the parallel path ~lines 893-899.)
561  }
562  }
563  }
564  // before cleaning up , we need to settle the position of the intersection points
565  // on the boundary edges
566  // this needs to be collective, so we should maybe wait something
567 #ifdef MOAB_HAVE_MPI
568  if( nullptr != parcomm )
569  {
570  MB_CHK_SET_ERR( resolve_intersection_sharing(), "can't resolve intersection sharing (correct position)" );
571  }
572 #endif
573 
574  this->clean();
575  return MB_SUCCESS;
576 }

References moab::Range::begin(), box_error, moab::AdaptiveKDTree::build_tree(), clean(), moab::Range::clear(), computeIntersectionBetweenTgtAndSrc(), counting, 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(), gid, moab::Interface::id_from_handle(), imaskTag, moab::index, moab::Interface::INTERSECT, max_edges_1, max_edges_2, MAXEDGES, mb, 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, 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(), and main().

◆ set_box_error()

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

Definition at line 132 of file Intx2Mesh.hpp.

133  {
134  box_error = berror;
135  }

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 115 of file Intx2Mesh.hpp.

116  {
117  epsilon_1 = eps;
118  epsilon_area = eps * sqrt( eps );
119  }

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 243 of file Intx2Mesh.hpp.

◆ box_error

◆ counting

◆ countTag

◆ epsilon_1

◆ epsilon_area

◆ extraNodesVec

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

◆ gid

◆ imaskTag

Tag moab::Intx2Mesh::imaskTag
protected

for coverage mesh, will store the original sender

Definition at line 220 of file Intx2Mesh.hpp.

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

◆ localEnts

Range moab::Intx2Mesh::localEnts
protected

Definition at line 247 of file Intx2Mesh.hpp.

◆ localRoot

EntityHandle moab::Intx2Mesh::localRoot
protected

Definition at line 246 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 193 of file Intx2Mesh.hpp.

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

◆ mbs2

EntityHandle moab::Intx2Mesh::mbs2
protected

◆ my_rank

unsigned int moab::Intx2Mesh::my_rank
protected

◆ neighTgtEdgeTag

◆ orgSendProcTag

Tag moab::Intx2Mesh::orgSendProcTag
protected

Definition at line 219 of file Intx2Mesh.hpp.

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

◆ outSet

◆ rs1

◆ rs2

◆ srcConn

◆ srcCoords

◆ srcCoords2D

◆ srcNeighTag

Tag moab::Intx2Mesh::srcNeighTag
protected

Definition at line 212 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

◆ srcParentTag

◆ tgtConn

◆ tgtCoords

◆ tgtCoords2D

◆ TgtEdges

◆ TgtFlagTag

Tag moab::Intx2Mesh::TgtFlagTag
protected

Definition at line 202 of file Intx2Mesh.hpp.

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

◆ tgtNeighTag

Tag moab::Intx2Mesh::tgtNeighTag
protected

Definition at line 214 of file Intx2Mesh.hpp.

Referenced by createTags(), and intersect_meshes().

◆ tgtParentTag


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