Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
Intx2MeshEdges.cpp
Go to the documentation of this file.
1 /*
2  * Intx2MeshEdges.cpp
3  *
4  * Created on: Aug. 4, 2025
5  * Author: iulian
6  */
7 
8 #include "moab/IntxMesh/Intx2MeshEdges.hpp"
9 #ifdef MOAB_HAVE_MPI
10 #include "moab/ParallelComm.hpp"
11 #endif
12 
13 namespace moab
14 {
15 
16 Intx2MeshEdges::Intx2MeshEdges( Interface* mbimpl, IntxAreaUtils::AreaMethod amethod )
17  : Intx2MeshOnSphere( mbimpl, amethod )
18 {
19 }
20 
21 Intx2MeshEdges::~Intx2MeshEdges()
22 {
23  // TODO Auto-generated destructor stub
24 }
25 
26 // local method
27 ErrorCode Intx2MeshEdges::orderSubEdges( std::vector< EntityHandle >& subEdges,
28  std::vector< EntityHandle >& VerticesSubEdges,
29  const EntityHandle* connEdge,
30  std::vector< EntityHandle >& chainVertices,
31  std::vector< int >& polygonIds,
32  Tag otherParentTag )
33 {
34  int numEdges = (int)subEdges.size();
35  if( numEdges == 1 ) return MB_SUCCESS; // nothing to do
36 
37  EntityHandle currentVertex = connEdge[0]; // start vertex
38  chainVertices.push_back( currentVertex );
39  EntityHandle endVertex = connEdge[1];
40  std::vector< EntityHandle > chain;
41  std::vector< int > markedEdge( numEdges, 0 ); // 0 means not found yet; -1 or +1 for orientation
42  // start finding the current vertex, until we close the chain; double loop, we could be smarter :)
43  for( int i = 0; i < numEdges; i++ )
44  {
45  for( int j = 0; j < numEdges; j++ )
46  {
47  if( 0 != markedEdge[j] ) continue; // do not use it anymore
48  if( VerticesSubEdges[j * 2] == currentVertex )
49  {
50  currentVertex = VerticesSubEdges[j * 2 + 1];
51  chainVertices.push_back( currentVertex );
52  chain.push_back( subEdges[j] );
53  markedEdge[j] = 1; // positive
54  break; // break the j loop
55  }
56  if( VerticesSubEdges[j * 2 + 1] == currentVertex )
57  {
58  currentVertex = VerticesSubEdges[j * 2];
59  chainVertices.push_back( currentVertex );
60  chain.push_back( subEdges[j] );
61  markedEdge[j] = -1; // reversed
62  break; // break the j loop
63  }
64  }
65  }
66  if( (int)chain.size() == numEdges && currentVertex == endVertex )
67  {
68  subEdges = chain; // reordered list, no orientation saved; maybe we should ?
69  // from chain, form the list of original polygons that contain each subedge
70  // get the parent tag of 2 intx polys connected to each edge in chain
71  for( int j = 0; j < (int)chain.size(); j++ )
72  {
73  EntityHandle sEdge = chain[j];
74  std::vector< EntityHandle > intxPolys;
75  ErrorCode rval = mb->get_adjacencies( &sEdge, 1, 2, false, intxPolys, Interface::UNION );MB_CHK_ERR( rval );
76  EntityHandle onePolygon = intxPolys[0];
77  int global_id = 0;
78  rval = mb->tag_get_data( otherParentTag, &onePolygon, 1, &global_id );MB_CHK_ERR( rval );
79  polygonIds.push_back( global_id );
80  }
81  return MB_SUCCESS;
82  }
83  else
84  return MB_FAILURE; // we did not find a chain, do not change anything
85 }
86 
87 ErrorCode Intx2MeshEdges::EdgeSplits( double areaTolerance )
88 {
89 
90  Tag parentTag, otherParentTag;
91  ErrorCode rval;
92 
93  rval = mb->tag_get_handle( "TargetParent", parentTag );MB_CHK_SET_ERR( rval, "can't get parent target tag in edge map" );
94  rval = mb->tag_get_handle( "SourceParent", otherParentTag );MB_CHK_SET_ERR( rval, "can't get parent source tag in edge map" );
95 
96  Tag fractionTag;
97  rval = mb->tag_get_handle( "EdgeRecoveryFraction", fractionTag );MB_CHK_SET_ERR( rval, "can't get tag for recovery fraction" );
98  Tag subTag;
99  rval = mb->tag_get_handle( "NumSubEnts", subTag );MB_CHK_SET_ERR( rval, "can't get tag NumSubEdges" );
100  Tag areaDiffTag;
101  rval = mb->tag_get_handle( "AreaDiff", areaDiffTag );MB_CHK_SET_ERR( rval, "can't get tag AreaDiff" );
102  Tag areaTag;
103  rval = mb->tag_get_handle( "Area", areaTag );MB_CHK_SET_ERR( rval, "can't get tag Area" );
104 
105  Range parentCells;
106  rval = mb->get_entities_by_dimension( mbs2, 2, parentCells );MB_CHK_SET_ERR( rval, "can't get parent cells" );
107  Range initialEdges;
108  // we will assume the edges are already included in the target mesh
109  // we do not want to create them
110  rval = mb->get_adjacencies( parentCells, 1, false, initialEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get parent edges" );
111  std::cout << " number of original input edges:" << initialEdges.size() << "\n";
112 
113  // get all polygons in the intx set
114  Range cells;
115  rval = mb->get_entities_by_dimension( outSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get intersection cells" );
116  // create all edges adjacent to the cells in the intx set
117  // some might be original edges from edge or target meshes
118  Range intxEdges;
119  rval = mb->get_adjacencies( cells, 1, true, intxEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't create adjacent intx edges " );
120  std::cout << " number of intx edges:" << intxEdges.size() << "\n";
121 
122  // first, identify input polygons that are recovered fully
123  // get global ids of the initial cells
124  std::vector< int > parentGids( parentCells.size() );
125  Tag gidTag = mb->globalId_tag();
126  rval = mb->tag_get_data( gidTag, parentCells, &parentGids[0] );MB_CHK_SET_ERR( rval, "can't get parent global ids" );
127  std::map< int, double > initAreas; // get areas of those initial cells
128  int i = 0;
129  std::vector< double > coords( 3 * 20 ); // enough vertices, at most 30, good for intx polygons too
130  int num_nodes = 0; // num nodes in cells
131  const EntityHandle* verts;
132  IntxAreaUtils areaAdaptor( IntxAreaUtils::lHuiller ); // GaussQuadrature , lHuiller
133  std::map< int, double > recoveredAreas; // get areas of from intersection cells
134  std::map< int, EntityHandle > mapFromGIDToParent;
135  std::map< int, std::vector< EntityHandle > > mapFromParentGIDToIntxCells;
136  for( Range::iterator it = parentCells.begin(); it != parentCells.end(); ++it, i++ )
137  {
138  EntityHandle parentCell = *it;
139  rval = mb->get_connectivity( parentCell, verts, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity of parent cell" );
140  // get coordinates
141  rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coords of parent cell" );
142 
143  double area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, 1. );
144  rval = mb->tag_set_data( areaTag, &parentCell, 1, &area );MB_CHK_SET_ERR( rval, "can't set area Tag on parent cell" );
145  int parentID = parentGids[i];
146  initAreas[parentID] = area;
147  recoveredAreas[parentID] = 0.;
148  mapFromGIDToParent[parentID] = parentCell;
149  mapFromParentGIDToIntxCells[parentID]; // just initialize it with empty vector
150  }
151 
152  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
153  {
154  EntityHandle cell = *it;
155  rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_SET_ERR( rval, "can't get connectivity of intx cell" );
156  // get coordinates
157  rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coods of intx cell" );
158  double intx_area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, 1. );
159  int parentID = 0;
160  rval = mb->tag_get_data( parentTag, &cell, 1, &parentID );MB_CHK_SET_ERR( rval, "can't get parent Tag" );
161  recoveredAreas[parentID] += intx_area;
162  mapFromParentGIDToIntxCells[parentID].push_back( cell );
163  }
164  int recovered = 0, notRecovered = 0;
165  //Range recoveredCells;
166  for( size_t j = 0; j < parentGids.size(); j++ )
167  {
168  int parentID = parentGids[j];
169  EntityHandle parentCell = parentCells[j];
170  double areaDiff = initAreas[parentID] - recoveredAreas[parentID];
171  rval = mb->tag_set_data( areaDiffTag, &parentCell, 1, &areaDiff );MB_CHK_SET_ERR( rval, "can't set diff area Tag" );
172  double numIntxcells = static_cast< double >( mapFromParentGIDToIntxCells[parentID].size() );
173  rval = mb->tag_set_data( subTag, &parentCell, 1, &numIntxcells );MB_CHK_SET_ERR( rval, "can't set num sub ents Tag" );
174  if( fabs( areaDiff ) < areaTolerance )
175  {
176  recovered++;
177  recoveredCells.insert( parentCells[j] ); // should we use a std::vector, that will be ordered already ?
178  }
179  else
180  notRecovered++;
181  }
182  std::cout << "recovered initial cells: " << recovered << " vs:" << notRecovered << " not recovered \n";
183  // initial edges that should be decomposable from intx edges
184  Range recoverableEdges;
185  rval = mb->get_adjacencies( recoveredCells, 1, false, recoverableEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get recoverable edges" );
186 
187  // now recover each edge, looking at intx polygons forming one of the adjacent cells, that is in initial recoverable set
188  int recoveredEdges = 0;
189  int unrecovered = 0;
190  int identity_edges = 0; // edges that are formed by one intx edge, itself, the original
191  std::map< EntityHandle, std::vector< EntityHandle > > mapEdges;
192  for( Range::iterator eit = recoverableEdges.begin(); eit != recoverableEdges.end(); ++eit )
193  {
194  EntityHandle initialEdge = *eit;
195  // if this edge is among intxEdges, we are done
196  int index = intxEdges.index( initialEdge );
197  if( index >= 0 )
198  {
199  mapEdges[initialEdge].push_back( initialEdge );
200  identity_edges++;
201  recoveredEdges++;
202  double numSubEdges = 1.;
203  rval = mb->tag_set_data( subTag, &initialEdge, 1, &numSubEdges );MB_CHK_SET_ERR( rval, "can't set num sub ents Tag" );
204  // get vertices and intx poly attached to it
205  int nve = 0;
206  const EntityHandle* connCell;
207  rval = mb->get_connectivity( initialEdge, connCell, nve );MB_CHK_SET_ERR( rval, "can't get connectivity of parent cell" );
208  edgeVertices[initialEdge].push_back( connCell[0] );
209  edgeVertices[initialEdge].push_back( connCell[1] );
210  // find cells in intx set adjacent to it, and get the other tag parent
211  Range adjPolys;
212  rval = mb->get_adjacencies( &initialEdge, 1, 2, false, adjPolys, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adj polys" );
213  adjPolys = subtract( adjPolys, parentCells );
214  if( adjPolys.empty() ) MB_CHK_SET_ERR( MB_FAILURE, "no adjacent intx cells" );
215  EntityHandle intxPoly = adjPolys[0];
216  // get its parent tag
217  int gid;
218  rval = mb->tag_get_data( otherParentTag, &intxPoly, 1, &gid );MB_CHK_SET_ERR( rval, "can't get global id from other tag" );
219  edgePolygons[initialEdge].push_back( gid );
220  continue; // no need to sweat it anymore
221  }
222  // get adjacent polygons; if there is an adj polygon in intx cells, we are done; if not, get one in the initial recoveredCells range
223  Range initialAdjCells;
224  rval = mb->get_adjacencies( &initialEdge, 1, 2, false, initialAdjCells, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent initial cells" );
225  if( initialAdjCells.size() == 0 ) MB_CHK_SET_ERR( MB_FAILURE, "no adjacent initial cells" );
226  // intersect with recoveredCells
227  Range adjRecoveredInitialCell = intersect( initialAdjCells, recoveredCells );
228  if( adjRecoveredInitialCell.size() == 0 )
229  MB_CHK_SET_ERR( MB_FAILURE, "no adjacent initial cells that are recovered" );
230  // get all edges adjacent to intersection polygons that form one of the recovered initial cell
231  EntityHandle recoveredCell = adjRecoveredInitialCell[0]; // first one
232  int nv = 0;
233  const EntityHandle* connCell;
234  rval = mb->get_connectivity( recoveredCell, connCell, nv );MB_CHK_SET_ERR( rval, "can't get connectivity of parent cell" );
235  int cellGlobalId = 0;
236  rval = mb->tag_get_data( gidTag, &recoveredCell, 1, &cellGlobalId );MB_CHK_SET_ERR( rval, "can't get id of parent cell" );
237  // list of intx polygons in it:
238  std::vector< EntityHandle >& listIntxCells = mapFromParentGIDToIntxCells[cellGlobalId];
239  std::vector< EntityHandle > candidateEdges;
240  rval =
241  mb->get_adjacencies( &listIntxCells[0], listIntxCells.size(), 1, false, candidateEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adj edges" );
242  // add all edges that have both vertices on the original edge
243  CartVect verticesOriginal[2];
244  const EntityHandle* connEdge;
245  int numVerts = 0;
246  rval = mb->get_connectivity( initialEdge, connEdge, numVerts );MB_CHK_SET_ERR( rval, "can't get connectivity of initial edge" );
247 
248  rval = mb->get_coords( connEdge, numVerts, verticesOriginal[0].array() );MB_CHK_SET_ERR( rval, "can't get coordinates of vertices of initial edge" );
249  // get gnomonic plane for the start of the edge
250  int gnomonicPlane = 0;
251  IntxUtils::decide_gnomonic_plane( verticesOriginal[0], gnomonicPlane );
252  double coords2D[6]; // coords in gnomonic plane, for interior point decision
253  rval = IntxUtils::gnomonic_projection( verticesOriginal[0], 1.0, gnomonicPlane, coords2D[0], coords2D[1] );MB_CHK_SET_ERR( rval, "can't get gnomonic coords" );
254  rval = IntxUtils::gnomonic_projection( verticesOriginal[1], 1.0, gnomonicPlane, coords2D[2], coords2D[3] );MB_CHK_SET_ERR( rval, "can't get gnomonic coords" );
255 
256  double edgeLength =
257  angle_robust( verticesOriginal[0], verticesOriginal[1] ); // distance on sphere of radius 1 is the angle
258  // loop now over candidateEdges
259  double recoveredLength = 0.;
260  std::vector< EntityHandle > VerticesSubEdges; // push all vertices here, so we can order the subedges later
261  for( size_t i = 0; i < candidateEdges.size(); i++ )
262  {
263  EntityHandle subedge = candidateEdges[i];
264  // get its vertex coordinates:
265  const EntityHandle* connEdge2;
266  rval = mb->get_connectivity( subedge, connEdge2, numVerts );MB_CHK_SET_ERR( rval, "can't get connectivity of candidate edge" );
267  CartVect verticesSubEdge[2];
268  rval = mb->get_coords( connEdge2, numVerts, verticesSubEdge[0].array() );MB_CHK_SET_ERR( rval, "can't get coordinates of vertices of initial edge" );
269  // decide if they are on the initial edge (form an angle)
270  bool onEdge = true;
271  for( int j = 0; j < 2 && onEdge; j++ )
272  {
273  rval =
274  IntxUtils::gnomonic_projection( verticesSubEdge[j], 1.0, gnomonicPlane, coords2D[4], coords2D[5] );MB_CHK_SET_ERR( rval, "can't get gnomonic coords of subedge" );
275  // area of triangle in gnomonic plane should be 0
276  double areaAbs = fabs( IntxUtils::area2D( &coords2D[0], &coords2D[2], &coords2D[4] ) );
277 
278  if( areaAbs > 1.e-14 ) onEdge = false;
279  }
280  if( onEdge )
281  {
282  double subEdgeLen = angle_robust( verticesSubEdge[0], verticesSubEdge[1] );
283  recoveredLength += subEdgeLen;
284  mapEdges[initialEdge].push_back( subedge );
285  VerticesSubEdges.push_back( connEdge2[0] );
286  VerticesSubEdges.push_back( connEdge2[1] );
287  }
288  }
289  double fraction = recoveredLength / edgeLength;
290  rval = mb->tag_set_data( fractionTag, &initialEdge, 1, &fraction );MB_CHK_SET_ERR( rval, "can't set fraction on initial edge" );
291  double numSubEdge = double( mapEdges[initialEdge].size() );
292  rval = mb->tag_set_data( subTag, &initialEdge, 1, &numSubEdge );MB_CHK_SET_ERR( rval, "can't set number of subedges on initial edge" );
293  // now , reorder the subedges to create a chain along the original edge
294  // if we cannot form a chain, it means we have a problem
295  // order subedges on the original edge, and find out their orientation
296  // set also the parent tag on the edge, either source or target parent
297  // in some cases, the parents can be both
298  rval = orderSubEdges( mapEdges[initialEdge], VerticesSubEdges, connEdge, edgeVertices[initialEdge],
299  edgePolygons[initialEdge], otherParentTag );
300  if( fabs( edgeLength - recoveredLength ) < 1.e-10 || rval == MB_SUCCESS )
301  recoveredEdges++;
302  else
303  {
304  mb->list_entity( initialEdge );
305  std::vector< EntityHandle >& listSubEdges = mapEdges[initialEdge];
306  double newCheckLength = 0;
307  CartVect vertices[2];
308  const EntityHandle* conn2;
309  int numVerts2 = 0;
310  rval = mb->get_connectivity( initialEdge, conn2, numVerts2 );MB_CHK_SET_ERR( rval, "can't get connectivity of initial edge " );
311  rval = mb->get_coords( conn2, numVerts2, vertices[0].array() );MB_CHK_SET_ERR( rval, "can't get coordinates of vertices of initial edge" );
312  double length2 = angle_robust( vertices[0], vertices[1] );
313  std::cout << std::setprecision( 14 );
314  std::cout << " initial edge:" << mb->id_from_handle( initialEdge ) << " v: " << conn2[0] << ", " << conn2[1]
315  << " len: " << length2 << "\n";
316  for( size_t j = 0; j < listSubEdges.size(); j++ )
317  {
318  EntityHandle subEdge = listSubEdges[j];
319  rval = mb->get_connectivity( subEdge, conn2, numVerts2 );MB_CHK_SET_ERR( rval, "can't get connectivity of subedge " );
320  rval = mb->get_coords( conn2, numVerts2, vertices[0].array() );MB_CHK_SET_ERR( rval, "can't get coordinates of vertices of subedge" );
321  length2 = angle_robust( vertices[0], vertices[1] );
322  newCheckLength += length2;
323  std::cout << " sub edge:" << mb->id_from_handle( subEdge ) << " v: " << conn2[0] << ", " << conn2[1]
324  << " len: " << length2 << "\n";
325  }
326  std::cout << " edge length:" << edgeLength << " newCheckLength:" << newCheckLength
327  << " diff:" << edgeLength - newCheckLength << " fraction:" << fraction
328  << " subedges:" << numSubEdge << "\n";
329  unrecovered++;
330  std::cout << std::setprecision( 7 );
331  }
332  }
333 
334  std::cout << " recoveredEdges:" << recoveredEdges << " identity edges:" << identity_edges
335  << " unrecovered edges:" << unrecovered << "\n";
336  return MB_SUCCESS;
337 }
338 
339 #ifdef MOAB_HAVE_PNETCDF
340 #ifdef MOAB_HAVE_MPI
341 ErrorCode Intx2MeshEdges::write_edge_map_parallel( const char* fileName )
342 {
343  // open for writing the pnetcdf nc file
344  int ncid; // file id
345  int ncell_dimid, max_edge_dimid, max_sub_edge_dimid, max_sub_edgeP1_dimid;
346  int retval; // return val for nc
347  //Below is an example code fragment that sets the file header hint to 1MB and pass it to PnetCDF when creating a file.
348 
349  if( ( retval = ncmpi_create( parcomm->comm(), fileName, NC_CLOBBER | NC_64BIT_DATA, MPI_INFO_NULL, &ncid ) ) )
350  ERR( retval );
351 
352  int num_cells_local;
353  Range polys;
354  ErrorCode rval = mb->get_entities_by_dimension( mbs2, 2, polys );MB_CHK_SET_ERR( rval, "Failed to get polygons" );
355  num_cells_local = (int)polys.size();
356 
357  int num_cells = 0;
358  // do mpi reduce all, to find out the total num_cells;
359  MPI_Allreduce( &num_cells_local, &num_cells, 1, MPI_INTEGER, MPI_SUM, parcomm->comm() );
360 
361  if( ( retval = ncmpi_def_dim( ncid, "num_cells", num_cells, &ncell_dimid ) ) ) ERR( retval );
362 
363  Tag gid = mb->globalId_tag();
364  std::vector< int > global_ids_polys( num_cells_local );
365  rval = mb->tag_get_data( gid, polys, global_ids_polys.data() );MB_CHK_SET_ERR( rval, "Failed to get ids of cells" );
366  Range localGidCells;
367  std::copy( global_ids_polys.rbegin(), global_ids_polys.rend(), range_inserter( localGidCells ) );
368 
369  // find max_edges and max subedges
370  int max_edge = -1, max_sub_edge = -1, max_subedge1 = -1;
371 
372  for( auto it = polys.begin(); it != polys.end(); ++it )
373  {
374  EntityHandle polygon = *it;
375  const EntityHandle* conn = nullptr;
376  int nv;
377  rval = mb->get_connectivity( polygon, conn, nv );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
378  if( max_edge < nv ) max_edge = nv;
379  }
380  int max_edgeg = 0;
381  // do mpi reduce all, to find out the max number of edges;
382  MPI_Allreduce( &max_edge, &max_edgeg, 1, MPI_INTEGER, MPI_MAX, parcomm->comm() );
383  if( ( retval = ncmpi_def_dim( ncid, "max_edges", max_edgeg, &max_edge_dimid ) ) ) ERR( retval );
384  max_edge = max_edgeg;
385  for( auto mapit = edgePolygons.begin(); mapit != edgePolygons.end(); ++mapit )
386  {
387  int nsb = (int)mapit->second.size();
388  if( max_sub_edge < nsb ) max_sub_edge = nsb;
389  }
390  int max_sub_edgeg = 0;
391  // do mpi reduce all, to find out the max number of subedges;
392  MPI_Allreduce( &max_sub_edge, &max_sub_edgeg, 1, MPI_INTEGER, MPI_MAX, parcomm->comm() );
393  if( ( retval = ncmpi_def_dim( ncid, "max_sub_edges", max_sub_edgeg, &max_sub_edge_dimid ) ) ) ERR( retval );
394  max_sub_edge = max_sub_edgeg;
395  max_subedge1 = max_sub_edge + 1;
396  if( ( retval = ncmpi_def_dim( ncid, "max_sub_edges1", max_subedge1, &max_sub_edgeP1_dimid ) ) ) ERR( retval );
397 
398  int dimids_nbs[2];
399 
400  dimids_nbs[0] = ncell_dimid;
401  dimids_nbs[1] = max_edge_dimid;
402  int varid_nsub;
403  if( ( retval = ncmpi_def_var( ncid, "nb_sub_edge", NC_INT, 2, dimids_nbs, &varid_nsub ) ) ) ERR( retval );
404 
405  int dimids_cell_assoc[3];
406  dimids_cell_assoc[0] = ncell_dimid;
407  dimids_cell_assoc[1] = max_edge_dimid;
408  dimids_cell_assoc[2] = max_sub_edge_dimid;
409  int varid_cell_assoc;
410  if( ( retval = ncmpi_def_var( ncid, "cells_assoc", NC_INT, 3, dimids_cell_assoc, &varid_cell_assoc ) ) )
411  ERR( retval );
412 
413  dimids_cell_assoc[2] = max_sub_edgeP1_dimid;
414  int varid_lat, varid_lon;
415  if( ( retval = ncmpi_def_var( ncid, "lat_sub_edge", NC_DOUBLE, 3, dimids_cell_assoc, &varid_lat ) ) ) ERR( retval );
416  if( ( retval = ncmpi_def_var( ncid, "lon_sub_edge", NC_DOUBLE, 3, dimids_cell_assoc, &varid_lon ) ) ) ERR( retval );
417 
418  ncmpi_enddef( ncid );
419 
420  // need to find out the owned local cells on this task
421  // write only those owned local cells, at the correct location, given by their global id
422  // we are repeating the shared edges anyway;
423  /// so basically, all local cells are owned; no need to worry, just write them at the correct location in the file, based on
424  // their global id
425  std::vector< int > nb_sub_edge_per_edge( num_cells_local * max_edge, -9999 );
426  std::vector< int > cells_assoc_per_edge( num_cells_local * max_edge * max_sub_edge, -9999 );
427 
428  std::vector< double > latvals( num_cells_local * max_edge * max_subedge1, -9999 );
429  std::vector< double > lonvals( num_cells_local * max_edge * max_subedge1, -9999 );
430 
431  for( auto it = recoveredCells.begin(); it != recoveredCells.end(); ++it )
432  {
433  EntityHandle polygon = *it;
434  int gidPoly = 0;
435  rval = mb->tag_get_data( gid, &polygon, 1, &gidPoly );MB_CHK_SET_ERR( rval, "Failed to get id of poly" );
436  int indexGidPoly = localGidCells.index( gidPoly );
437 
438  const EntityHandle* conn = nullptr;
439  int nv;
440  rval = mb->get_connectivity( polygon, conn, nv );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
441  for( int i = 0; i < nv; i++ )
442  {
443  EntityHandle v[2];
444  v[0] = conn[i];
445  v[1] = conn[( i + 1 ) % nv];
446  Range edges;
447  rval = mb->get_adjacencies( v, 2, 1, false, edges, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "Failed to get edge" );
448  EntityHandle edge = edges[0];
449  // reverse or not?
450  std::vector< int > poly_assoc = edgePolygons[edge];
451  nb_sub_edge_per_edge[indexGidPoly * max_edge + i] = (int)poly_assoc.size();
452  std::vector< EntityHandle > vertexEdges = edgeVertices[edge];
453  std::vector< CartVect > coords( vertexEdges.size() );
454  rval = mb->get_coords( &vertexEdges[0], vertexEdges.size(), &( coords[0][0] ) );MB_CHK_SET_ERR( rval, "can't get coordinates" );
455  // convert to lat/lon
456  std::vector< double > latv( vertexEdges.size() ), lonv( vertexEdges.size() );
457  for( int j = 0; j < (int)vertexEdges.size(); j++ )
458  {
459  IntxUtils::SphereCoords sph1 = IntxUtils::cart_to_spherical( coords[j] );
460  lonv[j] = sph1.lon;
461  latv[j] = sph1.lat;
462  }
463  // reversed edge or not?
464  const EntityHandle* edgeconn = nullptr;
465  int nve;
466  rval = mb->get_connectivity( edge, edgeconn, nve );MB_CHK_SET_ERR( rval, "Failed to get edge connectivity" );
467  bool reverse = false;
468  if( v[0] == edgeconn[1] ) reverse = true;
469  if( reverse )
470  {
471  // fill the arrays in reverse order
472  int sizep = (int)poly_assoc.size();
473  for( int j = 0; j < sizep; j++ )
474  {
475  cells_assoc_per_edge[indexGidPoly * max_edge * max_sub_edge + i * max_sub_edge + j] =
476  poly_assoc[sizep - 1 - j];
477  }
478  sizep = (int)vertexEdges.size();
479  for( int j = 0; j < sizep; j++ )
480  {
481  latvals[indexGidPoly * max_edge * max_subedge1 + i * max_subedge1 + j] = latv[sizep - 1 - j];
482  lonvals[indexGidPoly * max_edge * max_subedge1 + i * max_subedge1 + j] = lonv[sizep - 1 - j];
483  }
484  }
485  else
486  {
487  // max_sub_edges
488  for( int j = 0; j < (int)poly_assoc.size(); j++ )
489  {
490  cells_assoc_per_edge[indexGidPoly * max_edge * max_sub_edge + i * max_sub_edge + j] = poly_assoc[j];
491  }
492  for( int j = 0; j < (int)vertexEdges.size(); j++ )
493  {
494  latvals[indexGidPoly * max_edge * max_subedge1 + i * max_subedge1 + j] = latv[j];
495  lonvals[indexGidPoly * max_edge * max_subedge1 + i * max_subedge1 + j] = lonv[j];
496  }
497  }
498  }
499  }
500 
501  size_t idxReq = 0;
502  size_t nb_writes = localGidCells.psize();
503  std::vector< int > requests( nb_writes * 4 );
504  std::vector< int > statuss( nb_writes * 4 );
505 
506  size_t indexInArray1 = 0;
507  size_t indexInArray2 = 0;
508  size_t indexInArray3 = 0;
509 
510  for( Range::pair_iterator pair_iter = localGidCells.pair_begin(); pair_iter != localGidCells.pair_end();
511  ++pair_iter ) // these will be the size of nb_writes
512  {
513  EntityHandle starth = pair_iter->first;
514  EntityHandle endh = pair_iter->second;
515  MPI_Offset write_start[3], write_count[3]; // max dimension
516  write_start[0] = static_cast< MPI_Offset >( starth - 1 );
517  write_start[1] = static_cast< MPI_Offset >( 0 );
518  write_start[2] = static_cast< MPI_Offset >( 0 );
519  write_count[0] = static_cast< MPI_Offset >( endh - starth + 1 );
520  write_count[1] = static_cast< MPI_Offset >( max_edge );
521  write_count[2] = static_cast< MPI_Offset >( max_sub_edge );
522 
523  if( ( retval = ncmpi_iput_vara_int( ncid, varid_nsub, write_start, write_count, // first 2 are used
524  &nb_sub_edge_per_edge[indexInArray1], &requests[idxReq++] ) ) )
525  ERR( retval );
526  indexInArray1 += ( endh - starth + 1 ) * max_edge;
527 
528  if( ( retval = ncmpi_iput_vara_int( ncid, varid_cell_assoc, write_start, write_count, // all 3 are used
529  &cells_assoc_per_edge[indexInArray2], &requests[idxReq++] ) ) )
530  ERR( retval );
531  indexInArray2 += ( endh - starth + 1 ) * max_edge * max_sub_edge;
532 
533  write_count[2] = max_subedge1;
534 
535  if( ( retval = ncmpi_iput_vara_double( ncid, varid_lat, write_start, write_count, // all 3 are used
536  &latvals[indexInArray3], &requests[idxReq++] ) ) )
537  ERR( retval );
538 
539  if( ( retval = ncmpi_iput_vara_double( ncid, varid_lon, write_start, write_count, // all 3 are used
540  &lonvals[indexInArray3], &requests[idxReq++] ) ) )
541  ERR( retval );
542 
543  indexInArray3 += ( endh - starth + 1 ) * max_edge * max_subedge1;
544  }
545  // Wait outside the loop
546  if( ( retval = ncmpi_wait_all( ncid, requests.size(), &requests[0], &statuss[0] ) ) ) ERR( retval );
547  if( ( retval = ncmpi_close( ncid ) ) ) ERR( retval );
548  return MB_SUCCESS;
549 }
550 #endif
551 #endif
552 
553 #ifdef MOAB_HAVE_NETCDF
554 ErrorCode Intx2MeshEdges::write_edge_map( const char* filename )
555 /*Interface * mb, EntityHandle sf1,
556  std::map<EntityHandle, std::vector<EntityHandle>> & edgeVertices,
557  std::map<EntityHandle, std::vector<int>> & edgePolygons,
558  moab::Range & recoveredCells)*/
559 {
560  // open for writing the nc file
561  int ncid; // file id
562  int ncell_dimid, max_edge_dimid, max_sub_edge_dimid, max_sub_edgeP1_dimid;
563  int retval; // return val for nc
564  if( ( retval = nc_create( filename, NC_CLASSIC_MODEL | NC_CLOBBER, &ncid ) ) ) ERRN( retval );
565  int num_cells;
566  Range polys;
567  ErrorCode rval = mb->get_entities_by_dimension( mbs2, 2, polys );MB_CHK_SET_ERR( rval, "Failed to get polygons" );
568  num_cells = (int)polys.size();
569 
570  if( ( retval = nc_def_dim( ncid, "num_cells", num_cells, &ncell_dimid ) ) ) ERRN( retval );
571 
572  Tag gid = mb->globalId_tag();
573  // find max_edges and max subedges
574  int max_edge = -1, max_sub_edge = -1, max_subedge1 = -1;
575 
576  for( auto it = polys.begin(); it != polys.end(); ++it )
577  {
578  EntityHandle polygon = *it;
579  const EntityHandle* conn = nullptr;
580  int nv;
581  rval = mb->get_connectivity( polygon, conn, nv );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
582  if( max_edge < nv ) max_edge = nv;
583  }
584  if( ( retval = nc_def_dim( ncid, "max_edges", max_edge, &max_edge_dimid ) ) ) ERRN( retval );
585 
586  for( auto mapit = edgePolygons.begin(); mapit != edgePolygons.end(); ++mapit )
587  {
588  int nsb = (int)mapit->second.size();
589  if( max_sub_edge < nsb ) max_sub_edge = nsb;
590  }
591  if( ( retval = nc_def_dim( ncid, "max_sub_edges", max_sub_edge, &max_sub_edge_dimid ) ) ) ERRN( retval );
592  max_subedge1 = max_sub_edge + 1;
593  if( ( retval = nc_def_dim( ncid, "max_sub_edges1", max_subedge1, &max_sub_edgeP1_dimid ) ) ) ERRN( retval );
594 
595  int dimids_nbs[2];
596 
597  dimids_nbs[0] = ncell_dimid;
598  dimids_nbs[1] = max_edge_dimid;
599  int varid_nsub;
600  if( ( retval = nc_def_var( ncid, "nb_sub_edge", NC_INT, 2, dimids_nbs, &varid_nsub ) ) ) ERRN( retval );
601 
602  int dimids_cell_assoc[3];
603  dimids_cell_assoc[0] = ncell_dimid;
604  dimids_cell_assoc[1] = max_edge_dimid;
605  dimids_cell_assoc[2] = max_sub_edge_dimid;
606  int varid_cell_assoc;
607  if( ( retval = nc_def_var( ncid, "cells_assoc", NC_INT, 3, dimids_cell_assoc, &varid_cell_assoc ) ) )
608  ERRN( retval );
609 
610  dimids_cell_assoc[2] = max_sub_edgeP1_dimid;
611  int varid_lat, varid_lon;
612  if( ( retval = nc_def_var( ncid, "lat_sub_edge", NC_DOUBLE, 3, dimids_cell_assoc, &varid_lat ) ) ) ERRN( retval );
613  if( ( retval = nc_def_var( ncid, "lon_sub_edge", NC_DOUBLE, 3, dimids_cell_assoc, &varid_lon ) ) ) ERRN( retval );
614 
615  nc_enddef( ncid );
616 
617  std::vector< int > nb_sub_edge_per_edge( num_cells * max_edge, -9999 );
618  std::vector< int > cells_assoc_per_edge( num_cells * max_edge * max_sub_edge, -9999 );
619 
620  std::vector< double > latvals( num_cells * max_edge * max_subedge1, -9999 );
621  std::vector< double > lonvals( num_cells * max_edge * max_subedge1, -9999 );
622 
623  for( auto it = recoveredCells.begin(); it != recoveredCells.end(); ++it )
624  {
625  EntityHandle polygon = *it;
626  int gidPoly = 0;
627  rval = mb->tag_get_data( gid, &polygon, 1, &gidPoly );MB_CHK_SET_ERR( rval, "Failed to get id of poly" );
628  const EntityHandle* conn = nullptr;
629  int nv;
630  rval = mb->get_connectivity( polygon, conn, nv );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
631  for( int i = 0; i < nv; i++ )
632  {
633  EntityHandle v[2];
634  v[0] = conn[i];
635  v[1] = conn[( i + 1 ) % nv];
636  Range edges;
637  rval = mb->get_adjacencies( v, 2, 1, false, edges, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "Failed to get edge" );
638  EntityHandle edge = edges[0];
639  // reverse or not?
640  std::vector< int > poly_assoc = edgePolygons[edge];
641  nb_sub_edge_per_edge[( gidPoly - 1 ) * max_edge + i] = (int)poly_assoc.size();
642  std::vector< EntityHandle > vertexEdges = edgeVertices[edge];
643  std::vector< CartVect > coords( vertexEdges.size() );
644  rval = mb->get_coords( &vertexEdges[0], vertexEdges.size(), &( coords[0][0] ) );MB_CHK_SET_ERR( rval, "can't get coordinates" );
645  // convert to lat/lon
646  std::vector< double > latv( vertexEdges.size() ), lonv( vertexEdges.size() );
647  for( int j = 0; j < (int)vertexEdges.size(); j++ )
648  {
649  IntxUtils::SphereCoords sph1 = IntxUtils::cart_to_spherical( coords[j] );
650  lonv[j] = sph1.lon;
651  latv[j] = sph1.lat;
652  }
653  // reversed edge or not?
654  const EntityHandle* edgeconn = nullptr;
655  int nve;
656  rval = mb->get_connectivity( edge, edgeconn, nve );MB_CHK_SET_ERR( rval, "Failed to get edge connectivity" );
657  bool reverse = false;
658  if( v[0] == edgeconn[1] ) reverse = true;
659  if( reverse )
660  {
661  // fill the arrays in reverse order
662  int sizep = (int)poly_assoc.size();
663  for( int j = 0; j < sizep; j++ )
664  {
665  cells_assoc_per_edge[( gidPoly - 1 ) * max_edge * max_sub_edge + i * max_sub_edge + j] =
666  poly_assoc[sizep - 1 - j];
667  }
668  sizep = (int)vertexEdges.size();
669  for( int j = 0; j < sizep; j++ )
670  {
671  latvals[( gidPoly - 1 ) * max_edge * max_subedge1 + i * max_subedge1 + j] = latv[sizep - 1 - j];
672  lonvals[( gidPoly - 1 ) * max_edge * max_subedge1 + i * max_subedge1 + j] = lonv[sizep - 1 - j];
673  }
674  }
675  else
676  {
677  // max_sub_edges
678  for( int j = 0; j < (int)poly_assoc.size(); j++ )
679  {
680  cells_assoc_per_edge[( gidPoly - 1 ) * max_edge * max_sub_edge + i * max_sub_edge + j] =
681  poly_assoc[j];
682  }
683  for( int j = 0; j < (int)vertexEdges.size(); j++ )
684  {
685  latvals[( gidPoly - 1 ) * max_edge * max_subedge1 + i * max_subedge1 + j] = latv[j];
686  lonvals[( gidPoly - 1 ) * max_edge * max_subedge1 + i * max_subedge1 + j] = lonv[j];
687  }
688  }
689  }
690  }
691 
692  if( ( retval = nc_put_var_int( ncid, varid_nsub, &nb_sub_edge_per_edge[0] ) ) ) ERRN( retval );
693 
694  if( ( retval = nc_put_var_int( ncid, varid_cell_assoc, &cells_assoc_per_edge[0] ) ) ) ERRN( retval );
695 
696  if( ( retval = nc_put_var_double( ncid, varid_lat, &latvals[0] ) ) ) ERRN( retval );
697 
698  if( ( retval = nc_put_var_double( ncid, varid_lon, &lonvals[0] ) ) ) ERRN( retval );
699 
700  if( ( retval = nc_close( ncid ) ) ) ERRN( retval );
701  return MB_SUCCESS;
702 }
703 #endif
704 
705 } // namespace moab