Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
mbedgesplit.cpp
Go to the documentation of this file.
1 /*
2  * edge_maps_test.cpp
3  * will take source and target spherical meshes, and create edge maps using intersection file
4  * between them
5  * First, every edge in source and target will be decomposed in edges that are adjacent to
6  * the intersection polygons
7  *
8  */
9 #include <iostream>
10 #include <sstream>
11 #include <ctime>
12 #include <cstdlib>
13 #include <cstdio>
14 #include <cstring>
15 #include "moab/Core.hpp"
16 #ifdef MOAB_HAVE_MPI
17 #include "moab/ParallelComm.hpp"
18 #endif
19 #include "moab/IntxMesh/Intx2MeshEdges.hpp"
20 #include "moab/IntxMesh/IntxUtils.hpp"
21 #include "moab/ProgOptions.hpp"
22 #include <cmath>
23 
24 using namespace moab;
25 
26 int main( int argc, char* argv[] )
27 {
28 
29  std::string sourceFile, targetFile, intersectionFile, mapEdgeTargetFile;
30  //sourceFile =
31  // "../sandbox/MeshFiles/e3sm/edge_maps/source_1.h5m"; // it also has data associated to edges
32  //targetFile =
33  // "../sandbox/MeshFiles/e3sm/edge_maps/target_1.h5m"; //
34  intersectionFile = "intx_edges.h5m";
35  mapEdgeTargetFile = "target_edge.nc";
36 
37  ProgOptions opts;
38  opts.addOpt< std::string >( "source,s", "first mesh filename (source)", &sourceFile );
39  opts.addOpt< std::string >( "target,t", "second mesh filename (target)", &targetFile );
40  opts.addOpt< std::string >( "intersectionFile,i", "output intersection file", &intersectionFile );
41  opts.addOpt< std::string >( "edgeTarget,p", "output map edge target file", &mapEdgeTargetFile );
42 
43  double R = 1.; // input
44  double epsrel = 1.e-12;
45  double boxeps = 1.e-4;
46  double areaTolerance = 5.e-12;
47  intersectionFile = "intx.h5m";
48  opts.addOpt< double >( "radius,R", "radius for model intx", &R );
49  opts.addOpt< double >( "epsilon,e", "relative error in intx", &epsrel );
50  opts.addOpt< double >( "boxerror,b", "relative error for box boundaries", &boxeps );
51  opts.addOpt< double >( "areaTol,a", "area recovery tolerance", &areaTolerance );
52 
53  opts.addOpt< void >( "writeFiles,w", "write files of interest" );
54  opts.addOpt< void >( "kdtreeOption,k", "use kd tree for intersection" );
55 
56  opts.parseCommandLine( argc, argv );
57 
58  bool write_files_rank = opts.numOptSet( "writeFiles" ) > 0;
59  bool brute_force = opts.numOptSet( "kdtreeOption" ) > 0;
60 
61  int rank = 0, size = 1;
62 #ifdef MOAB_HAVE_MPI
63  MPI_Init( &argc, &argv );
64  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
65  MPI_Comm_size( MPI_COMM_WORLD, &size );
66 
67  std::string optsRead = ( size == 1 ? ""
68  : std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
69  std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) );
70 #else
71  std::string optsRead;
72 #endif
73  // read meshes in 2 file sets
74  ErrorCode rval;
75  Core* mb = new Core();
76  EntityHandle sf1, sf2, outputSet;
77 
78  // create meshsets and load files
79 
80  rval = mb->create_meshset( MESHSET_SET, sf1 );MB_CHK_ERR( rval );
81  rval = mb->create_meshset( MESHSET_SET, sf2 );MB_CHK_ERR( rval );
82  if( 0 == rank ) std::cout << "Loading mesh file " << sourceFile << "\n";
83  rval = mb->load_file( sourceFile.c_str(), &sf1, optsRead.c_str() );MB_CHK_ERR( rval );
84  if( 0 == rank ) std::cout << "Loading mesh file " << targetFile << "\n";
85  rval = mb->load_file( targetFile.c_str(), &sf2, optsRead.c_str() );MB_CHK_ERR( rval );
86 
87  if( 0 == rank )
88  {
89  std::cout << "Radius: " << R << "\n";
90  std::cout << "relative eps: " << epsrel << "\n";
91  std::cout << "box eps: " << boxeps << "\n";
92  if( brute_force )
93  std::cout << " use kd tree for intersection \n";
94  else
95  std::cout << " use advancing front for intersection \n";
96 
97  std::cout << " area tolerance:" << areaTolerance << "\n";
98  std::cout << " target edge file: " << mapEdgeTargetFile << "\n";
99  }
100  rval = mb->create_meshset( MESHSET_SET, outputSet );MB_CHK_ERR( rval );
101 
102  // fix radius of both meshes, to be consistent with input R
103  rval = moab::IntxUtils::ScaleToRadius( mb, sf1, R );MB_CHK_ERR( rval );
104  rval = moab::IntxUtils::ScaleToRadius( mb, sf2, R );MB_CHK_ERR( rval );
105 
106 #ifdef MOAB_HAVE_MPI
107  ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
108 #endif
109  Intx2MeshEdges worker( mb );
110  IntxAreaUtils areaAdaptor;
111 
112  worker.set_error_tolerance( R * epsrel );
113  worker.set_box_error( boxeps );
114 #ifdef MOAB_HAVE_MPI
115  worker.set_parallel_comm( pcomm );
116 #endif
117  // worker.SetEntityType(moab::MBQUAD);
118  worker.set_radius_source_mesh( R );
119  worker.set_radius_destination_mesh( R );
120  // worker.enable_debug();
121 
122  rval = worker.FindMaxEdges( sf1, sf2 );MB_CHK_ERR( rval );
123 
124  EntityHandle covering_set;
125 #ifdef MOAB_HAVE_MPI
126  if( size > 1 )
127  {
128  Range local_verts;
129  rval = worker.build_processor_euler_boxes( sf2, local_verts );MB_CHK_ERR( rval ); // output also the local_verts
130  if( write_files_rank )
131  {
132  std::stringstream outf;
133  outf << "second_mesh" << rank << ".h5m";
134  rval = mb->write_file( outf.str().c_str(), 0, 0, &sf2, 1 );MB_CHK_ERR( rval );
135  }
136  }
137  if( size > 1 )
138  {
139  double elapsed = MPI_Wtime();
140  rval = mb->create_meshset( moab::MESHSET_SET, covering_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
141  bool gnomonic = true;
142  int order = 0; // we should not need ghost layers here
143  bool include_edges = true; // this is by default false; make it true for this case, for edge maps computation
144  rval = worker.construct_covering_set( sf1, covering_set, gnomonic, order, include_edges );MB_CHK_ERR( rval ); // lots of communication if mesh is distributed very differently
145  elapsed = MPI_Wtime() - elapsed;
146  if( 0 == rank ) std::cout << "\nTime to communicate the mesh = " << elapsed << std::endl;
147  if( write_files_rank )
148  {
149  std::stringstream cof;
150  cof << "covering_mesh" << rank << ".h5m";
151  rval = mb->write_file( cof.str().c_str(), 0, 0, &covering_set, 1 );MB_CHK_ERR( rval );
152  }
153  }
154  else
155 #endif
156  covering_set = sf1;
157 
158  if( 0 == rank ) std::cout << "Computing intersections ..\n";
159 #ifdef MOAB_HAVE_MPI
160  double elapsed = MPI_Wtime();
161 #endif
162  if( brute_force )
163  {
164  rval = worker.intersect_meshes_kdtree( covering_set, sf2, outputSet );MB_CHK_SET_ERR( rval, "failed to intersect meshes with slow method" );
165  }
166  else
167  {
168  rval = worker.intersect_meshes( covering_set, sf2, outputSet );MB_CHK_SET_ERR( rval, "failed to intersect meshes" );
169  }
170 #ifdef MOAB_HAVE_MPI
171  elapsed = MPI_Wtime() - elapsed;
172  if( 0 == rank ) std::cout << "\nTime to compute the intersection between meshes = " << elapsed << std::endl;
173 #endif
174  if( write_files_rank )
175  {
176  std::stringstream outf;
177  outf << "intersect" << rank << ".h5m";
178  rval = mb->write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
179  }
180  double intx_area = areaAdaptor.area_on_sphere( mb, outputSet, R );
181  double arrival_area = areaAdaptor.area_on_sphere( mb, sf2, R );
182  std::cout << "On rank : " << rank << " arrival area: " << arrival_area << " intersection area:" << intx_area
183  << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
184 
185 #ifdef MOAB_HAVE_MPI
186 #ifdef MOAB_HAVE_HDF5_PARALLEL
187  std::ostringstream intx_str;
188  intx_str << "p" << pcomm->size() << "_" << intersectionFile;
189  rval = mb->write_file( intx_str.str().c_str(), 0, "PARALLEL=WRITE_PART", &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
190  if( 0 == rank ) std::cout << " Wrote intx file: " << intx_str.str() << "\n";
191 #else
192  // write intx set on rank 0, in serial; we cannot write in parallel
193  if( 0 == rank )
194  {
195  rval = mb->write_file( intersectionFile.c_str(), 0, 0, &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
196  }
197 #endif
198 
199 #else
200  rval = mb->write_file( intersectionFile.c_str(), 0, 0, &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
201 #endif
202  moab::Tag fractionTag;
203  moab::Tag numSubTag;
204  moab::Tag areaDiffTag;
205  moab::Tag areaTag;
206  double defVal = 0.;
207  rval = mb->tag_get_handle( "EdgeRecoveryFraction", 1, MB_TYPE_DOUBLE, fractionTag, MB_TAG_DENSE | MB_TAG_CREAT,
208  &defVal );MB_CHK_SET_ERR( rval, "can't create edge recovery fraction tag" );
209  rval = mb->tag_get_handle( "NumSubEnts", 1, MB_TYPE_DOUBLE, numSubTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );MB_CHK_SET_ERR( rval, "can't create NumSubEnts tag" );
210  rval = mb->tag_get_handle( "AreaDiff", 1, MB_TYPE_DOUBLE, areaDiffTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );MB_CHK_SET_ERR( rval, "can't create AreaDiff tag" );
211  rval = mb->tag_get_handle( "Area", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );MB_CHK_SET_ERR( rval, "can't create Area tag" );
212  /* std::map<EntityHandle, std::vector<EntityHandle>> edgeVertices; // for each recovered edge, the chain of vertices that form subedges
213  std::map<EntityHandle, std::vector<int>> edgePolygons; // for each recovered edge, the list of intersected polygons;
214  moab::Range recoveredPolys;*/
215 
216  //bool sourceEdgeMap = false;
217  rval = worker.EdgeSplits( areaTolerance );MB_CHK_SET_ERR( rval, "failed to compute edge splits for target" );
218 #ifdef MOAB_HAVE_MPI
219 #ifdef MOAB_HAVE_HDF5_PARALLEL
220  std::ostringstream h5mFile;
221  h5mFile << "p" << pcomm->size() << "_targetWithEdges.h5m";
222  rval = mb->write_file( h5mFile.str().c_str(), 0, "PARALLEL=WRITE_PART", &sf2, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
223  if( 0 == rank ) std::cout << " Wrote file with edge mapping info: " << h5mFile.str() << "\n";
224 #endif
225 #endif
226 
227 #ifdef MOAB_HAVE_PNETCDF
228 #ifdef MOAB_HAVE_MPI
229  std::ostringstream file_str;
230  file_str << "p" << pcomm->size() << "_" << mapEdgeTargetFile;
231  rval = worker.write_edge_map_parallel( file_str.str().c_str() );MB_CHK_SET_ERR( rval, "failed to write edge map for target" );
232  if( 0 == rank ) std::cout << " Wrote netcdf file with edge mapping info: " << file_str.str() << "\n";
233 #else
234 #ifdef MOAB_HAVE_NETCDF
235  rval = worker.write_edge_map( mapEdgeTargetFile.c_str() );MB_CHK_SET_ERR( rval, "failed to write edge map for target" );
236 #endif
237 #endif
238 #endif
239 
240  delete mb;
241 #ifdef MOAB_HAVE_MPI
242  MPI_Finalize();
243 #endif
244  return 0;
245 }