1 /*
2 * Usage: MOAB-Tempest tool
3 *
4 * Generate a Cubed-Sphere mesh: ./mbtempest -t 0 -res 25 -f cubed_sphere_mesh.exo
5 * Generate a RLL mesh: ./mbtempest -t 1 -res 25 -f rll_mesh.exo
6 * Generate a Icosahedral-Sphere mesh: ./mbtempest -t 2 -res 25 <-dual> -f icosa_mesh.exo
7 *
8 * Now you can compute the intersections between the meshes too!
9 *
10 * Generate the overlap mesh: ./mbtempest -t 3 -l cubed_sphere_mesh.exo -l rll_mesh.exo -f overlap_mesh.exo
11 *
12 */
13
14 #include <iostream>
15 #include <cstdlib>
16 #include <vector>
17 #include <string>
18 #include <sstream>
19 #include <cassert>
20
21 #include "moab/Core.hpp"
22 #include "moab/IntxMesh/IntxUtils.hpp"
23 #include "moab/Remapping/TempestRemapper.hpp"
24 #include "moab/Remapping/TempestOfflineMap.hpp"
25 #include "moab/ProgOptions.hpp"
26 #include "moab/CpuTimer.hpp"
27 #include "DebugOutput.hpp"
28
29 //#ifndef MOAB_HAVE_MPI
30 // #error mbtempest tool requires MPI configuration
31 //#endif
32
33 #ifdef MOAB_HAVE_MPI
34 // MPI includes
35 #include "moab_mpi.h"
36 #include "moab/ParallelComm.hpp"
37 #include "MBParallelConventions.h"
38 #endif
39
40 struct ToolContext
41 {
42 moab::Interface* mbcore;
43 #ifdef MOAB_HAVE_MPI
44 moab::ParallelComm* pcomm;
45 #endif
46 int blockSize;
47 std::vector<std::string> inFilenames;
48 std::vector<Mesh*> meshes;
49 std::vector<moab::EntityHandle> meshsets;
50 std::vector<int> disc_orders;
51 std::vector<std::string> disc_methods;
52 std::vector<std::string> doftag_names;
53 std::string outFilename;
54 moab::TempestRemapper::TempestMeshType meshType;
55 const int proc_id, n_procs;
56 bool computeDual;
57 bool computeWeights;
58 int ensureMonotonicity;
59 bool fNoConservation;
60 bool fVolumetric;
61
62 #ifdef MOAB_HAVE_MPI
ToolContextToolContext63 ToolContext ( moab::Interface* icore, moab::ParallelComm* p_pcomm ) :
64 mbcore(icore), pcomm(p_pcomm),
65 proc_id ( pcomm->rank() ), n_procs ( pcomm->size() ),
66 #else
67 ToolContext ( moab::Interface* icore ) :
68 mbcore(icore),
69 proc_id ( 0 ), n_procs ( 1 ),
70 #endif
71 blockSize ( 5 ), outFilename ( "output.exo" ), meshType ( moab::TempestRemapper::DEFAULT ),
72 computeDual ( false ), computeWeights ( false ), ensureMonotonicity ( 0 ),
73 fNoConservation ( false ), fVolumetric ( false )
74 {
75 inFilenames.resize ( 2 );
76 doftag_names.resize( 2 );
77 timer = new moab::CpuTimer();
78 }
79
~ToolContextToolContext80 ~ToolContext()
81 {
82 inFilenames.clear();
83 outFilename.clear();
84 delete timer;
85 }
86
clearToolContext87 void clear()
88 {
89 meshes.clear();
90 }
91
timer_pushToolContext92 void timer_push ( std::string operation )
93 {
94 timer_ops = timer->time_since_birth();
95 opName = operation;
96 }
97
timer_popToolContext98 void timer_pop()
99 {
100 double locElapsed=timer->time_since_birth() - timer_ops, avgElapsed=0, maxElapsed=0;
101 #ifdef MOAB_HAVE_MPI
102 MPI_Reduce(&locElapsed, &maxElapsed, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm());
103 MPI_Reduce(&locElapsed, &avgElapsed, 1, MPI_DOUBLE, MPI_SUM, 0, pcomm->comm());
104 #else
105 maxElapsed = locElapsed;
106 avgElapsed = locElapsed;
107 #endif
108 if (!proc_id) {
109 avgElapsed /= n_procs;
110 std::cout << "[LOG] Time taken to " << opName.c_str() << ": max = " << maxElapsed << ", avg = " << avgElapsed << "\n";
111 }
112 // std::cout << "\n[LOG" << proc_id << "] Time taken to " << opName << " = " << timer->time_since_birth() - timer_ops << std::endl;
113 opName.clear();
114 }
115
ParseCLOptionsToolContext116 void ParseCLOptions ( int argc, char* argv[] )
117 {
118 ProgOptions opts;
119 int imeshType = 0;
120 std::string expectedFName = "output.exo";
121 std::string expectedMethod = "fv";
122 std::string expectedDofTagName = "GLOBAL_ID";
123 int expectedOrder = 1;
124
125 opts.addOpt<int> ( "res,r", "Resolution of the mesh (default=5)", &blockSize );
126 opts.addOpt<int> ( "type,t", "Type of mesh (default=CS; Choose from [CS=0, RLL=1, ICO=2, OVERLAP_FILES=3, OVERLAP_MEMORY=4, OVERLAP_MOAB=5])", &imeshType );
127 opts.addOpt<std::string> ( "file,f", "Output mesh filename (default=output.exo)", &outFilename );
128 opts.addOpt<void> ( "dual,d", "Output the dual of the mesh (generally relevant only for ICO mesh)", &computeDual );
129 opts.addOpt<void> ( "weights,w", "Compute and output the weights using the overlap mesh (generally relevant only for OVERLAP mesh)", &computeWeights );
130 opts.addOpt<void> ( "noconserve,c", "Do not apply conservation to the resultant weights (relevant only when computing weights)", &fNoConservation );
131 opts.addOpt<void> ( "volumetric,v", "Apply a volumetric projection to compute the weights (relevant only when computing weights)", &fVolumetric );
132 opts.addOpt<int> ( "monotonic,n", "Ensure monotonicity in the weight generation", &ensureMonotonicity );
133 opts.addOpt<std::string> ( "load,l", "Input mesh filenames (a source and target mesh)", &expectedFName );
134 opts.addOpt<int> ( "order,o", "Discretization orders for the source and target solution fields", &expectedOrder );
135 opts.addOpt<std::string> ( "method,m", "Discretization method for the source and target solution fields", &expectedMethod );
136 opts.addOpt<std::string> ( "global_id,i", "Tag name that contains the global DoF IDs for source and target solution fields", &expectedDofTagName );
137
138 opts.parseCommandLine ( argc, argv );
139
140 switch ( imeshType )
141 {
142 case 0:
143 meshType = moab::TempestRemapper::CS;
144 break;
145
146 case 1:
147 meshType = moab::TempestRemapper::RLL;
148 break;
149
150 case 2:
151 meshType = moab::TempestRemapper::ICO;
152 break;
153
154 case 3:
155 meshType = moab::TempestRemapper::OVERLAP_FILES;
156 break;
157
158 case 4:
159 meshType = moab::TempestRemapper::OVERLAP_MEMORY;
160 break;
161
162 case 5:
163 meshType = moab::TempestRemapper::OVERLAP_MOAB;
164 break;
165
166 default:
167 meshType = moab::TempestRemapper::DEFAULT;
168 break;
169 }
170
171 if ( meshType > moab::TempestRemapper::ICO )
172 {
173 opts.getOptAllArgs ( "load,l", inFilenames );
174 opts.getOptAllArgs ( "order,o", disc_orders );
175 opts.getOptAllArgs ( "method,m", disc_methods );
176 opts.getOptAllArgs ( "global_id,i", doftag_names );
177
178 if ( disc_orders.size() == 0 )
179 { disc_orders.resize ( 2, 1 ); }
180
181 if ( disc_orders.size() == 1 )
182 { disc_orders.push_back ( 1 ); }
183
184 if ( disc_methods.size() == 0 )
185 { disc_methods.resize ( 2, "fv" ); }
186
187 if ( disc_methods.size() == 1 )
188 { disc_methods.push_back ( "fv" ); }
189
190 if ( doftag_names.size() == 0 )
191 { doftag_names.resize ( 2, "GLOBAL_ID" ); }
192
193 if ( doftag_names.size() == 1 )
194 { doftag_names.push_back ( "GLOBAL_ID" ); }
195
196 assert ( inFilenames.size() == 2 );
197 assert ( disc_orders.size() == 2 );
198 assert ( disc_methods.size() == 2 );
199 }
200
201 expectedFName.clear();
202 }
203 private:
204 moab::CpuTimer* timer;
205 double timer_ops;
206 std::string opName;
207 };
208
209 // Forward declare some methods
210 moab::ErrorCode CreateTempestMesh ( ToolContext&, moab::TempestRemapper& remapper, Mesh* );
211
main(int argc,char * argv[])212 int main ( int argc, char* argv[] )
213 {
214 moab::ErrorCode rval;
215 NcError error ( NcError::verbose_nonfatal );
216 std::stringstream sstr;
217
218 int proc_id = 0, nprocs = 1;
219 #ifdef MOAB_HAVE_MPI
220 MPI_Init ( &argc, &argv );
221 MPI_Comm_rank ( MPI_COMM_WORLD, &proc_id );
222 MPI_Comm_size ( MPI_COMM_WORLD, &nprocs );
223 #endif
224
225 moab::Interface* mbCore = new ( std::nothrow ) moab::Core;
226
227 if ( NULL == mbCore ) { return 1; }
228
229 #ifdef MOAB_HAVE_MPI
230 moab::ParallelComm* pcomm = new moab::ParallelComm ( mbCore, MPI_COMM_WORLD, 0 );
231
232 ToolContext ctx ( mbCore, pcomm );
233 #else
234 ToolContext ctx ( mbCore );
235 #endif
236 ctx.ParseCLOptions ( argc, argv );
237
238 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
239 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
240
241 #ifdef MOAB_HAVE_MPI
242 moab::TempestRemapper remapper ( mbCore, pcomm );
243 #else
244 moab::TempestRemapper remapper ( mbCore );
245 #endif
246 remapper.meshValidate = true;
247 remapper.constructEdgeMap = true;
248 remapper.initialize();
249
250 Mesh* tempest_mesh = new Mesh();
251 ctx.timer_push ( "create Tempest mesh" );
252 rval = CreateTempestMesh ( ctx, remapper, tempest_mesh ); MB_CHK_ERR ( rval );
253 ctx.timer_pop();
254
255 // Some constant parameters
256 const double epsrel = 1.e-15;
257 const double boxeps = 1e-6;
258
259 if ( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
260 {
261 // Compute intersections with MOAB
262 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
263 assert ( ctx.meshes.size() == 3 );
264
265 #ifdef MOAB_HAVE_MPI
266 rval = pcomm->check_all_shared_handles(); MB_CHK_ERR ( rval );
267 #endif
268
269 // Load the meshes and validate
270 rval = remapper.ConvertTempestMesh ( moab::Remapper::SourceMesh ); MB_CHK_ERR ( rval );
271 rval = remapper.ConvertTempestMesh ( moab::Remapper::TargetMesh ); MB_CHK_ERR ( rval );
272 remapper.SetMeshType ( moab::Remapper::IntersectedMesh, moab::TempestRemapper::OVERLAP_FILES );
273 rval = remapper.ConvertTempestMesh ( moab::Remapper::IntersectedMesh ); MB_CHK_ERR ( rval );
274 rval = mbCore->write_mesh ( "tempest_intersection.h5m", &ctx.meshsets[2], 1 ); MB_CHK_ERR ( rval );
275
276 // print verbosely about the problem setting
277 {
278 moab::Range rintxverts, rintxelems;
279 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[0], 0, rintxverts ); MB_CHK_ERR ( rval );
280 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[0], 2, rintxelems ); MB_CHK_ERR ( rval );
281 printf ( "The red set contains %lu vertices and %lu elements \n", rintxverts.size(), rintxelems.size() );
282
283 moab::Range bintxverts, bintxelems;
284 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[1], 0, bintxverts ); MB_CHK_ERR ( rval );
285 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[1], 2, bintxelems ); MB_CHK_ERR ( rval );
286 printf ( "The blue set contains %lu vertices and %lu elements \n", bintxverts.size(), bintxelems.size() );
287 }
288
289 moab::EntityHandle intxset; // == remapper.GetMeshSet(moab::Remapper::IntersectedMesh);
290
291 // Compute intersections with MOAB
292 {
293 // Create the intersection on the sphere object
294 ctx.timer_push ( "setup the intersector" );
295
296 moab::Intx2MeshOnSphere* mbintx = new moab::Intx2MeshOnSphere ( mbCore );
297 mbintx->set_error_tolerance ( epsrel );
298 mbintx->set_box_error ( boxeps );
299 mbintx->set_radius_source_mesh ( radius_src );
300 mbintx->set_radius_destination_mesh ( radius_dest );
301 #ifdef MOAB_HAVE_MPI
302 mbintx->set_parallel_comm ( pcomm );
303 #endif
304 rval = mbintx->FindMaxEdges ( ctx.meshsets[0], ctx.meshsets[1] ); MB_CHK_ERR ( rval );
305
306 #ifdef MOAB_HAVE_MPI
307 moab::Range local_verts;
308 rval = mbintx->build_processor_euler_boxes ( ctx.meshsets[1], local_verts ); MB_CHK_ERR ( rval );
309
310 ctx.timer_pop();
311
312 moab::EntityHandle covering_set;
313 ctx.timer_push ( "communicate the mesh" );
314 rval = mbintx->construct_covering_set ( ctx.meshsets[0], covering_set ); MB_CHK_ERR ( rval ); // lots of communication if mesh is distributed very differently
315 ctx.timer_pop();
316 #else
317 moab::EntityHandle covering_set = ctx.meshsets[0];
318 #endif
319 // Now let's invoke the MOAB intersection algorithm in parallel with a
320 // source and target mesh set representing two different decompositions
321 ctx.timer_push ( "compute intersections with MOAB" );
322 rval = mbCore->create_meshset ( moab::MESHSET_SET, intxset ); MB_CHK_SET_ERR ( rval, "Can't create new set" );
323 rval = mbintx->intersect_meshes ( covering_set, ctx.meshsets[1], intxset ); MB_CHK_SET_ERR ( rval, "Can't compute the intersection of meshes on the sphere" );
324 ctx.timer_pop();
325
326 // free the memory
327 delete mbintx;
328 }
329
330 {
331 moab::Range intxelems, intxverts;
332 rval = mbCore->get_entities_by_dimension ( intxset, 2, intxelems ); MB_CHK_ERR ( rval );
333 rval = mbCore->get_entities_by_dimension ( intxset, 0, intxverts, true ); MB_CHK_ERR ( rval );
334 printf ( "The intersection set contains %lu elements and %lu vertices \n", intxelems.size(), intxverts.size() );
335
336 double initial_area = area_on_sphere_lHuiller ( mbCore, ctx.meshsets[0], radius_src ); // use the target to compute the initial area
337 double area_method1 = area_on_sphere_lHuiller ( mbCore, intxset, radius_src );
338 double area_method2 = area_on_sphere ( mbCore, intxset, radius_src );
339
340 printf ( "initial area: %12.10f\n", initial_area );
341 printf ( " area with l'Huiller: %12.10f with Girard: %12.10f\n", area_method1, area_method2 );
342 printf ( " relative difference areas = %12.10e\n", fabs ( area_method1 - area_method2 ) / area_method1 );
343 printf ( " relative error = %12.10e\n", fabs ( area_method1 - initial_area ) / area_method1 );
344 }
345
346 // Write out our computed intersection file
347 rval = mbCore->write_mesh ( "moab_intersection.h5m", &intxset, 1 ); MB_CHK_ERR ( rval );
348
349 if ( ctx.computeWeights )
350 {
351 ctx.timer_push ( "compute weights with the Tempest meshes" );
352 // Call to generate an offline map with the tempest meshes
353 OfflineMap weightMap;
354 int err = GenerateOfflineMapWithMeshes ( weightMap, *ctx.meshes[0], *ctx.meshes[1], *ctx.meshes[2],
355 "", "", // std::string strInputMeta, std::string strOutputMeta,
356 ctx.disc_methods[0], ctx.disc_methods[1], // std::string strInputType, std::string strOutputType,
357 ctx.disc_orders[0], ctx.disc_orders[1] // int nPin=4, int nPout=4,
358 );
359 ctx.timer_pop();
360
361 if ( err ) { rval = moab::MB_FAILURE; }
362 else { weightMap.Write ( "outWeights.nc" ); }
363 }
364 }
365 else if ( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
366 {
367 // Usage: mpiexec -n 2 tools/mbtempest -t 5 -l mycs_2.h5m -l myico_2.h5m -f myoverlap_2.h5m
368 #ifdef MOAB_HAVE_MPI
369 rval = pcomm->check_all_shared_handles(); MB_CHK_ERR ( rval );
370 #endif
371 // print verbosely about the problem setting
372 {
373 moab::Range rintxverts, rintxelems;
374 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[0], 0, rintxverts ); MB_CHK_ERR ( rval );
375 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[0], 2, rintxelems ); MB_CHK_ERR ( rval );
376 rval = fix_degenerate_quads ( mbCore, ctx.meshsets[0] ); MB_CHK_ERR ( rval );
377 rval = positive_orientation ( mbCore, ctx.meshsets[0], radius_src ); MB_CHK_ERR ( rval );
378 printf ( "The red set contains %lu vertices and %lu elements \n", rintxverts.size(), rintxelems.size() );
379
380 moab::Range bintxverts, bintxelems;
381 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[1], 0, bintxverts ); MB_CHK_ERR ( rval );
382 rval = mbCore->get_entities_by_dimension ( ctx.meshsets[1], 2, bintxelems ); MB_CHK_ERR ( rval );
383 rval = fix_degenerate_quads ( mbCore, ctx.meshsets[1] ); MB_CHK_ERR ( rval );
384 rval = positive_orientation ( mbCore, ctx.meshsets[1], radius_dest ); MB_CHK_ERR ( rval );
385 printf ( "The blue set contains %lu vertices and %lu elements \n", bintxverts.size(), bintxelems.size() );
386 }
387
388 // Compute intersections with MOAB
389 ctx.timer_push ( "setup and compute mesh intersections" );
390 rval = remapper.ComputeOverlapMesh ( epsrel ); MB_CHK_ERR ( rval );
391 ctx.timer_pop();
392
393 // Write out our computed intersection file
394 if ( ctx.n_procs == 1 )
395 {
396 printf ( "writing out the intersection mesh file to %s\n", "moab_intersection.h5m" );
397 // rval = mbCore->add_entities ( ctx.meshsets[2], &ctx.meshsets[0], 2 ); MB_CHK_ERR ( rval );
398 rval = mbCore->write_file ( "moab_intersection.h5m", NULL, "PARALLEL=WRITE_PART", &ctx.meshsets[0], 3 ); MB_CHK_ERR ( rval );
399 }
400
401 {
402 double local_areas[3], global_areas[3]; // Array for Initial area, and through Method 1 and Method 2
403 // local_areas[0] = area_on_sphere_lHuiller ( mbCore, ctx.meshsets[1], radius_src );
404 local_areas[0] = area_on_sphere ( mbCore, ctx.meshsets[1], radius_src );
405 local_areas[1] = area_on_sphere_lHuiller ( mbCore, ctx.meshsets[2], radius_src );
406 local_areas[2] = area_on_sphere ( mbCore, ctx.meshsets[2], radius_src );
407
408 #ifdef MOAB_HAVE_MPI
409 MPI_Allreduce ( &local_areas[0], &global_areas[0], 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
410 #else
411 global_areas[0] = local_areas[0];
412 global_areas[1] = local_areas[1];
413 global_areas[2] = local_areas[2];
414 #endif
415 if ( !proc_id )
416 {
417 printf ( "initial area: %12.10f\n", global_areas[0] );
418 printf ( " area with l'Huiller: %12.10f with Girard: %12.10f\n", global_areas[1], global_areas[2] );
419 printf ( " relative difference areas = %12.10e\n", fabs ( global_areas[1] - global_areas[2] ) / global_areas[1] );
420 printf ( " relative error = %12.10e\n", fabs ( global_areas[1] - global_areas[0] ) / global_areas[1] );
421 }
422 }
423
424 if ( ctx.computeWeights )
425 {
426 ctx.meshes[2] = remapper.GetMesh ( moab::Remapper::IntersectedMesh );
427
428 ctx.timer_push ( "setup computation of weights" );
429 // Call to generate an offline map with the tempest meshes
430 moab::TempestOfflineMap* weightMap = new moab::TempestOfflineMap ( &remapper );
431 ctx.timer_pop();
432
433 ctx.timer_push ( "compute weights with TempestRemap" );
434
435 rval = weightMap->GenerateOfflineMap ( ctx.disc_methods[0], ctx.disc_methods[1], // std::string strInputType, std::string strOutputType,
436 ctx.disc_orders[0], ctx.disc_orders[1], // int nPin=4, int nPout=4,
437 false, ctx.ensureMonotonicity, // bool fBubble=false, int fMonotoneTypeID=0,
438 ctx.fVolumetric, ctx.fNoConservation, false, // bool fVolumetric=false, bool fNoConservation=false, bool fNoCheck=false,
439 ctx.doftag_names[0], ctx.doftag_names[1],
440 "", //"", // std::string strVariables="", std::string strOutputMap="",
441 "", "", // std::string strInputData="", std::string strOutputData="",
442 "", false, // std::string strNColName="", bool fOutputDouble=false,
443 "", false, 0.0, // std::string strPreserveVariables="", bool fPreserveAll=false, double dFillValueOverride=0.0,
444 false, false // bool fInputConcave = false, bool fOutputConcave = false
445 );MB_CHK_ERR ( rval );
446 ctx.timer_pop();
447
448 #if 0
449 // OfflineMap* weightMap;
450 // weightMap = GenerateOfflineMapWithMeshes( NULL, *ctx.meshes[0], *ctx.meshes[1], *ctx.meshes[2],
451 // "", "", // std::string strInputMeta, std::string strOutputMeta,
452 // ctx.disc_methods[0], ctx.disc_methods[1],// std::string strInputType, std::string strOutputType,
453 // ctx.disc_orders[0], ctx.disc_orders[1], // int nPin=4, int nPout=4,
454 // false, 0, // bool fBubble=false, int fMonotoneTypeID=0,
455 // false, false, (nprocs>1) // bool fVolumetric=false, bool fNoConservation=false, bool fNoCheck=false,
456 // // std::string strVariables="", std::string strOutputMap="",
457 // // std::string strInputData="", std::string strOutputData="",
458 // // std::string strNColName="", bool fOutputDouble=false,
459 // // std::string strPreserveVariables="", bool fPreserveAll=false, double dFillValueOverride=0.0
460 // );
461
462 // weightMap->m_vecSourceDimSizes.resize(ctx.meshes[0]->faces.size());
463 // weightMap->m_vecTargetDimSizes.resize(ctx.meshes[1]->faces.size());
464
465 #endif
466
467 /*
468 ctx.timer_push ( "apply weights onto a vector" );
469 rval = weightMap->ApplyWeights ( srcVals, tgtvals, false);MB_CHK_ERR ( rval );
470 ctx.timer_pop();
471 */
472
473 /*
474 * the file can be written in parallel, and it will contain additional tags defined by the user
475 * we may extend the method to write only desired tags to the file
476 */
477 if (nprocs == 1) {
478 // free allocated data
479 char outputFileTgt[] = "fIntxTarget.h5m";
480 #ifdef MOAB_HAVE_MPI
481 const char *writeOptions = (nprocs > 1 ? "PARALLEL=WRITE_PART" : "");
482 #else
483 const char *writeOptions = "";
484 #endif
485
486 rval = mbCore->write_file ( outputFileTgt, NULL, writeOptions, &ctx.meshsets[2], 1 ); MB_CHK_ERR ( rval );
487
488 }
489
490 #ifdef VERBOSE
491 sstr.str("");
492 sstr << "outWeights_" << proc_id << ".nc";
493 weightMap->Write(sstr.str().c_str());
494 sstr.str("");
495 #endif
496
497 // sstr << "newoutWeights_" << proc_id << ".nc";
498 // weightMap->WriteParallelWeightsToFile(sstr.str().c_str());
499
500 delete weightMap;
501 }
502 }
503
504 // Clean up
505 ctx.clear();
506 delete mbCore;
507
508 #ifdef MOAB_HAVE_MPI
509 MPI_Finalize();
510 #endif
511 exit ( 0 );
512 }
513
514
CreateTempestMesh(ToolContext & ctx,moab::TempestRemapper & remapper,Mesh * tempest_mesh)515 moab::ErrorCode CreateTempestMesh ( ToolContext& ctx, moab::TempestRemapper& remapper, Mesh* tempest_mesh )
516 {
517 moab::ErrorCode rval = moab::MB_SUCCESS;
518 int err;
519
520 if ( !ctx.proc_id ) { printf ( "Creating TempestRemap Mesh object ...\n" ); }
521
522 if ( ctx.meshType == moab::TempestRemapper::OVERLAP_FILES )
523 {
524 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
525 err = GenerateOverlapMesh ( ctx.inFilenames[0], ctx.inFilenames[1], *tempest_mesh, ctx.outFilename, "exact", true );
526
527 if ( err ) { rval = moab::MB_FAILURE; }
528 else
529 {
530 ctx.meshes.push_back ( tempest_mesh );
531 }
532 }
533 else if ( ctx.meshType == moab::TempestRemapper::OVERLAP_MEMORY )
534 {
535 // Load the meshes and validate
536 ctx.meshsets.resize ( 3 );
537 ctx.meshes.resize ( 3 );
538 ctx.meshsets[0] = remapper.GetMeshSet ( moab::Remapper::SourceMesh );
539 ctx.meshsets[1] = remapper.GetMeshSet ( moab::Remapper::TargetMesh );
540 ctx.meshsets[2] = remapper.GetMeshSet ( moab::Remapper::IntersectedMesh );
541
542 // First the source
543 rval = remapper.LoadMesh ( moab::Remapper::SourceMesh, ctx.inFilenames[0], moab::TempestRemapper::DEFAULT ); MB_CHK_ERR ( rval );
544 ctx.meshes[0] = remapper.GetMesh ( moab::Remapper::SourceMesh );
545
546 // Next the target
547 rval = remapper.LoadMesh ( moab::Remapper::TargetMesh, ctx.inFilenames[1], moab::TempestRemapper::DEFAULT ); MB_CHK_ERR ( rval );
548 ctx.meshes[1] = remapper.GetMesh ( moab::Remapper::TargetMesh );
549
550 // Now let us construct the overlap mesh, by calling TempestRemap interface directly
551 // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
552 err = GenerateOverlapWithMeshes ( *ctx.meshes[0], *ctx.meshes[1], *tempest_mesh, "" /*ctx.outFilename*/, "exact", false );
553
554 if ( err ) { rval = moab::MB_FAILURE; }
555 else
556 {
557 remapper.SetMesh ( moab::Remapper::IntersectedMesh, tempest_mesh );
558 ctx.meshes[2] = remapper.GetMesh ( moab::Remapper::IntersectedMesh );
559 // ctx.meshes.push_back(*tempest_mesh);
560 }
561 }
562 else if ( ctx.meshType == moab::TempestRemapper::OVERLAP_MOAB )
563 {
564 ctx.meshsets.resize ( 3 );
565 ctx.meshes.resize ( 3 );
566 ctx.meshsets[0] = remapper.GetMeshSet ( moab::Remapper::SourceMesh );
567 ctx.meshsets[1] = remapper.GetMeshSet ( moab::Remapper::TargetMesh );
568 ctx.meshsets[2] = remapper.GetMeshSet ( moab::Remapper::IntersectedMesh );
569
570 const double radius_src = 1.0 /*2.0*acos(-1.0)*/;
571 const double radius_dest = 1.0 /*2.0*acos(-1.0)*/;
572
573 // Load the source mesh and validate
574 rval = remapper.LoadNativeMesh ( ctx.inFilenames[0], ctx.meshsets[0], 0 ); MB_CHK_ERR ( rval );
575 // Rescale the radius of both to compute the intersection
576 rval = ScaleToRadius(ctx.mbcore, ctx.meshsets[0], radius_src);MB_CHK_ERR ( rval );
577 rval = remapper.ConvertMeshToTempest ( moab::Remapper::SourceMesh ); MB_CHK_ERR ( rval );
578 ctx.meshes[0] = remapper.GetMesh ( moab::Remapper::SourceMesh );
579
580 // Load the target mesh and validate
581 rval = remapper.LoadNativeMesh ( ctx.inFilenames[1], ctx.meshsets[1], 0 ); MB_CHK_ERR ( rval );
582 rval = ScaleToRadius(ctx.mbcore, ctx.meshsets[1], radius_dest);MB_CHK_ERR ( rval );
583 rval = remapper.ConvertMeshToTempest ( moab::Remapper::TargetMesh ); MB_CHK_ERR ( rval );
584 ctx.meshes[1] = remapper.GetMesh ( moab::Remapper::TargetMesh );
585
586 // Set the references for the overlap mesh
587 // ctx.meshes[2] = remapper.GetMesh(moab::Remapper::IntersectedMesh);
588 }
589 else if ( ctx.meshType == moab::TempestRemapper::ICO )
590 {
591 err = GenerateICOMesh ( *tempest_mesh, ctx.blockSize, ctx.computeDual, ctx.outFilename );
592
593 if ( err ) { rval = moab::MB_FAILURE; }
594 else
595 {
596 ctx.meshes.push_back ( tempest_mesh );
597 }
598 }
599 else if ( ctx.meshType == moab::TempestRemapper::RLL )
600 {
601 err = GenerateRLLMesh ( *tempest_mesh, // Mesh& meshOut,
602 ctx.blockSize * 2, ctx.blockSize, // int nLongitudes, int nLatitudes,
603 0.0, 360.0, // double dLonBegin, double dLonEnd,
604 -90.0, 90.0, // double dLatBegin, double dLatEnd,
605 false, false, // bool fFlipLatLon, bool fForceGlobal,
606 "" /*ctx.inFilename*/, ctx.outFilename, // std::string strInputFile, std::string strOutputFile,
607 true // bool fVerbose
608 );
609
610 if ( err ) { rval = moab::MB_FAILURE; }
611 else
612 {
613 ctx.meshes.push_back ( tempest_mesh );
614 }
615 }
616 else // default
617 {
618 err = GenerateCSMesh ( *tempest_mesh, ctx.blockSize, true, ctx.outFilename );
619
620 if ( err ) { rval = moab::MB_FAILURE; }
621 else
622 {
623 ctx.meshes.push_back ( tempest_mesh );
624 }
625 }
626
627 if ( ctx.meshType != moab::TempestRemapper::OVERLAP_MOAB && !tempest_mesh )
628 {
629 std::cout << "Tempest Mesh is not a complete object; Quitting...";
630 exit ( -1 );
631 }
632
633 return rval;
634 }
635