1
2 #include "Mesquite.hpp"
3 #include "MsqIBase.hpp"
4 #include "MsqIGeom.hpp"
5 #include "MsqIMesh.hpp"
6 #include "MBiMesh.hpp"
7 #include "MeshImpl.hpp"
8 #include "moab/Core.hpp"
9 #include "moab/Skinner.hpp"
10 #include "moab/LloydSmoother.hpp"
11 #include "FacetModifyEngine.hpp"
12 #include "MsqError.hpp"
13 #include "InstructionQueue.hpp"
14 #include "PatchData.hpp"
15 #include "TerminationCriterion.hpp"
16 #include "QualityAssessor.hpp"
17 #include "SphericalDomain.hpp"
18 #include "PlanarDomain.hpp"
19 #include "MeshWriter.hpp"
20
21 #include "moab/Core.hpp"
22 #include "moab/ProgOptions.hpp"
23 #ifdef MOAB_HAVE_MPI
24 #include "moab/ParallelComm.hpp"
25 #endif
26
27 // algorithms
28 #include "IdealWeightInverseMeanRatio.hpp"
29 #include "TMPQualityMetric.hpp"
30 #include "AspectRatioGammaQualityMetric.hpp"
31 #include "ConditionNumberQualityMetric.hpp"
32 #include "VertexConditionNumberQualityMetric.hpp"
33 #include "MultiplyQualityMetric.hpp"
34 #include "LPtoPTemplate.hpp"
35 #include "LInfTemplate.hpp"
36 #include "PMeanPTemplate.hpp"
37 #include "SteepestDescent.hpp"
38 #include "FeasibleNewton.hpp"
39 #include "QuasiNewton.hpp"
40 #include "ConjugateGradient.hpp"
41 #include "SmartLaplacianSmoother.hpp"
42 #include "Randomize.hpp"
43
44 #include "ReferenceMesh.hpp"
45 #include "RefMeshTargetCalculator.hpp"
46 #include "TShapeB1.hpp"
47 #include "TQualityMetric.hpp"
48 #include "IdealShapeTarget.hpp"
49
50 #include <sys/stat.h>
51 #include <iostream>
52 using std::cerr;
53 using std::cout;
54 using std::endl;
55
56 #include "iBase.h"
57
58 using namespace MBMesquite;
59
print_iGeom_error(const char * desc,int err,iGeom_Instance geom,const char * file,int line)60 static int print_iGeom_error( const char* desc,
61 int err,
62 iGeom_Instance geom,
63 const char* file,
64 int line )
65 {
66 char buffer[1024];
67 iGeom_getDescription( geom, buffer, sizeof(buffer) );
68 buffer[sizeof(buffer)-1] = '\0';
69
70 std::cerr << "ERROR: " << desc << std::endl
71 << " Error code: " << err << std::endl
72 << " Error desc: " << buffer << std::endl
73 << " At : " << file << ':' << line << std::endl
74 ;
75
76 return -1; // must always return false or CHECK macro will break
77 }
78
print_iMesh_error(const char * desc,int err,iMesh_Instance mesh,const char * file,int line)79 static int print_iMesh_error( const char* desc,
80 int err,
81 iMesh_Instance mesh,
82 const char* file,
83 int line )
84 {
85 char buffer[1024];
86 iMesh_getDescription( mesh, buffer, sizeof(buffer) );
87 buffer[sizeof(buffer)-1] = '\0';
88
89 std::cerr << "ERROR: " << desc << std::endl
90 << " Error code: " << err << std::endl
91 << " Error desc: " << buffer << std::endl
92 << " At : " << file << ':' << line << std::endl
93 ;
94
95 return -1; // must always return false or CHECK macro will break
96 }
97
98 #define CHECK_IGEOM( STR ) if (err != iBase_SUCCESS) return print_iGeom_error( STR, err, geom, __FILE__, __LINE__ )
99
100 #define CHECK_IMESH( STR ) if (err != iBase_SUCCESS) return print_iMesh_error( STR, err, instance, __FILE__, __LINE__ )
101
102 const std::string default_file_name = std::string(MESH_DIR) + std::string("/surfrandomtris-4part.h5m");
103 const std::string default_geometry_file_name = std::string(MESH_DIR) + std::string("/surfrandom.facet");
104
105 std::vector<double> solution_indicator;
106
107 int write_vtk_mesh( Mesh* mesh, const char* filename);
108
109 // Construct a MeshTSTT from the file
110 int get_imesh_mesh( MBMesquite::Mesh** , const char* file_name, int dimension );
111
112 // Construct a MeshImpl from the file
113 int get_native_mesh( MBMesquite::Mesh** , const char* file_name, int dimension );
114
115 int get_itaps_domain(MeshDomain**, const char*);
116 int get_mesquite_domain(MeshDomain**, const char*);
117
118 // Run FeasibleNewton solver
119 int run_global_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value=1e-4 );
120
121 // Run SmoothLaplacian solver
122 int run_local_smoother( MeshDomainAssoc& mesh, MsqError& err, double OF_value=1e-3 );
123 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value=1e-3 );
124
125 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
126
127 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err );
128
file_exists(const std::string & name)129 bool file_exists (const std::string& name) {
130 struct stat buffer;
131 return (stat (name.c_str(), &buffer) == 0);
132 }
133
main(int argc,char * argv[])134 int main(int argc, char* argv[])
135 {
136 MBMesquite::MsqPrintError err(cout);
137 // command line arguments
138 std::string file_name, geometry_file_name;
139 bool use_native = false;
140 int dimension = 2;
141
142 #ifdef MOAB_HAVE_MPI
143 // MPI_Init(&argc, &argv);
144 #endif
145 ProgOptions opts;
146
147 opts.addOpt<void>(std::string("native,N"),
148 std::string("Use native representation (default=false)"), &use_native);
149 opts.addOpt<int>(std::string("dim,d"),
150 std::string("Topological dimension of the mesh (default=2)"), &dimension);
151 opts.addOpt<std::string>(std::string("input_geo,i"),
152 std::string("Input file name for the geometry (pattern=*.stl, *.facet)"), &geometry_file_name);
153 opts.addOpt<std::string>(std::string("input_mesh,m"),
154 std::string("Input file name for the mesh (pattern=*.vtk, *.h5m)"), &file_name);
155
156 opts.parseCommandLine(argc, argv);
157
158 if (!geometry_file_name.length())
159 {
160 file_name = default_file_name;
161 geometry_file_name = default_geometry_file_name;
162 cout << "No file specified: Using defaults.\n";
163 }
164 cout << "\t Mesh filename = " << file_name << endl;
165 cout << "\t Geometry filename = " << geometry_file_name << endl;
166
167 // Vector3D pnt(0,0,0);
168 // Vector3D s_norm(0,0,1);
169 // PlanarDomain plane(s_norm, pnt);
170
171 // PlanarDomain plane( PlanarDomain::XY );
172 MeshDomain* plane;
173 int ierr;
174 if (!file_exists(geometry_file_name)) geometry_file_name = "";
175 ierr = get_itaps_domain(&plane, geometry_file_name.c_str());//MB_CHK_ERR(ierr);
176
177 // Try running a global smoother on the mesh
178 Mesh* mesh=0;
179 if(use_native) {
180 ierr = get_native_mesh(&mesh, file_name.c_str(), dimension);//MB_CHK_ERR(ierr);
181 }
182 else {
183 ierr = get_imesh_mesh(&mesh, file_name.c_str(), dimension);//MB_CHK_ERR(ierr);
184 }
185
186 if (!mesh) {
187 std::cerr << "Failed to load input file. Aborting." << std::endl;
188 return 1;
189 }
190
191 MeshDomainAssoc mesh_and_domain(mesh, plane);
192
193 // ierr = write_vtk_mesh( mesh, "original.vtk");MB_CHK_ERR(ierr);
194 // cout << "Wrote \"original.vtk\"" << endl;
195
196 // run_global_smoother( mesh_and_domain, err );
197 // if (err) return 1;
198
199 // Try running a local smoother on the mesh
200 // Mesh* meshl=0;
201 // if(use_native)
202 // ierr = get_native_mesh(&meshl, file_name.c_str(), dimension);
203 // else
204 // ierr = get_imesh_mesh(&meshl, file_name.c_str(), dimension);
205 // if (!mesh || ierr) {
206 // std::cerr << "Failed to load input file. Aborting." << std::endl;
207 // return 1;
208 // }
209
210 // MeshDomainAssoc mesh_and_domain_local(meshl, plane);
211
212 // run_solution_mesh_optimizer( mesh_and_domain, err );
213 // if (err) return 1;
214
215 run_local_smoother( mesh_and_domain, err, 1e-4 );//MB_CHK_ERR(err);
216 if (err) return 1;
217
218 double reps = 5e-2;
219 for (int iter=0; iter < 5; iter++) {
220
221 if (!(iter % 2)) {
222 run_local_smoother2( mesh_and_domain, err, reps*10 );//CHECK_IMESH("local smoother2 failed");
223 if (err) return 1;
224 }
225
226 // run_global_smoother( mesh_and_domain, err, reps );MB_CHK_ERR(ierr);
227
228 run_solution_mesh_optimizer( mesh_and_domain, err );//CHECK_IMESH("solution mesh optimizer failed");
229 if (err) return 1;
230
231 reps *= 0.01;
232 }
233
234 run_local_smoother2( mesh_and_domain, err, 1e-4 );//CHECK_IMESH("local smoother2 failed");
235 if (err) return 1;
236
237 // run_quality_optimizer( mesh_and_domain, err );MB_CHK_ERR(ierr);
238
239 // run_local_smoother( mesh_and_domain, err );MB_CHK_ERR(ierr);
240
241 // Delete MOAB instance
242 delete mesh;
243 delete plane;
244
245 #ifdef MOAB_HAVE_MPI
246 //MPI_Finalize();
247 #endif
248
249 return 0;
250 }
251
252
run_global_smoother(MeshDomainAssoc & mesh_and_domain,MsqError & err,double OF_value)253 int run_global_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
254 {
255 // double OF_value = 1e-6;
256
257 Mesh* mesh = mesh_and_domain.get_mesh();
258 MeshDomain* domain = mesh_and_domain.get_domain();
259
260 // creates an intruction queue
261 InstructionQueue queue1;
262
263 // creates a mean ratio quality metric ...
264 IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
265 // ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
266 // TMPQualityMetric* mean_ratio = new TMPQualityMetric();
267
268 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
269 if (err) return 1;
270 mean_ratio->set_averaging_method(QualityMetric::RMS, err);
271 if (err) return 1;
272
273 // ... and builds an objective function with it
274 LPtoPTemplate* obj_func = new LPtoPTemplate(mean_ratio, 1, err);
275 // LInfTemplate* obj_func = new LInfTemplate(mean_ratio);
276 if (err) return 1;
277
278 // creates the feas newt optimization procedures
279 // ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
280 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
281 SteepestDescent* pass1 = new SteepestDescent( obj_func );
282 pass1->use_element_on_vertex_patch();
283 pass1->do_block_coordinate_descent_optimization();
284 pass1->use_global_patch();
285 if (err) return 1;
286
287 QualityAssessor stop_qa( mean_ratio );
288
289 // **************Set stopping criterion****************
290 TerminationCriterion tc_inner;
291 tc_inner.add_absolute_vertex_movement( OF_value );
292 if (err) return 1;
293 TerminationCriterion tc_outer;
294 tc_inner.add_iteration_limit( 10 );
295 tc_outer.add_iteration_limit( 5 );
296 tc_outer.set_debug_output_level(3);
297 tc_inner.set_debug_output_level(3);
298 pass1->set_inner_termination_criterion(&tc_inner);
299 pass1->set_outer_termination_criterion(&tc_outer);
300
301 queue1.add_quality_assessor(&stop_qa, err);
302 if (err) return 1;
303
304 // adds 1 pass of pass1 to mesh_set1
305 queue1.set_master_quality_improver(pass1, err);
306 if (err) return 1;
307
308 queue1.add_quality_assessor(&stop_qa, err);
309 if (err) return 1;
310
311 // launches optimization on mesh_set
312 if (domain) {
313 queue1.run_instructions(&mesh_and_domain, err);
314 }
315 else {
316 queue1.run_instructions(mesh, err);
317 }
318 if (err) return 1;
319
320
321 // Construct a MeshTSTT from the file
322 int ierr = write_vtk_mesh( mesh, "feasible-newton-result.vtk");
323 if (ierr) return 1;
324 // MeshWriter::write_vtk(mesh, "feasible-newton-result.vtk", err);
325 // if (err) return 1;
326 cout << "Wrote \"feasible-newton-result.vtk\"" << endl;
327
328 //print_timing_diagnostics( cout );
329 return 0;
330 }
331
write_vtk_mesh(Mesh * mesh,const char * filename)332 int write_vtk_mesh( Mesh* mesh, const char* filename )
333 {
334 moab::Interface* mbi = reinterpret_cast<MBiMesh*>(dynamic_cast<MsqIMesh*>(mesh)->get_imesh_instance())->mbImpl;
335
336 mbi->write_file(filename);
337
338 return 0;
339 }
340
341 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value );
run_local_smoother(MeshDomainAssoc & mesh_and_domain,MsqError & err,double OF_value)342 int run_local_smoother( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
343 {
344 Mesh* mesh = mesh_and_domain.get_mesh();
345 moab::Interface* mbi = reinterpret_cast<MBiMesh*>(dynamic_cast<MsqIMesh*>(mesh)->get_imesh_instance())->mbImpl;
346
347
348 moab::Tag fixed;
349 moab::ErrorCode rval = mbi->tag_get_handle("fixed", 1, moab::MB_TYPE_INTEGER, fixed); MB_CHK_SET_ERR(rval, "Getting tag handle failed");
350 moab::Range cells;
351 rval = mbi->get_entities_by_dimension(0, 2, cells); MB_CHK_SET_ERR(rval, "Querying elements failed");
352
353 moab::LloydSmoother lloyd(mbi, 0, cells, 0, 0/*fixed*/);
354
355 lloyd.perform_smooth();
356
357 // run_local_smoother2(mesh_and_domain, err, OF_value);
358
359 // Construct a MeshTSTT from the file
360 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk");
361 if (ierr) return 1;
362 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
363 // if (err) return 1;
364 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
365 return 0;
366 }
367
run_local_smoother2(MeshDomainAssoc & mesh_and_domain,MsqError & err,double OF_value)368 int run_local_smoother2( MeshDomainAssoc& mesh_and_domain, MsqError& err, double OF_value )
369 {
370 // double OF_value = 1e-5;
371
372 Mesh* mesh = mesh_and_domain.get_mesh();
373 MeshDomain* domain = mesh_and_domain.get_domain();
374
375 // creates an intruction queue
376 InstructionQueue queue1;
377
378 // creates a mean ratio quality metric ...
379 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
380 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
381 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
382 if (err) return 1;
383 // mean_ratio->set_gradient_type(QualityMetric::NUMERICAL_GRADIENT);
384 // mean_ratio->set_hessian_type(QualityMetric::NUMERICAL_HESSIAN);
385 mean_ratio->set_averaging_method(QualityMetric::RMS, err);
386 if (err) return 1;
387
388 // ... and builds an objective function with it
389 LPtoPTemplate* obj_func = new LPtoPTemplate(mean_ratio, 1, err);
390 if (err) return 1;
391
392 if (false)
393 {
394 InstructionQueue qtmp;
395 Randomize rand (-0.005) ;
396 TerminationCriterion sc_rand;
397 sc_rand.add_iteration_limit( 2 ) ;
398 rand.set_outer_termination_criterion(&sc_rand);
399 qtmp.set_master_quality_improver(&rand, err);
400 if (err) return 1;
401 if (domain) {
402 qtmp.run_instructions(&mesh_and_domain, err);
403 }
404 else {
405 qtmp.run_instructions(mesh, err);
406 }
407 if (err) return 1;
408 }
409
410 // creates the smart laplacian optimization procedures
411 SmartLaplacianSmoother* pass1 = new SmartLaplacianSmoother( obj_func );
412 // SteepestDescent* pass1 = new SteepestDescent( obj_func );
413
414 QualityAssessor stop_qa( mean_ratio );
415
416 // **************Set stopping criterion****************
417 TerminationCriterion tc_inner;
418 tc_inner.add_absolute_vertex_movement( OF_value );
419 TerminationCriterion tc_outer;
420 tc_outer.add_iteration_limit( 10 );
421 pass1->set_inner_termination_criterion(&tc_inner);
422 pass1->set_outer_termination_criterion(&tc_outer);
423
424 queue1.add_quality_assessor(&stop_qa, err);
425 if (err) return 1;
426
427 // adds 1 pass of pass1 to mesh_set
428 queue1.set_master_quality_improver(pass1, err);
429 if (err) return 1;
430
431 queue1.add_quality_assessor(&stop_qa, err);
432 if (err) return 1;
433
434 // launches optimization on mesh_set
435 if (domain) {
436 queue1.run_instructions(&mesh_and_domain, err);
437 }
438 else {
439 queue1.run_instructions(mesh, err);
440 }
441 if (err) return 1;
442
443 // Construct a MeshTSTT from the file
444 int ierr = write_vtk_mesh( mesh, "smart-laplacian-result.vtk");
445 if (ierr) return 1;
446 // MeshWriter::write_vtk(mesh, "smart-laplacian-result.vtk", err);
447 // if (err) return 1;
448 cout << "Wrote \"smart-laplacian-result.vtk\"" << endl;
449
450 //print_timing_diagnostics( cout );
451 return 0;
452 }
453
454
run_quality_optimizer(MeshDomainAssoc & mesh_and_domain,MsqError & err)455 int run_quality_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
456 {
457 Mesh* mesh = mesh_and_domain.get_mesh();
458 MeshDomain* domain = mesh_and_domain.get_domain();
459
460 // creates an intruction queue
461 InstructionQueue q;
462
463 // do optimization
464 const double eps = 0.01;
465 IdealShapeTarget w;
466 TShapeB1 mu;
467 TQualityMetric metric( &w, &mu );
468 PMeanPTemplate func( 1.0, &metric );
469
470 SteepestDescent solver( &func );
471 solver.use_element_on_vertex_patch();
472 solver.do_jacobi_optimization();
473
474 TerminationCriterion inner, outer;
475 inner.add_absolute_vertex_movement( 0.5*eps );
476 outer.add_absolute_vertex_movement( 0.5*eps );
477
478 QualityAssessor qa( &metric );
479
480 q.add_quality_assessor( &qa, err );
481 if (err) return 1;
482 q.set_master_quality_improver( &solver, err );
483 if (err) return 1;
484 q.add_quality_assessor( &qa, err );
485 if (err) return 1;
486
487 // launches optimization on mesh_set
488 if (domain) {
489 q.run_instructions(&mesh_and_domain, err);
490 }
491 else {
492 q.run_instructions(mesh, err);
493 }
494 if (err) return 1;
495
496 // Construct a MeshTSTT from the file
497 int ierr = write_vtk_mesh( mesh, "ideal-shape-result.vtk");
498 if (ierr) return 1;
499 // MeshWriter::write_vtk(mesh, "ideal-shape-result.vtk", err);
500 // if (err) return 1;
501 cout << "Wrote \"ideal-shape-result.vtk\"" << endl;
502
503 print_timing_diagnostics( cout );
504 return 0;
505 }
506
507
run_solution_mesh_optimizer(MeshDomainAssoc & mesh_and_domain,MsqError & err)508 int run_solution_mesh_optimizer( MeshDomainAssoc& mesh_and_domain, MsqError& err )
509 {
510 double OF_value = 0.01;
511
512 Mesh* mesh = mesh_and_domain.get_mesh();
513 MeshDomain* domain = mesh_and_domain.get_domain();
514
515 // creates an intruction queue
516 InstructionQueue queue1;
517
518 // creates a mean ratio quality metric ...
519 // IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
520 ConditionNumberQualityMetric* mean_ratio = new ConditionNumberQualityMetric();
521 // VertexConditionNumberQualityMetric* mean_ratio = new VertexConditionNumberQualityMetric();
522 // AspectRatioGammaQualityMetric* mean_ratio = new AspectRatioGammaQualityMetric();
523
524 //ElementSolIndQM* solindqm = new ElementSolIndQM(solution_indicator);
525 //MultiplyQualityMetric* mean_ratio = new MultiplyQualityMetric(new VertexConditionNumberQualityMetric(), solindqm, err);
526 // ElementSolIndQM* mean_ratio = solindqm;
527
528 // LocalSizeQualityMetric* mean_ratio = new LocalSizeQualityMetric();
529
530 mean_ratio->set_averaging_method(QualityMetric::SUM_OF_RATIOS_SQUARED, err);
531 if (err) return 1;
532
533 // ... and builds an objective function with it
534 LPtoPTemplate* obj_func = new LPtoPTemplate(mean_ratio, 1, err);
535 if (err) return 1;
536
537 // // creates the feas newt optimization procedures
538 ConjugateGradient* pass1 = new ConjugateGradient( obj_func, err );
539 // QuasiNewton* pass1 = new QuasiNewton( obj_func );
540 // FeasibleNewton* pass1 = new FeasibleNewton( obj_func );
541 pass1->use_global_patch();
542
543 QualityAssessor stop_qa( mean_ratio );
544
545 // **************Set stopping criterion****************
546 TerminationCriterion tc_inner;
547 tc_inner.add_absolute_vertex_movement( OF_value );
548 if (err) return 1;
549 TerminationCriterion tc_outer;
550 tc_inner.add_iteration_limit( 20 );
551 tc_outer.add_iteration_limit( 5 );
552 pass1->set_inner_termination_criterion(&tc_inner);
553 pass1->set_outer_termination_criterion(&tc_outer);
554 pass1->set_debugging_level(3);
555
556 queue1.add_quality_assessor(&stop_qa, err);
557 if (err) return 1;
558
559 // adds 1 pass of pass1 to mesh_set1
560 queue1.set_master_quality_improver(pass1, err);
561 if (err) return 1;
562
563 queue1.add_quality_assessor(&stop_qa, err);
564 if (err) return 1;
565
566 // launches optimization on mesh_set
567 if (domain) {
568 queue1.run_instructions(&mesh_and_domain, err);
569 }
570 else {
571 queue1.run_instructions(mesh, err);
572 }
573 if (err) return 1;
574
575 // Construct a MeshTSTT from the file
576 int ierr = write_vtk_mesh( mesh, "solution-mesh-result.vtk");
577 if (ierr) return 1;
578 // MeshWriter::write_vtk(mesh, "solution-mesh-result.vtk", err);
579 // if (err) return 1;
580 cout << "Wrote \"solution-mesh-result.vtk\"" << endl;
581
582 print_timing_diagnostics( cout );
583 return 0;
584 }
585
586
587
get_imesh_mesh(MBMesquite::Mesh ** mesh,const char * file_name,int dimension)588 int get_imesh_mesh( MBMesquite::Mesh** mesh, const char* file_name, int dimension )
589 {
590 int err;
591 iMesh_Instance instance = 0;
592 iMesh_newMesh( NULL, &instance, &err, 0 ); CHECK_IMESH("Creation of mesh instance failed");
593
594 iBase_EntitySetHandle root_set;
595 iMesh_getRootSet( instance, &root_set, &err );CHECK_IMESH("Could not get root set");
596
597 iMesh_load( instance, root_set, file_name, 0, &err, strlen(file_name), 0 );CHECK_IMESH("Could not load mesh from file");
598
599 iBase_TagHandle fixed_tag;
600 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen("fixed") );
601 // if (iBase_SUCCESS != err)
602 {
603 // get the skin vertices of those cells and mark them as fixed; we don't want to fix the vertices on a
604 // part boundary, but since we exchanged a layer of ghost cells, those vertices aren't on the skin locally
605 // ok to mark non-owned skin vertices too, I won't move those anyway
606 // use MOAB's skinner class to find the skin
607
608 // get all vertices and cells
609 // make tag to specify fixed vertices, since it's input to the algorithm; use a default value of non-fixed
610 // so we only need to set the fixed tag for skin vertices
611 moab::Interface* mbi = reinterpret_cast<MBiMesh*>(instance)->mbImpl;
612 moab::EntityHandle currset=0;
613 moab::Tag fixed;
614 int def_val = 0;
615 err=0;
616 moab::ErrorCode rval = mbi->tag_get_handle("fixed", 1, moab::MB_TYPE_INTEGER, fixed, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val); MB_CHK_SET_ERR(rval, "Getting tag handle failed");
617 moab::Range verts, cells, skin_verts;
618 rval = mbi->get_entities_by_type(currset, moab::MBVERTEX, verts); MB_CHK_SET_ERR(rval, "Querying vertices failed");
619 rval = mbi->get_entities_by_dimension(currset, dimension, cells); MB_CHK_SET_ERR(rval, "Querying elements failed");
620 std::cout << "Found " << verts.size() << " vertices and " << cells.size() << " elements" << std::endl;
621
622 moab::Skinner skinner(mbi);
623 rval = skinner.find_skin(currset, cells, true, skin_verts); MB_CHK_SET_ERR(rval, "Finding the skin of the mesh failed"); // 'true' param indicates we want vertices back, not cells
624
625 std::vector<int> fix_tag(skin_verts.size(), 1); // initialized to 1 to indicate fixed
626 rval = mbi->tag_set_data(fixed, skin_verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Setting tag data failed");
627 std::cout << "Found " << skin_verts.size() << " vertices on the skin of the domain." << std::endl;
628
629 // fix_tag.resize(verts.size(),0);
630 // rval = mbi->tag_get_data(fixed, verts, &fix_tag[0]); MB_CHK_SET_ERR(rval, "Getting tag data failed");
631
632 iMesh_getTagHandle( instance, "fixed", &fixed_tag, &err, strlen("fixed") );CHECK_IMESH("Getting tag handle (fixed) failed");
633
634 // Set some arbitrary solution indicator
635 moab::Tag solindTag;
636 double def_val_dbl=0.0;
637 rval = mbi->tag_get_handle("solution_indicator", 1, moab::MB_TYPE_DOUBLE, solindTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE, &def_val_dbl); MB_CHK_SET_ERR(rval, "Getting tag handle failed");
638 solution_indicator.resize(cells.size(),0.01);
639 for (unsigned i=0; i < cells.size()/4; i++)
640 solution_indicator[i]=0.1;
641 for (unsigned i=cells.size()/4; i < 2*cells.size()/4; i++)
642 solution_indicator[i]=0.5;
643 for (unsigned i=2*cells.size()/4; i < 3*cells.size()/4; i++)
644 solution_indicator[i]=0.5;
645 for (unsigned i=3*cells.size()/4; i < cells.size(); i++)
646 solution_indicator[i]=0.5;
647
648 rval = mbi->tag_set_data(solindTag, cells, &solution_indicator[0]); MB_CHK_SET_ERR(rval, "Setting tag data failed");
649
650 }
651
652
653 MsqError ierr;
654 MBMesquite::MsqIMesh* imesh = new MBMesquite::MsqIMesh( instance, root_set, (dimension == 3 ? iBase_REGION : iBase_FACE), ierr, &fixed_tag );
655 if (MSQ_CHKERR(ierr)) {
656 delete imesh;
657 cerr << err << endl;
658 err = iBase_FAILURE;
659 CHECK_IMESH("Creation of MsqIMesh instance failed");
660 return 0;
661 }
662
663 *mesh = imesh;
664 return iBase_SUCCESS;
665 }
666
667
668
get_native_mesh(MBMesquite::Mesh ** mesh,const char * file_name,int)669 int get_native_mesh( MBMesquite::Mesh** mesh, const char* file_name, int )
670 {
671 MsqError err;
672 MBMesquite::MeshImpl* imesh = new MBMesquite::MeshImpl();
673 imesh->read_vtk( file_name, err );
674 if (err)
675 {
676 cerr << err << endl;
677 exit(3);
678 }
679 *mesh = imesh;
680
681 return iBase_SUCCESS;
682 }
683
684
get_itaps_domain(MeshDomain ** odomain,const char * filename)685 int get_itaps_domain(MeshDomain** odomain, const char* filename)
686 {
687
688 if (filename == 0 || strlen(filename) == 0) {
689 *odomain=new PlanarDomain( PlanarDomain::XY );
690 return 0;
691 }
692
693 int err;
694 iGeom_Instance geom;
695 iGeom_newGeom( "", &geom, &err, 0 ); CHECK_IGEOM("ERROR: iGeom creation failed");
696
697 #ifdef MOAB_HAVE_CGM_FACET
698 FacetModifyEngine::set_modify_enabled(CUBIT_TRUE);
699 #endif
700
701 iGeom_load( geom, filename, 0, &err, strlen(filename), 0 );
702 CHECK_IGEOM( "Cannot load the geometry" );
703
704 iBase_EntitySetHandle root_set;
705 iGeom_getRootSet( geom, &root_set, &err );
706 CHECK_IGEOM( "getRootSet failed!" );
707
708 // print out the number of entities
709 std::cout << "Model contents: " << std::endl;
710 const char *gtype[] = {"vertices: ", "edges: ", "faces: ", "regions: "};
711 int nents[4];
712 for (int i = 0; i <= 3; ++i) {
713 iGeom_getNumOfType( geom, root_set, i, &nents[i], &err );
714 CHECK_IGEOM( "Error: problem getting entities after gLoad." );
715 std::cout << gtype[i] << nents[i] << std::endl;
716 }
717
718 iBase_EntityHandle* hd_geom_ents;
719 int csize=0, sizealloc=0;
720 if (nents[3] > 0) {
721 hd_geom_ents = (iBase_EntityHandle*)malloc(sizeof(iBase_EntityHandle)*nents[2]);
722 csize = nents[2];
723 iGeom_getEntities(geom, root_set, 2, &hd_geom_ents, &csize, &sizealloc, &err);
724 }
725 else {
726 hd_geom_ents = (iBase_EntityHandle*)malloc(sizeof(iBase_EntityHandle)*nents[1]);
727 csize = nents[1];
728 iGeom_getEntities(geom, root_set, 1, &hd_geom_ents, &csize, &sizealloc, &err);
729 }
730 CHECK_IGEOM( "ERROR: Could not get entities" );
731
732 *odomain = new MsqIGeom( geom, hd_geom_ents[0] );
733 return iBase_SUCCESS;
734 }
735
736