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