1 /* *****************************************************************
2 MESQUITE -- The Mesh Quality Improvement Toolkit
3
4 Copyright 2010 Sandia National Laboratories. Developed at the
5 University of Wisconsin--Madison under SNL contract number
6 624796. The U.S. Government and the University of Wisconsin
7 retain certain rights to this software.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 (lgpl.txt) along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 (2011) kraftche@cae.wisc.edu
24
25 ***************************************************************** */
26
27
28 /** \file voshell.cpp
29 * \brief Implement some of the examples from N. Voshell.
30 * \author Jason Kraftcheck
31 * \author Nick Voshell
32 */
33
34 #include "TestUtil.hpp"
35 #include "ShapeImprover.hpp"
36 #include "UntangleWrapper.hpp"
37 #include "MeshImpl.hpp"
38 #include "MsqError.hpp"
39 #include "QualityAssessor.hpp"
40
41 #include "PlanarDomain.hpp"
42 #include "CylinderDomain.hpp"
43 #include "SphericalDomain.hpp"
44 #include "MeshDomain1D.hpp"
45 #include "DomainClassifier.hpp"
46
47 #include "HexLagrangeShape.hpp"
48 #include "TestUtil.hpp"
49
50 #include <iostream>
51 #include <string>
52 #include <cstdlib>
53
54 using namespace MBMesquite;
55
56 std::string get_homogenious_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
57
58 std::string get_part_example_tri( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
59
60 std::string get_part_example_quad( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
61
62 std::string get_sphere_cube_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
63
64 std::string get_cut_cube_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
65
66 std::string get_sphere_cylinder_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
67
68 std::string get_hex_3d_part_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
69
70
71 typedef std::string (*example_setup_t)( DomainClassifier& geom, MeshImpl& mesh, MsqError& err );
72 struct Example {
73 char flag;
74 example_setup_t func;
75 const char* desc;
76 };
77
78 int run_example( const Example& e, bool write_output_file );
79
80
81 const Example examples[] = {
82 { 'H', &get_homogenious_example, "homogenous mesh (curvey surface) example" },
83 { 't', &get_part_example_tri, "part mesh with triangles" },
84 { 'q', &get_part_example_quad, "part mesh with quads" },
85 { 'c', &get_sphere_cube_example, "cube with subtracted hemisphere" },
86 { 'l', &get_cut_cube_example, "cube with slot removed" },
87 { 's', &get_sphere_cylinder_example, "cylinder with subtracted sphere" },
88 { 'x', &get_hex_3d_part_example, "rectangular prism with two cylindrical slots" }
89 };
90 const int num_examples = sizeof(examples)/sizeof(examples[0]);
91
92
93
usage(const char * argv0,bool help=false)94 void usage( const char* argv0, bool help = false )
95 {
96 std::cerr << "Usage: " << argv0 << "[-w] ";
97 for (int i = 0; i < num_examples; ++i)
98 std::cerr << " [-" << examples[i].flag << "]";
99 std::cerr << std::endl;
100
101 if (!help) {
102 std::cerr << " " << argv0 << " -h" << std::endl;
103 std::exit(1);
104 }
105 std::cerr << " -w : Write final meshes to VTK file" << std::endl;
106 std::cerr << "NOTE: If no options are specified, all examples will be run" << std::endl;
107 for (int i = 0; i < num_examples; ++i)
108 std::cerr << " -" << examples[i].flag << " : Run " << examples[i].desc << std::endl;
109 std::cerr << std::endl;
110 std::exit(0);
111 }
112
113
114
main(int argc,char * argv[])115 int main( int argc, char* argv[] )
116 {
117 bool write_final_meshes = false;
118 std::vector<Example> list;
119 for (int i = 1; i < argc; ++i) {
120 if (argv[i][0] != '-') {
121 std::cerr << "Invalid argument: \"" << argv[i] << '"' << std::endl;
122 usage(argv[0]);
123 }
124 for (int j = 1; argv[i][j]; ++j) {
125 if (argv[i][j] == 'w') {
126 write_final_meshes = true;
127 }
128 else {
129 int idx;
130 for (idx = 0; idx < num_examples; ++idx)
131 if (examples[idx].flag == argv[i][j])
132 break;
133 if (idx == num_examples) {
134 if (argv[i][j] == 'h')
135 usage(argv[0],true);
136 else {
137 std::cerr << "Invalid flag: '" << argv[i][j] << "'" << std::endl;
138 usage(argv[0]);
139 }
140 }
141 list.push_back(examples[idx]);
142 }
143 }
144 }
145
146 const Example* exlist;
147 int exlistlen;
148 if (list.empty()) {
149 exlist = examples;
150 exlistlen = num_examples;
151 }
152 else {
153 exlist = &list[0];
154 exlistlen = list.size();
155 }
156
157 int result = 0;
158 for (int i = 0; i < exlistlen; ++i)
159 result += run_example( exlist[i], write_final_meshes );
160
161 return result;
162 }
163
164
165
run_example(const Example & e,bool write_output_file)166 int run_example( const Example& e, bool write_output_file )
167 {
168 MsqPrintError err(std::cerr);
169 MeshImpl mesh;
170 DomainClassifier domain;
171 HexLagrangeShape hex27;
172
173 std::cout << std::endl
174 << "--------------------------------------------------------------------"
175 << std::endl
176 << e.desc << std::endl
177 << "--------------------------------------------------------------------"
178 << std::endl;
179
180 std::string name = e.func( domain, mesh, err );
181 if (MSQ_CHKERR(err)) return 2;
182 std::cout << "Loaded mesh from: " << name << std::endl;
183
184 UntangleWrapper untangler;
185 untangler.set_slaved_ho_node_mode( Settings::SLAVE_NONE );
186 untangler.set_mapping_function( &hex27 );
187 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc(&mesh, &domain, false, true);
188 untangler.run_instructions( &mesh_and_domain, err );
189 if (MSQ_CHKERR(err)) return 1;
190 ShapeImprover smoother;
191 smoother.set_slaved_ho_node_mode( Settings::SLAVE_NONE );
192 smoother.set_mapping_function( &hex27 );
193 smoother.set_vertex_movement_limit_factor( 0.05 );
194 smoother.run_instructions( &mesh_and_domain, err );
195 if (MSQ_CHKERR(err)) return 1;
196
197 if (write_output_file) {
198 size_t idx = name.find( ".vtk" );
199 if (idx != std::string::npos) {
200 std::string newname( name.substr(0, idx) );
201 newname += ".out";
202 newname += name.substr(idx);
203 name.swap(newname);
204 }
205 else {
206 name += ".out";
207 }
208
209 mesh.write_vtk( name.c_str(), err ); MSQ_CHKERR(err);
210 std::cout << "Write mesh to file: " << name << std::endl;
211 }
212
213 return smoother.quality_assessor().invalid_elements();
214 }
215
get_planar_example(const char * filename,DomainClassifier & geom,MeshImpl & mesh,MsqError & err)216 void get_planar_example( const char* filename,
217 DomainClassifier& geom,
218 MeshImpl& mesh,
219 MsqError& err )
220 {
221 static PlanarDomain z(PlanarDomain::XY);
222 mesh.read_vtk(filename, err); MSQ_ERRRTN(err);
223 mesh.mark_skin_fixed(err); MSQ_ERRRTN(err);
224 DomainClassifier::DomainSet set(&z);
225 mesh.get_all_vertices( set.vertices, err ); MSQ_ERRRTN(err);
226 mesh.get_all_elements( set.elements, err ); MSQ_ERRRTN(err);
227 DomainClassifier::classify_by_handle( geom, &mesh, &set, 1, err ); MSQ_ERRRTN(err);
228 }
229
get_homogenious_example(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)230 std::string get_homogenious_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
231 {
232 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/homogeneousPart.vtk";
233 get_planar_example( filename.c_str(), geom, mesh, err ); MSQ_CHKERR(err);
234 return filename;
235 }
236
237
get_part_example_tri(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)238 std::string get_part_example_tri( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
239 {
240 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/triPart.vtk";
241 get_planar_example( filename.c_str(), geom, mesh, err ); MSQ_CHKERR(err);
242 return filename;
243 }
244
245
get_part_example_quad(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)246 std::string get_part_example_quad( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
247 {
248 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/quadPart.vtk";
249 get_planar_example( filename.c_str(), geom, mesh, err ); MSQ_CHKERR(err);
250 return filename;
251 }
252
253
get_sphere_cube_example(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)254 std::string get_sphere_cube_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
255 {
256 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/sphereCube.vtk";
257
258 const Vector3D vec_i(5,0,0), vec_ni(-5,0,0);
259 const Vector3D vec_j(0,5,0), vec_nj(0,-5,0);
260 const Vector3D vec_k(0,0,5), vec_nk(0,0,-5);
261
262 const Vector3D vec_tfl(5,-5,5), vec_tfr(5,5,5);
263 const Vector3D vec_tbl(-5,-5,5), vec_tbr(-5,5,5);
264 const Vector3D vec_bfl(5,-5,-5), vec_bfr(5,5,-5);
265 const Vector3D vec_bbl(-5,-5,-5), vec_bbr(-5,5,-5);
266
267 const Vector3D vec_ts(5,0,2), vec_bs(5,0,-2);
268
269
270 //++ 0D domains ++
271
272 static PointDomain pt_tfl(vec_tfl), pt_tfr(vec_tfr);
273 static PointDomain pt_tbl(vec_tbl), pt_tbr(vec_tbr);
274 static PointDomain pt_bfl(vec_bfl), pt_bfr(vec_bfr);
275 static PointDomain pt_bbl(vec_bbl), pt_bbr(vec_bbr);
276
277
278 //++ 1D domains ++
279
280 //top square
281 static LineDomain ln_tf(vec_tfl,vec_j), lin_tb(vec_tbr,vec_nj);
282 static LineDomain ln_tl(vec_tfl,vec_ni), lin_tr(vec_tbr,vec_i);
283 //bottom square
284 static LineDomain ln_bf(vec_bfl,vec_j), lin_bb(vec_bbr,vec_nj);
285 static LineDomain ln_bl(vec_bfl,vec_ni), lin_br(vec_bbr,vec_i);
286 //sides
287 static LineDomain ln_lf(vec_bfl,vec_k), lin_rf(vec_bfr,vec_k);
288 static LineDomain ln_lb(vec_bbl,vec_k), lin_rb(vec_bbr,vec_k);
289
290 //top sphere
291 //CircleDomain cr_tf(vec_ts,vec_i,1.0), cr_tn(vec_ts,vec_k,1.0);
292 //bottom sphere
293 //CircleDomain cr_bf(vec_bs,vec_i,1.0), cr_bn(vec_bs,vec_k,1.0);
294 static CircleDomain cr_ct(vec_k,vec_k,1.0);
295
296 //++ 2D domains ++
297
298 //cube
299 static PlanarDomain sf_i(vec_i,vec_i), sf_ni(vec_ni,vec_ni);
300 static PlanarDomain sf_j(vec_j,vec_j), sf_nj(vec_nj,vec_nj);
301 static PlanarDomain sf_k(vec_k,vec_k), sf_nk(vec_nk,vec_nk);
302 //cut
303 //CylinderDomain cy(1.0,vec_k,vec_bs);
304 static SphericalDomain sp_t(vec_k,1.0);//, sp_b(vec_bs,1.0);
305
306 MeshDomain* base_domains[] = {
307 &pt_tfl, &pt_tfr,
308 &pt_tbl, &pt_tbr,
309 &pt_bfl, &pt_bfr,
310 &pt_bbl, &pt_bbr,
311
312 &ln_tf, &lin_tb,
313 &ln_tl, &lin_tr,
314 &ln_bf, &lin_bb,
315 &ln_bl, &lin_br,
316 &ln_lf, &lin_rf,
317 &ln_lb, &lin_rb,
318
319 // &cr_tf, &cr_tn,
320 // &cr_bf, &cr_bn,
321 &cr_ct,
322
323 &sf_i, &sf_ni,
324 &sf_j, &sf_nj,
325 &sf_k, &sf_nk,
326 //&cy,
327 &sp_t//, &sp_b
328 };
329 const int NDOM = sizeof(base_domains)/sizeof(base_domains[0]);
330
331 int dim_array[NDOM] = {
332 0,0,0,0,0,0,0,0,
333 1,1,1,1,1,1,1,1,
334 1,1,1,1,1,//1,1,1,
335 2,2,2,2,2,2,2//,2,2
336 };
337
338 //++ MESH & DOMAIN ++
339
340 mesh.read_vtk(filename.c_str(), err); MSQ_ERRZERO(err);
341 DomainClassifier::classify_skin_geometrically (geom, &mesh, 0.1, base_domains, dim_array, NDOM, err);
342 MSQ_ERRZERO(err);
343 mesh.set_skin_flags( false, false, true, err ); MSQ_ERRZERO(err);
344
345 return filename;
346 }
347
348
get_cut_cube_example(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)349 std::string get_cut_cube_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
350 {
351 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/cutCube.vtk";
352
353 const Vector3D vec_i(5,0,0), vec_ni(-5,0,0);
354 const Vector3D vec_j(0,5,0), vec_nj(0,-5,0);
355 const Vector3D vec_k(0,0,5), vec_nk(0,0,-5);
356
357 const Vector3D vec_tfl(5,-5,5), vec_tfr(5,5,5);
358 const Vector3D vec_tbl(-5,-5,5), vec_tbr(-5,5,5);
359 const Vector3D vec_bfl(5,-5,-5), vec_bfr(5,5,-5);
360 const Vector3D vec_bbl(-5,-5,-5), vec_bbr(-5,5,-5);
361
362 const Vector3D vec_ts(5,0,2), vec_bs(5,0,-2);
363
364
365 //++ 0D domains ++
366
367 static PointDomain pt_tfl(vec_tfl), pt_tfr(vec_tfr);
368 static PointDomain pt_tbl(vec_tbl), pt_tbr(vec_tbr);
369 static PointDomain pt_bfl(vec_bfl), pt_bfr(vec_bfr);
370 static PointDomain pt_bbl(vec_bbl), pt_bbr(vec_bbr);
371
372
373 //++ 1D domains ++
374
375 //top square
376 static LineDomain ln_tf(vec_tfl,vec_j), lin_tb(vec_tbr,vec_nj);
377 static LineDomain ln_tl(vec_tfl,vec_ni), lin_tr(vec_tbr,vec_i);
378 //bottom square
379 static LineDomain ln_bf(vec_bfl,vec_j), lin_bb(vec_bbr,vec_nj);
380 static LineDomain ln_bl(vec_bfl,vec_ni), lin_br(vec_bbr,vec_i);
381 //sides
382 static LineDomain ln_lf(vec_bfl,vec_k), lin_rf(vec_bfr,vec_k);
383 static LineDomain ln_lb(vec_bbl,vec_k), lin_rb(vec_bbr,vec_k);
384
385 //top sphere
386 static CircleDomain cr_tf(vec_ts,vec_i,1.0), cr_tn(vec_ts,vec_k,1.0);
387 //bottom sphere
388 static CircleDomain cr_bf(vec_bs,vec_i,1.0), cr_bn(vec_bs,vec_k,1.0);
389
390
391 //++ 2D domains ++
392
393 //cube
394 static PlanarDomain sf_i(vec_i,vec_i), sf_ni(vec_ni,vec_ni);
395 static PlanarDomain sf_j(vec_j,vec_j), sf_nj(vec_nj,vec_nj);
396 static PlanarDomain sf_k(vec_k,vec_k), sf_nk(vec_nk,vec_nk);
397 //cut
398 static CylinderDomain cy(1.0,vec_k,vec_bs);
399 static SphericalDomain sp_t(vec_ts,1.0), sp_b(vec_bs,1.0);
400
401
402 MeshDomain* base_domains[] = {
403 &pt_tfl, &pt_tfr,
404 &pt_tbl, &pt_tbr,
405 &pt_bfl, &pt_bfr,
406 &pt_bbl, &pt_bbr,
407
408 &ln_tf, &lin_tb,
409 &ln_tl, &lin_tr,
410 &ln_bf, &lin_bb,
411 &ln_bl, &lin_br,
412 &ln_lf, &lin_rf,
413 &ln_lb, &lin_rb,
414 &cr_tf, &cr_tn,
415 &cr_bf, &cr_bn,
416
417 &sf_i, &sf_ni,
418 &sf_j, &sf_nj,
419 &sf_k, &sf_nk,
420 &cy,
421 &sp_t, &sp_b
422 };
423 const int NDOM = sizeof(base_domains)/sizeof(base_domains[0]);
424
425 int dim_array[NDOM] = {
426 0,0,0,0,0,0,0,0,
427 1,1,1,1,1,1,1,1,
428 1,1,1,1,1,1,1,1,
429 2,2,2,2,2,2,2,2,2
430 };
431
432 //++ MESH & DOMAIN ++
433
434 mesh.read_vtk(filename.c_str(), err); MSQ_ERRZERO(err);
435 DomainClassifier::classify_skin_geometrically (geom, &mesh, 0.1, base_domains, dim_array, NDOM, err);
436 MSQ_ERRZERO(err);
437 mesh.set_skin_flags( false, false, true, err ); MSQ_ERRZERO(err);
438
439 return filename;
440 }
441
442
get_sphere_cylinder_example(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)443 std::string get_sphere_cylinder_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
444 {
445 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/sphereCylinder_1194_inv.vtk";
446
447 const Vector3D vec_k(0,0,8), vec_nk(0,0,-8);
448 const Vector3D vec_c(0,0,5);
449
450 //++ 0D domains ++
451
452 //++ 1D domains ++
453
454 //top circle
455 static CircleDomain cr_to(vec_k,vec_k,8.0), cr_ti(vec_k,vec_k,4.0);
456 //bottom circle
457 static CircleDomain cr_bo(vec_nk,vec_nk,8.0);
458
459 //++ 2D domains ++
460
461 static PlanarDomain sf_t(vec_k,vec_k), sf_b(vec_nk,vec_nk);
462 static CylinderDomain cy(8.0,vec_k,vec_k);
463 static SphericalDomain sp(vec_c,5.0);
464
465 MeshDomain* base_domains[] = {
466 &cr_to, &cr_ti, &cr_bo,
467 &sf_t, &sf_b, &cy, &sp
468 };
469 const int NDOM = sizeof(base_domains)/sizeof(base_domains[0]);
470
471 int dim_array[NDOM] = {
472 1, 1, 1,
473 2, 2, 2, 2
474 };
475
476 //++ MESH & DOMAIN ++
477
478 mesh.read_vtk(filename.c_str(), err); MSQ_ERRZERO(err);
479 DomainClassifier::classify_skin_geometrically (geom, &mesh, 0.001, base_domains, dim_array, NDOM, err);
480 MSQ_ERRZERO(err);
481 mesh.set_skin_flags( false, false, true, err ); MSQ_ERRZERO(err);
482
483 return filename;
484 }
485
486
get_hex_3d_part_example(DomainClassifier & geom,MeshImpl & mesh,MsqError & err)487 std::string get_hex_3d_part_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
488 {
489 std::string filename = std::string ( STRINGIFY(SRCDIR) ) + "/hex3Dpart.vtk";
490
491 //2D domains
492
493 const Vector3D vec_i(1,0,0), vec_j(0,1,0), vec_k(0,0,1);
494 const Vector3D vec_ni(-1,0,0), vec_nj(0,-1,0), vec_nk(0,0,-1);
495
496 const Vector3D vec_left(-11.5,0,0), vec_right(11.5,0,0);
497 const Vector3D vec_top(0,5,0), vec_bottom(0,-5,0);
498 const Vector3D vec_front(0,0,5), vec_back(0,0,-5);
499
500 //1D domains
501
502 const Vector3D pt_top_front(0,5,5),pt_top_back(0,5,-5);
503 const Vector3D pt_bottom_front(0,-5,5),pt_bottom_back(0,-5,-5);
504
505 const Vector3D pt_top_left(-11.5,5,0),pt_top_right(11.5,5,0);
506 const Vector3D pt_bottom_left(-11.5,-5,0),pt_bottom_right(11.5,-5,0);
507
508 const Vector3D pt_left_front(-11.5,0,5),pt_left_back(-11.5,0,-5);
509 const Vector3D pt_right_front(11.5,0,5),pt_right_back(11.5,0,-5);
510
511
512 const Vector3D cpt_top_left(-1.5,5,0),cpt_top_right(1.5,5,0);
513 const Vector3D cpt_bottom_left(-1.5,-5,0),cpt_bottom_right(1.5,-5,0);
514
515 //0D domains
516
517 const Vector3D pt_tlf(-11.5,5,5), pt_tlb(-11.5,5,-5);
518 const Vector3D pt_trf(11.5,5,5), pt_trb(11.5,5,-5);
519 const Vector3D pt_blf(-11.5,-5,5), pt_blb(-11.5,-5,-5);
520 const Vector3D pt_brf(11.5,-5,5), pt_brb(11.5,-5,-5);
521
522 const Vector3D pt_c_tlf(-1.5,5,5), pt_c_tlb(-1.5,5,-5);
523 const Vector3D pt_c_trf(1.5,5,5), pt_c_trb(1.5,5,-5);
524 const Vector3D pt_c_blf(-1.5,-5,5), pt_c_blb(-1.5,-5,-5);
525 const Vector3D pt_c_brf(1.5,-5,5), pt_c_brb(1.5,-5,-5);
526
527 static PointDomain p00(pt_tlf);
528 static PointDomain p01(pt_tlb);
529 static PointDomain p02(pt_trf);
530 static PointDomain p03(pt_trb);
531 static PointDomain p04(pt_blf);
532 static PointDomain p05(pt_blb);
533 static PointDomain p06(pt_brf);
534 static PointDomain p07(pt_brb);
535
536 static PointDomain p08(pt_c_tlf);
537 static PointDomain p09(pt_c_tlb);
538 static PointDomain p10(pt_c_trf);
539 static PointDomain p11(pt_c_trb);
540 static PointDomain p12(pt_c_blf);
541 static PointDomain p13(pt_c_blb);
542 static PointDomain p14(pt_c_brf);
543 static PointDomain p15(pt_c_brb);
544
545 static CircleDomain c00(pt_top_front,vec_k,1.5);
546 static CircleDomain c01(pt_top_back,vec_k,1.5);
547 static CircleDomain c02(pt_bottom_front,vec_k,1.5);
548 static CircleDomain c03(pt_bottom_back,vec_k,1.5);
549
550 static LineDomain l00(cpt_top_left,vec_k);
551 static LineDomain l01(cpt_top_right,vec_k);
552 static LineDomain l02(cpt_bottom_left,vec_k);
553 static LineDomain l03(cpt_bottom_right,vec_k);
554
555 static LineDomain l04(pt_top_front,vec_i);
556 static LineDomain l05(pt_top_back,vec_i);
557 static LineDomain l06(pt_bottom_front,vec_i);
558 static LineDomain l07(pt_bottom_back,vec_i);
559 static LineDomain l08(pt_top_left,vec_k);
560 static LineDomain l09(pt_top_right,vec_k);
561 static LineDomain l10(pt_bottom_left,vec_k);
562 static LineDomain l11(pt_bottom_right,vec_k);
563 static LineDomain l12(pt_left_front,vec_j);
564 static LineDomain l13(pt_left_back,vec_j);
565 static LineDomain l14(pt_right_front,vec_j);
566 static LineDomain l15(pt_right_back,vec_j);
567
568 static CylinderDomain C00(1.5,vec_k,vec_top,false);
569 static CylinderDomain C01(1.5,vec_k,vec_bottom,false);
570 static PlanarDomain P00(vec_ni,vec_left);
571 static PlanarDomain P01(vec_i,vec_right);
572 static PlanarDomain P02(vec_j,vec_top);
573 static PlanarDomain P03(vec_nj,vec_bottom);
574 static PlanarDomain P04(vec_k,vec_front);
575 static PlanarDomain P05(vec_nk,vec_back);
576
577 const int NDOM = 44;
578 MeshDomain* base_domains[NDOM] = {
579 &p00,
580 &p01,
581 &p02,
582 &p03,
583 &p04,
584 &p05,
585 &p06,
586 &p07,
587
588 &p08,
589 &p09,
590 &p10,
591 &p11,
592 &p12,
593 &p13,
594 &p14,
595 &p15,
596
597 &c00,
598 &c01,
599 &c02,
600 &c03,
601
602 &l00,
603 &l01,
604 &l02,
605 &l03,
606
607 &l04,
608 &l05,
609 &l06,
610 &l07,
611 &l08,
612 &l09,
613 &l10,
614 &l11,
615 &l12,
616 &l13,
617 &l14,
618 &l15,
619
620 &C00,
621 &C01,
622 &P00,
623 &P01,
624 &P02,
625 &P03,
626 &P04,
627 &P05
628 };
629
630 int dim_array[NDOM] = {
631 0, 0, 0, 0, 0, 0, 0, 0,
632 0, 0, 0, 0, 0, 0, 0, 0,
633 1, 1, 1, 1,
634 1, 1, 1, 1,
635 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
636 2, 2, 2, 2, 2, 2, 2, 2
637 };
638
639
640 //++ MESH & DOMAIN ++
641
642 mesh.read_vtk(filename.c_str(), err); MSQ_ERRZERO(err);
643 DomainClassifier::classify_skin_geometrically (geom, &mesh, 0.1, base_domains, dim_array, NDOM, err);
644 MSQ_ERRZERO(err);
645 mesh.set_skin_flags( false, false, true, err ); MSQ_ERRZERO(err);
646
647 return filename;
648 }
649
650
651
652