1 #include <iostream>
2 #include <stdlib.h>
3 #include <vector>
4 #include <string>
5 #include <stdio.h>
6 #include <iomanip>
7 #include <fstream>
8 #include <sys/time.h>
9 #include <time.h>
10 #include <math.h>
11 #include <assert.h>
12 #include <float.h>
13
14 #include "moab/Core.hpp"
15 #include "moab/Interface.hpp"
16 #include "moab/verdict/VerdictWrapper.hpp"
17 #include "moab/NestedRefine.hpp"
18
19 #ifdef MOAB_HAVE_MPI
20 #include "moab_mpi.h"
21 #include "moab/ParallelComm.hpp"
22 #include "MBParallelConventions.h"
23 #endif
24
25 /* Exit values */
26 #define SUCCESS 0
27 #define USAGE_ERROR 1
28 #define NOT_IMPLEMENTED 2
29
30 using namespace moab;
31
print_usage(const char * name,std::ostream & stream)32 static void print_usage( const char* name, std::ostream& stream )
33 {
34 stream << "Usage: " << name
35 << " <options> <input_file> [<input_file2> ...]" << std::endl
36 << "Options: " << std::endl
37 << "\t-h - Print this help text and exit." << std::endl
38 << "\t-d - Dimension of the mesh."<<std::endl
39 << "\t-n - Exact or a maximum number of levels for the hierarchy. Default 1." << std::endl
40 << "\t-L - Degree of refinement for each level. Pass an array or a number. It is mandatory to pass dimension and num_levels before to use this option. If this flag is not used, a default of deg 2 refinement is used. " << std::endl
41 << "\t-V - Pass a desired volume (absolute) . This will generate a hierarchy such that the maximum volume is reduced to the given amount approximately. The length of the hierarchy can be constrained further if a maximum number of levels is passed. It is mandatory to pass the dimension for this option. " << std::endl
42 <<"\t-q - Prints out the maximum volume of the mesh and exits the program. This option can be used as a guide to volume constrainted mesh hierarchy later. "<<std::endl
43 << "\t-t - Print out the time taken to generate hierarchy." << std::endl
44 <<"\t-s - Print out the mesh sizes of each level of the generated hierarchy." << std::endl
45 << "\t-o - Specify true for output files for the mesh levels of the hierarchy." << std::endl
46 //<< "\t-O option - Specify read option." << std::endl
47 #ifdef MOAB_HAVE_MPI
48 << "\t-p[0|1|2] - Read in parallel[0], optionally also doing resolve_shared_ents (1) and exchange_ghosts (2)" << std::endl
49 #endif
50 ;
51 exit(USAGE_ERROR);
52 }
53
usage_error(const char * name)54 static void usage_error( const char* name )
55 {
56 print_usage( name, std::cerr );
57 #ifdef MOAB_HAVE_MPI
58 MPI_Finalize();
59 #endif
60 exit(USAGE_ERROR);
61 }
62
63 bool parse_id_list(const char* string, int dim, int nval, std::vector<int> &results );
64
65 bool make_opts_string( std::vector<std::string> options, std::string& opts );
66
67 ErrorCode get_degree_seq(Core &mb, EntityHandle fileset, int dim, double desired_vol, int &num_levels, std::vector<int> &level_degs);
68
69 ErrorCode get_max_volume(Core &mb, EntityHandle fileset, int dim, double &vmax);
70
main(int argc,char * argv[])71 int main(int argc, char* argv[])
72 {
73 int proc_id = 0, size = 1;
74 #ifdef MOAB_HAVE_MPI
75 MPI_Init(&argc,&argv);
76 MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
77 MPI_Comm_size( MPI_COMM_WORLD, &size );
78 #endif
79
80 int num_levels = 0, dim=0;
81 std::vector<int> level_degrees;
82 bool optimize = false;
83 bool do_flag = true;
84 bool print_times = false, print_size=false, output = false;
85 bool parallel = false, resolve_shared = false, exchange_ghosts = false;
86 bool printusage = false, parselevels = true;
87 bool qc_vol = false, only_quality = false; double cvol = 0 ;
88 std::string infile;
89
90 int i=1;
91 while ( i<argc)
92 {
93 if (!argv[i][0] && 0==proc_id)
94 {
95 usage_error(argv[0]);
96 #ifdef MOAB_HAVE_MPI
97 MPI_Finalize();
98 #endif
99 }
100
101 if (do_flag && argv[i][0] == '-')
102 {
103 switch ( argv[i][1] )
104 {
105 // do flag arguments:
106 case '-': do_flag = false; break;
107 case 't': print_times = true; break;
108 case 's': print_size = true; break;
109 case 'h':
110 case 'H': print_usage( argv[0], std::cerr ); printusage = true; break;
111 case 'd': dim = atol(argv[i+1]); ++i; break;
112 case 'n': num_levels = atol(argv[i+1]); ++i; break;
113 case 'L': if (dim !=0 && num_levels !=0)
114 { parselevels = parse_id_list(argv[i+1], dim, num_levels, level_degrees); ++i;}
115 else
116 {print_usage( argv[0], std::cerr ); printusage = true;}
117 break;
118 case 'V': qc_vol = true; cvol = strtod(argv[i+1], NULL); ++i; break;
119 case 'q': only_quality = true; break;
120 case 'o': output = true; break;
121 case 'O': optimize = true; break;
122 #ifdef MOAB_HAVE_MPI
123 case 'p':
124 parallel = true;
125 resolve_shared = true;
126 if (argv[i][1] == '1') exchange_ghosts = true;
127 break;
128 #endif
129 default:
130 ++i;
131 switch ( argv[i-1][1] )
132 {
133 //case 'O': read_opts.push_back(argv[i]); break;
134 default: std::cerr << "Invalid option: " << argv[i] << std::endl;
135 }
136 }
137 ++i;
138 }
139 // do file names
140 else {
141 infile = argv[i] ;
142 ++i;
143 }
144 }
145
146 if (infile.empty() && !printusage)
147 print_usage(argv[0], std::cerr);
148
149 if (!parselevels)
150 exit(USAGE_ERROR);
151
152 #ifdef MOAB_HAVE_MPI
153 parallel = true;
154 resolve_shared = true;
155 #endif
156
157 ErrorCode error;
158
159 Core *moab = new Core;
160 Interface *mb = (Interface*)moab;
161 EntityHandle fileset;
162
163 //Create a fileset
164 error = mb->create_meshset(MESHSET_SET, fileset);MB_CHK_ERR(error);
165
166 //Set the read options for parallel file loading
167 std::vector<std::string> read_opts, write_opts;
168 std::string read_options, write_options;
169
170 if (parallel && size > 1){
171 read_opts.push_back("PARALLEL=READ_PART");
172 read_opts.push_back("PARTITION=PARALLEL_PARTITION");
173 if (resolve_shared) read_opts.push_back("PARALLEL_RESOLVE_SHARED_ENTS");
174 if (exchange_ghosts) read_opts.push_back("PARALLEL_GHOSTS=3.0.1");
175
176 /* if (output)
177 {
178 write_opts.push_back(";;PARALLEL=WRITE_PART");
179 std::cout<<"Here"<<std::endl;
180 }*/
181 }
182
183 if (!make_opts_string( read_opts, read_options ) || !make_opts_string( write_opts, write_options ))
184 {
185 #ifdef MOAB_HAVE_MPI
186 MPI_Finalize();
187 #endif
188 return USAGE_ERROR;
189 }
190
191 //Load file
192 std::cout<<"READ OPTS="<<(char*)read_options.c_str()<<std::endl;
193 error = mb->load_file(infile.c_str(), &fileset, read_options.c_str());MB_CHK_ERR(error);
194
195 //Create the nestedrefine instance
196
197 #ifdef MOAB_HAVE_MPI
198 ParallelComm *pc = new ParallelComm(mb, MPI_COMM_WORLD);
199 NestedRefine * uref = new NestedRefine(moab, pc, fileset);
200 #else
201 NestedRefine * uref = new NestedRefine(moab);
202 #endif
203
204 std::vector<EntityHandle> lsets;
205
206 // std::cout<<"Read input file"<<std::endl;
207
208 if (only_quality)
209 {
210 double vmax;
211 error = get_max_volume(*moab, fileset, dim, vmax);MB_CHK_ERR(error);
212 #ifdef MOAB_HAVE_MPI
213 int rank = 0;
214 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
215 if (rank==0)
216 std::cout<<"Maximum global volume = "<<vmax<<std::endl;
217 #else
218 std::cout<<"Maximum volume = "<<vmax<<std::endl;
219 #endif
220 exit(SUCCESS);
221 }
222
223 //If a volume constraint is given, find an optimal degree sequence to reach the desired volume constraint.
224 if (qc_vol){
225 error = get_degree_seq(*moab, fileset, dim, cvol, num_levels, level_degrees);MB_CHK_ERR(error);
226
227 if (dim==0)
228 print_usage(argv[0], std::cerr);
229 }
230
231 if (num_levels == 0)
232 num_levels = 1;
233
234 int *ldeg = new int[num_levels];
235 //std::cout<<"level_degrees.size = "<<level_degrees.size()<<std::endl;
236 if (level_degrees.empty()){
237 for (int l=0; l<num_levels; l++)
238 ldeg[l] = 2;
239 }
240 else
241 {
242 for (int l=0; l<num_levels; l++)
243 ldeg[l] = level_degrees[l];
244 }
245
246 std::cout<<"Starting hierarchy generation"<<std::endl;
247
248 std::cout<<"opt = "<<optimize<<std::endl;
249
250 error = uref->generate_mesh_hierarchy( num_levels, ldeg, lsets, optimize);MB_CHK_ERR(error);
251
252 if (print_times)
253 {
254 std::cout<<"Finished hierarchy generation in "<<uref->timeall.tm_total<<" secs"<<std::endl;
255 if (parallel)
256 {
257 std::cout<<"Time taken for refinement "<<uref->timeall.tm_refine<<" secs"<<std::endl;
258 std::cout<<"Time taken for resolving shared interface "<<uref->timeall.tm_resolve<<" secs"<<std::endl;
259 }
260 }
261 else
262 std::cout<<"Finished hierarchy generation "<<std::endl;
263
264 if (print_size)
265 {
266 Range all_ents, ents[4];
267 error = mb->get_entities_by_handle(fileset, all_ents); MB_CHK_ERR(error);
268
269 for (int k=0; k<4; k++)
270 ents[k] = all_ents.subset_by_dimension(k);
271
272 std::cout<<std::endl;
273 if (qc_vol)
274 {
275 double volume;
276 error = get_max_volume(*moab, fileset, dim, volume);MB_CHK_ERR(error);
277 std::cout<<"Mesh size for level 0"<<" :: nverts = "<<ents[0].size()<<", nedges = "<<ents[1].size()<<", nfaces = "<<ents[2].size()<<", ncells = "<<ents[3].size()<<" :: Vmax = "<<volume<<std::endl;
278 }
279 else
280 std::cout<<"Mesh size for level 0"<<" :: nverts = "<<ents[0].size()<<", nedges = "<<ents[1].size()<<", nfaces = "<<ents[2].size()<<", ncells = "<<ents[3].size()<<std::endl;
281
282 for (int l=0; l<num_levels; l++)
283 {
284 all_ents.clear();
285 ents[0].clear(); ents[1].clear(); ents[2].clear(); ents[3].clear();
286 error = mb->get_entities_by_handle(lsets[l+1], all_ents); MB_CHK_ERR(error);
287
288 for (int k=0; k<4; k++)
289 ents[k] = all_ents.subset_by_dimension(k);
290
291 // std::cout<<std::endl;
292
293 if (qc_vol)
294 {
295 double volume;
296 error = get_max_volume(*moab, lsets[l+1], dim, volume);MB_CHK_ERR(error);
297 std::cout<<"Mesh size for level "<<l+1<<" :: nverts = "<<ents[0].size()<<", nedges = "<<ents[1].size()<<", nfaces = "<<ents[2].size()<<", ncells = "<<ents[3].size()<<" :: Vmax = "<<volume<<std::endl;
298 }
299 else
300 std::cout<<"Mesh size for level "<<l+1<<" :: nverts = "<<ents[0].size()<<", nedges = "<<ents[1].size()<<", nfaces = "<<ents[2].size()<<", ncells = "<<ents[3].size()<<std::endl;
301 }
302 }
303
304
305 if (output){
306 for (int l=0; l<num_levels; l++){
307 std::string::size_type idx1 = infile.find_last_of("\\/");
308 std::string::size_type idx2 = infile.find_last_of(".");
309 std::string file = infile.substr(idx1+1,idx2-idx1-1);
310 std::stringstream out;
311 if (parallel)
312 out << "_ML_" <<l+1<<".h5m";
313 else
314 out << "_ML_" <<l+1<<".vtk";
315 file = file + out.str();
316 const char* output_file = file.c_str();
317 #ifdef MOAB_HAVE_MPI
318 error = mb->write_file(output_file, 0, ";;PARALLEL=WRITE_PART",&lsets[l+1], 1);MB_CHK_ERR(error);
319 #else
320 error = mb->write_file(output_file, 0, NULL, &lsets[l+1], 1); MB_CHK_ERR(error);
321 #endif
322 // const char* output_file = file.c_str();
323 // error = mb->write_file(output_file, 0, write_options.c_str(), &lsets[l+1], 1);MB_CHK_ERR(error);
324 // mb->list_entity(lsets[l+1]);
325 // mb->write_file(output_file, 0, "PARALLEL=WRITE_PART", &lsets[l+1], 1);
326 }
327 }
328
329 delete uref;
330 #ifdef MOAB_HAVE_MPI
331 delete pc;
332 #endif
333
334 delete [] ldeg;
335 delete moab;
336
337 #ifdef MOAB_HAVE_MPI
338 MPI_Finalize();
339 #endif
340
341 exit(SUCCESS);
342 }
343
get_degree_seq(Core & mb,EntityHandle fileset,int dim,double desired_vol,int & num_levels,std::vector<int> & level_degs)344 ErrorCode get_degree_seq(Core &mb, EntityHandle fileset, int dim, double desired_vol, int &num_levels, std::vector<int> &level_degs)
345 {
346 //Find max volume
347 double vmax_global;
348 ErrorCode error = get_max_volume(mb, fileset, dim, vmax_global);MB_CHK_ERR(error);
349
350 int init_nl = num_levels;
351 num_levels = 0; level_degs.clear();
352
353 //Find a sequence that reduces the maximum volume to desired.
354 //desired_vol should be a fraction or absolute value ?
355
356 //double remV = vmax_global*desired_vol/dim;
357 double Vdesired = desired_vol;
358 double try_x;
359 double remV = vmax_global;
360 int degs[3][3] = {{5,3,2},{25,9,4},{0,27,8}};
361
362 if (dim==1 || dim == 2)
363 {
364 while (remV -Vdesired>= 0){
365 try_x = degs[dim-1][0];
366 if ((remV/try_x - Vdesired)>=0)
367 {
368 level_degs.push_back(5);
369 num_levels += 1;
370 remV /= try_x;
371 }
372 else
373 {
374 try_x = degs[dim-1][1];
375 if ((remV/try_x - Vdesired)>=0)
376 {
377 level_degs.push_back(3);
378 num_levels += 1;
379 remV /= try_x;
380 }
381 else {
382 try_x = degs[dim-1][2];
383 if ((remV/try_x - Vdesired)>=0)
384 {
385 level_degs.push_back(2);
386 num_levels += 1;
387 remV /= try_x;
388 }
389 else
390 break;
391 }
392 }
393 }
394 }
395 else
396 {
397 while (remV -Vdesired>= 0){
398 try_x = degs[dim-1][1];
399 if ((remV/try_x - Vdesired)>=0)
400 {
401 level_degs.push_back(3);
402 num_levels += 1;
403 remV /= try_x;
404 }
405 else
406 {
407 try_x = degs[dim-1][2];
408 if ((remV/try_x - Vdesired)>=0)
409 {
410 level_degs.push_back(2);
411 num_levels += 1;
412 remV /= try_x;
413 }
414 else
415 break;
416 }
417 }
418 }
419
420 if (init_nl != 0 && init_nl < num_levels)
421 {
422 for (int i=level_degs.size(); i>= init_nl; i--)
423 level_degs.pop_back();
424 num_levels = init_nl;
425 }
426
427 return MB_SUCCESS;
428 }
429
get_max_volume(Core & mb,EntityHandle fileset,int dim,double & vmax)430 ErrorCode get_max_volume(Core &mb, EntityHandle fileset, int dim, double &vmax)
431 {
432 ErrorCode error;
433 VerdictWrapper vw(&mb);
434 QualityType q;
435
436 switch (dim) {
437 case 1: q = MB_LENGTH; break;
438 case 2: q = MB_AREA; break;
439 case 3: q = MB_VOLUME; break;
440 default: return MB_FAILURE; break;
441 }
442
443 //Get all entities of the highest dimension which is passed as a command line argument.
444 Range allents, owned;
445 error = mb.get_entities_by_handle(fileset, allents);MB_CHK_ERR(error);
446 owned = allents.subset_by_dimension(dim);MB_CHK_ERR(error);
447
448 //Get all owned entities
449 #ifdef MOAB_HAVE_MPI
450 int size = 1;
451 MPI_Comm_size( MPI_COMM_WORLD, &size );
452 int mpi_err;
453 if (size>1)
454 {
455 // filter the entities not owned, so we do not process them more than once
456 ParallelComm* pcomm = moab::ParallelComm::get_pcomm(&mb, 0);
457 Range current = owned;
458 owned.clear();
459 error = pcomm->filter_pstatus(current, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned);
460 if (error != MB_SUCCESS)
461 {
462 MPI_Finalize();
463 return MB_FAILURE;
464 }
465 }
466 #endif
467
468 double vmax_local=0;
469 //Find the maximum volume of an entity in the owned mesh
470 for (Range::iterator it=owned.begin(); it != owned.end(); it++)
471 {
472 double volume;
473 error = vw.quality_measure(*it, q, volume);MB_CHK_ERR(error);
474 if (volume >vmax_local)
475 vmax_local = volume;
476 }
477
478 //Get the global maximum
479 double vmax_global = vmax_local;
480 #ifdef MOAB_HAVE_MPI
481 mpi_err = MPI_Reduce(&vmax_local, &vmax_global, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
482 if (mpi_err)
483 {
484 MPI_Finalize();
485 return MB_FAILURE;
486 }
487 #endif
488
489 vmax = vmax_global;
490
491 return MB_SUCCESS;
492 }
493
494
parse_id_list(const char * string,int dim,int nval,std::vector<int> & results)495 bool parse_id_list( const char* string, int dim, int nval, std::vector<int> &results )
496 {
497 bool okay = true;
498 char* mystr = strdup( string );
499 for (const char* ptr = strtok(mystr, ","); ptr; ptr = strtok(0,","))
500 {
501 char* endptr;
502 int val = strtol( ptr, &endptr, 0 );
503
504 if (dim==1 || dim == 2){
505 if(val != 2 && val != 3 && val != 5){
506 std::cerr <<"Not a valid degree for the passed dimension" << std::endl;
507 okay = false;
508 break;
509 }
510 }
511 else if (dim==3){
512 if(val != 2 && val != 3){
513 std::cerr <<"Not a valid degree for the passed dimension" << std::endl;
514 okay = false;
515 break;
516 }
517 }
518
519 if (endptr == ptr || val <= 0) {
520 std::cerr << "Not a valid id: " << ptr << std::endl;
521 okay = false;
522 break;
523 }
524
525 results.push_back(val);
526 }
527
528 if ((int)results.size() < nval)
529 {
530 for (int i=results.size(); i<=nval-1; i++)
531 results.push_back(results[0]);
532 }
533 else if ((int)results.size() > nval)
534 {
535 for (int i=results.size(); i>nval; i--)
536 results.pop_back();
537 }
538
539 free( mystr );
540 return okay;
541 }
542
make_opts_string(std::vector<std::string> options,std::string & opts)543 bool make_opts_string( std::vector<std::string> options, std::string& opts )
544 {
545 opts.clear();
546 if (options.empty())
547 return true;
548
549 // choose a separator character
550 std::vector<std::string>::const_iterator i;
551 char separator = '\0';
552 const char* alt_separators = ";+,:\t\n";
553 for (const char* sep_ptr = alt_separators; *sep_ptr; ++sep_ptr) {
554 bool seen = false;
555 for (i = options.begin(); i != options.end(); ++i)
556 if (i->find( *sep_ptr, 0 ) != std::string::npos) {
557 seen = true;
558 break;
559 }
560 if (!seen) {
561 separator = *sep_ptr;
562 break;
563 }
564 }
565 if (!separator) {
566 std::cerr << "Error: cannot find separator character for options string" << std::endl;
567 return false;
568 }
569 if (separator != ';') {
570 opts = ";";
571 opts += separator;
572 }
573
574 // concatenate options
575 i = options.begin();
576 opts += *i;
577 for (++i; i != options.end(); ++i) {
578 opts += separator;
579 opts += *i;
580 }
581
582 return true;
583 }
584
585