9 #include <boost/lexical_cast.hpp>
17 cerr <<
"Usage "<<binary_name<<
" meshtype splits nx ny sx sy nummodes output_file.xml\n";
18 cerr <<
"where 'meshtype' is\n";
19 cerr <<
" = 0 for regular diamond mesh\n";
20 cerr <<
" = 1 for regular diamond mesh with even vertical stripes of diamonds not being split\n";
21 cerr <<
" = 2 for regular diamond mesh with even horizontal stripes of diamonds shifted one quad right\n";
22 cerr <<
" = 3 for a single quad split into triangles symmetrically relative to its central point\n";
23 cerr <<
" 'splits' is:\n";
24 cerr <<
" - for 'meshtype' = 0 - a number of points that splits every edge on each even horizontal line,\n";
25 cerr <<
" - for 'meshtype' = 1 - a number of points that splits every odd edge on each even horizontal line,\n";
26 cerr <<
" - for 'meshtype' = 2 - a number of points that splits every interior edge on each even horizontal line,\n";
27 cerr <<
" - for 'meshtype' = 3 - a number of points splitting every symmetrical part of boundary,\n";
28 cerr <<
" nx is the number of points in x direction that forms quadrilateral grid skeleton,\n";
29 cerr <<
" ny is the number of points in y direction that forms quadrilateral grid skeleton,\n";
30 cerr <<
" sx is the coordinate scaling in x direction,\n";
31 cerr <<
" sy is the coordinate scaling in y direction,\n";
32 cerr <<
" nummodes is the number of boundary modes.\n";
33 cerr <<
"For meshtype in {0,1,2} it generates regular mesh with triangles filling quadrilateral grid (aka skeleton)\n";
34 cerr <<
"in the way that forms vertices of different valence.\n";
35 cerr <<
"All vertex coordinates of quadrilateral grid skeleton are evenly distributed within\n";
36 cerr <<
"unit square but then scaled by means of XSCALE and YSCALE attributes of VERTEX section.\n";
53 void init(
int id,
double x,
double y)
103 vector<Segment> south, north, east,
west;
120 assert(xc.size() > 2);
121 assert(yc.size() > 2);
123 assert(xc.size() % 2 == 1);
124 assert(yc.size() % 2 == 1);
126 double x_split_inc = (xc[1]-xc[0])/(splits+1);
137 vector<vector<Vertex> > verts(yc.size());
138 for (j = 0; j < yc.size(); j++)
140 verts[j].resize(xc.size());
141 for (i = 0; i < xc.size(); i++)
143 verts[j][i].init (vertex_id++, xc[i], yc[j]);
144 mesh.
vert.push_back( verts[j][i] );
149 vector<vector<Segment> > hseg(yc.size());
150 vector<vector<Segment> > vseg(yc.size()-1);
152 for (j = 0; j < verts.size(); j++)
154 hseg[j].resize(xc.size()-1);
158 for (j = 0; j < verts.size(); j+=2)
160 for (i = 0; i < verts[j].size()-1; i++)
162 hseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+0][i+1]);
163 mesh.
seg.push_back( hseg[j][i] );
167 mesh.
south.push_back( hseg[j][i] );
169 if (j == verts[j].size()-1)
171 mesh.
north.push_back( hseg[j][i] );
176 for (j = 0; j < verts.size()-1; j++)
178 vseg[j].resize(xc.size());
179 for (i = 0; i < verts[j].size(); i++)
181 vseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+1][i+0]);
182 mesh.
seg.push_back( vseg[j][i] );
186 mesh.
west.push_back( vseg[j][i] );
188 if (i == verts[j].size()-1)
190 mesh.
east.push_back( vseg[j][i] );
197 for (j = 0; j < yc.size()-1; j+=2)
199 for (i = 0; i < xc.size()-1; i+=2)
202 vector<Vertex> v_cli, v_cri;
203 v_cli.push_back( verts[j+1][i+0] );
204 v_cri.push_back( verts[j+1][i+2] );
205 for (k = 0; k < splits; k++)
208 v_cli.push_back(
Vertex (vertex_id++, xc[i+0]+(k+1)*x_split_inc, yc[j+1]) );
209 v_cri.push_back(
Vertex (vertex_id++, xc[i+2]-(k+1)*x_split_inc, yc[j+1]) );
211 v_cli.push_back( verts[j+1][i+1] );
212 v_cri.push_back( verts[j+1][i+1] );
215 mesh.
vert.insert( mesh.
vert.end(), v_cli.begin()+1, v_cli.end()-1 );
216 mesh.
vert.insert( mesh.
vert.end(), v_cri.begin()+1, v_cri.end()-1 );
219 vector<Segment> s_uld, s_urd, s_lld, s_lrd;
221 vector<Segment> s_clh, s_crh;
223 for (k = 0; k < splits+1; k++)
225 s_uld.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cli[k]) );
226 s_lld.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cli[k]) );
227 s_urd.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cri[k]) );
228 s_lrd.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cri[k]) );
232 for (k = 0; k < splits+1; k++)
234 s_clh.push_back(
Segment(segment_id++, v_cli[k], v_cli[k+1]) );
235 s_crh.push_back(
Segment(segment_id++, v_cri[k], v_cri[k+1]) );
239 mesh.
seg.insert( mesh.
seg.end(), s_uld.begin(), s_uld.end() );
240 mesh.
seg.insert( mesh.
seg.end(), s_lld.begin(), s_lld.end() );
241 mesh.
seg.insert( mesh.
seg.end(), s_urd.begin(), s_urd.end() );
242 mesh.
seg.insert( mesh.
seg.end(), s_lrd.begin(), s_lrd.end() );
243 mesh.
seg.insert( mesh.
seg.end(), s_clh.begin(), s_clh.end() );
244 mesh.
seg.insert( mesh.
seg.end(), s_crh.begin(), s_crh.end() );
247 mesh.
tri.push_back(
Triangle( hseg[j+2][i+0], vseg[j+1][i+0], s_uld[0]) );
248 mesh.
tri.push_back(
Triangle( hseg[j+0][i+0], s_lld[0], vseg[j+0][i+0]) );
249 mesh.
tri.push_back(
Triangle( hseg[j+2][i+1], s_urd[0], vseg[j+1][i+2]) );
250 mesh.
tri.push_back(
Triangle( hseg[j+0][i+1], vseg[j+0][i+2], s_lrd[0]) );
253 for (k = 0; k < splits; k++)
255 mesh.
tri.push_back(
Triangle( s_uld[k+0], s_clh[k+0], s_uld[k+1]) );
256 mesh.
tri.push_back(
Triangle( s_lld[k+0], s_lld[k+1], s_clh[k+0]) );
257 mesh.
tri.push_back(
Triangle( s_urd[k+0], s_urd[k+1], s_crh[k+0]) );
258 mesh.
tri.push_back(
Triangle( s_lrd[k+1], s_lrd[k+0], s_crh[k+0]) );
262 mesh.
tri.push_back(
Triangle( s_clh[splits], vseg[j+1][i+1], s_uld[splits]) );
263 mesh.
tri.push_back(
Triangle( s_clh[splits], s_lld[splits], vseg[j+0][i+1]) );
264 mesh.
tri.push_back(
Triangle( s_crh[splits], s_urd[splits], vseg[j+1][i+1]) );
265 mesh.
tri.push_back(
Triangle( s_crh[splits], vseg[j+0][i+1], s_lrd[splits]) );
282 assert(xc.size() > 2);
283 assert(yc.size() > 2);
285 assert(xc.size() % 2 == 1);
286 assert(yc.size() % 2 == 1);
288 double x_split_inc = (xc[1]-xc[0])/(splits+1);
299 vector<vector<Vertex> > verts(yc.size());
300 for (j = 0; j < yc.size(); j++)
302 verts[j].resize(xc.size());
303 for (i = 0; i < xc.size(); i++)
305 verts[j][i].init (vertex_id++, xc[i], yc[j]);
306 mesh.
vert.push_back( verts[j][i] );
311 vector<vector<Segment> > hseg(yc.size());
312 vector<vector<Segment> > vseg(yc.size()-1);
314 for (j = 0; j < verts.size(); j++)
316 hseg[j].resize(xc.size()-1);
320 for (j = 0; j < verts.size(); j+=2)
322 for (i = 0; i < verts[j].size()-1; i++)
324 hseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+0][i+1]);
325 mesh.
seg.push_back( hseg[j][i] );
329 mesh.
south.push_back( hseg[j][i] );
331 if (j == verts[j].size()-1)
333 mesh.
north.push_back( hseg[j][i] );
338 for (j = 0; j < verts.size()-1; j++)
340 vseg[j].resize(xc.size());
341 for (i = 0; i < verts[j].size(); i++)
343 vseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+1][i+0]);
344 mesh.
seg.push_back( vseg[j][i] );
348 mesh.
west.push_back( vseg[j][i] );
350 if (i == verts[j].size()-1)
352 mesh.
east.push_back( vseg[j][i] );
359 for (j = 0; j < yc.size()-1; j+=2)
361 for (i = 0; i < xc.size()-1; i+=2)
366 int this_diamond_splits = splits;
369 this_diamond_splits = 0;
373 vector<Vertex> v_cli, v_cri;
374 v_cli.push_back( verts[j+1][i+0] );
375 v_cri.push_back( verts[j+1][i+2] );
376 for (k = 0; k < this_diamond_splits; k++)
379 v_cli.push_back(
Vertex (vertex_id++, xc[i+0]+(k+1)*x_split_inc, yc[j+1]) );
380 v_cri.push_back(
Vertex (vertex_id++, xc[i+2]-(k+1)*x_split_inc, yc[j+1]) );
382 v_cli.push_back( verts[j+1][i+1] );
383 v_cri.push_back( verts[j+1][i+1] );
386 mesh.
vert.insert( mesh.
vert.end(), v_cli.begin()+1, v_cli.end()-1 );
387 mesh.
vert.insert( mesh.
vert.end(), v_cri.begin()+1, v_cri.end()-1 );
390 vector<Segment> s_uld, s_urd, s_lld, s_lrd;
392 vector<Segment> s_clh, s_crh;
394 for (k = 0; k < this_diamond_splits+1; k++)
396 s_uld.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cli[k]) );
397 s_lld.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cli[k]) );
398 s_urd.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cri[k]) );
399 s_lrd.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cri[k]) );
403 for (k = 0; k < this_diamond_splits+1; k++)
405 s_clh.push_back(
Segment(segment_id++, v_cli[k], v_cli[k+1]) );
406 s_crh.push_back(
Segment(segment_id++, v_cri[k], v_cri[k+1]) );
410 mesh.
seg.insert( mesh.
seg.end(), s_uld.begin(), s_uld.end() );
411 mesh.
seg.insert( mesh.
seg.end(), s_lld.begin(), s_lld.end() );
412 mesh.
seg.insert( mesh.
seg.end(), s_urd.begin(), s_urd.end() );
413 mesh.
seg.insert( mesh.
seg.end(), s_lrd.begin(), s_lrd.end() );
414 mesh.
seg.insert( mesh.
seg.end(), s_clh.begin(), s_clh.end() );
415 mesh.
seg.insert( mesh.
seg.end(), s_crh.begin(), s_crh.end() );
418 mesh.
tri.push_back(
Triangle( hseg[j+2][i+0], vseg[j+1][i+0], s_uld[0]) );
419 mesh.
tri.push_back(
Triangle( hseg[j+0][i+0], s_lld[0], vseg[j+0][i+0]) );
420 mesh.
tri.push_back(
Triangle( hseg[j+2][i+1], s_urd[0], vseg[j+1][i+2]) );
421 mesh.
tri.push_back(
Triangle( hseg[j+0][i+1], vseg[j+0][i+2], s_lrd[0]) );
424 for (k = 0; k < this_diamond_splits; k++)
426 mesh.
tri.push_back(
Triangle( s_uld[k+0], s_clh[k+0], s_uld[k+1]) );
427 mesh.
tri.push_back(
Triangle( s_lld[k+0], s_lld[k+1], s_clh[k+0]) );
428 mesh.
tri.push_back(
Triangle( s_urd[k+0], s_urd[k+1], s_crh[k+0]) );
429 mesh.
tri.push_back(
Triangle( s_lrd[k+1], s_lrd[k+0], s_crh[k+0]) );
433 mesh.
tri.push_back(
Triangle( s_clh[this_diamond_splits], vseg[j+1][i+1], s_uld[this_diamond_splits]) );
434 mesh.
tri.push_back(
Triangle( s_clh[this_diamond_splits], s_lld[this_diamond_splits], vseg[j+0][i+1]) );
435 mesh.
tri.push_back(
Triangle( s_crh[this_diamond_splits], s_urd[this_diamond_splits], vseg[j+1][i+1]) );
436 mesh.
tri.push_back(
Triangle( s_crh[this_diamond_splits], vseg[j+0][i+1], s_lrd[this_diamond_splits]) );
448 assert(xc.size() > 2);
449 assert(yc.size() > 2);
451 assert(xc.size() % 2 == 1);
452 assert(yc.size() % 2 == 1);
454 double x_split_inc = (xc[1]-xc[0])/(splits+1);
465 vector<vector<Vertex> > verts(yc.size());
466 for (j = 0; j < yc.size(); j++)
468 verts[j].resize(xc.size());
469 for (i = 0; i < xc.size(); i++)
471 verts[j][i].init (vertex_id++, xc[i], yc[j]);
472 mesh.
vert.push_back( verts[j][i] );
477 vector<vector<Segment> > hseg(yc.size());
478 vector<vector<Segment> > vseg(yc.size()-1);
480 for (j = 0; j < verts.size(); j++)
482 hseg[j].resize(xc.size()-1);
486 for (j = 0; j < verts.size(); j+=2)
488 for (i = 0; i < verts[j].size()-1; i++)
490 hseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+0][i+1]);
491 mesh.
seg.push_back( hseg[j][i] );
495 mesh.
south.push_back( hseg[j][i] );
497 if (j == verts.size()-1)
499 mesh.
north.push_back( hseg[j][i] );
504 for (j = 0; j < verts.size()-1; j++)
506 vseg[j].resize(xc.size());
507 for (i = 0; i < verts[j].size(); i++)
509 vseg[j][i].init (segment_id++, verts[j+0][i+0], verts[j+1][i+0]);
510 mesh.
seg.push_back( vseg[j][i] );
514 mesh.
west.push_back( vseg[j][i] );
516 if (i == verts[j].size()-1)
518 mesh.
east.push_back( vseg[j][i] );
525 for (j = 0; j < yc.size()-1; j+=2)
527 for (i = ((j % 4) >= 2); i < xc.size()-2; i+=2)
530 vector<Vertex> v_cli, v_cri;
531 v_cli.push_back( verts[j+1][i+0] );
532 v_cri.push_back( verts[j+1][i+2] );
533 for (k = 0; k < splits; k++)
536 v_cli.push_back(
Vertex (vertex_id++, xc[i+0]+(k+1)*x_split_inc, yc[j+1]) );
537 v_cri.push_back(
Vertex (vertex_id++, xc[i+2]-(k+1)*x_split_inc, yc[j+1]) );
539 v_cli.push_back( verts[j+1][i+1] );
540 v_cri.push_back( verts[j+1][i+1] );
543 mesh.
vert.insert( mesh.
vert.end(), v_cli.begin()+1, v_cli.end()-1 );
544 mesh.
vert.insert( mesh.
vert.end(), v_cri.begin()+1, v_cri.end()-1 );
547 vector<Segment> s_uld, s_urd, s_lld, s_lrd;
549 vector<Segment> s_clh, s_crh;
551 for (k = 0; k < splits+1; k++)
553 s_uld.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cli[k]) );
554 s_lld.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cli[k]) );
555 s_urd.push_back(
Segment(segment_id++, verts[j+2][i+1], v_cri[k]) );
556 s_lrd.push_back(
Segment(segment_id++, verts[j+0][i+1], v_cri[k]) );
560 for (k = 0; k < splits+1; k++)
562 s_clh.push_back(
Segment(segment_id++, v_cli[k], v_cli[k+1]) );
563 s_crh.push_back(
Segment(segment_id++, v_cri[k], v_cri[k+1]) );
567 mesh.
seg.insert( mesh.
seg.end(), s_uld.begin(), s_uld.end() );
568 mesh.
seg.insert( mesh.
seg.end(), s_lld.begin(), s_lld.end() );
569 mesh.
seg.insert( mesh.
seg.end(), s_urd.begin(), s_urd.end() );
570 mesh.
seg.insert( mesh.
seg.end(), s_lrd.begin(), s_lrd.end() );
571 mesh.
seg.insert( mesh.
seg.end(), s_clh.begin(), s_clh.end() );
572 mesh.
seg.insert( mesh.
seg.end(), s_crh.begin(), s_crh.end() );
575 mesh.
tri.push_back(
Triangle( hseg[j+2][i+0], vseg[j+1][i+0], s_uld[0]) );
576 mesh.
tri.push_back(
Triangle( hseg[j+0][i+0], s_lld[0], vseg[j+0][i+0]) );
577 mesh.
tri.push_back(
Triangle( hseg[j+2][i+1], s_urd[0], vseg[j+1][i+2]) );
578 mesh.
tri.push_back(
Triangle( hseg[j+0][i+1], vseg[j+0][i+2], s_lrd[0]) );
581 for (k = 0; k < splits; k++)
583 mesh.
tri.push_back(
Triangle( s_uld[k+0], s_clh[k+0], s_uld[k+1]) );
584 mesh.
tri.push_back(
Triangle( s_lld[k+0], s_lld[k+1], s_clh[k+0]) );
585 mesh.
tri.push_back(
Triangle( s_urd[k+0], s_urd[k+1], s_crh[k+0]) );
586 mesh.
tri.push_back(
Triangle( s_lrd[k+1], s_lrd[k+0], s_crh[k+0]) );
590 mesh.
tri.push_back(
Triangle( s_clh[splits], vseg[j+1][i+1], s_uld[splits]) );
591 mesh.
tri.push_back(
Triangle( s_clh[splits], s_lld[splits], vseg[j+0][i+1]) );
592 mesh.
tri.push_back(
Triangle( s_crh[splits], s_urd[splits], vseg[j+1][i+1]) );
593 mesh.
tri.push_back(
Triangle( s_crh[splits], vseg[j+0][i+1], s_lrd[splits]) );
602 int imax = xc.size()-3;
605 Segment s_clh (segment_id++, verts[j+1][imin+0], verts[j+1][imin+1]);
606 Segment s_crh (segment_id++, verts[j+1][imax+2], verts[j+1][imax+1]);
609 Segment s_uld(segment_id++, verts[j+0][imin+0], verts[j+1][imin+1] );
610 Segment s_urd(segment_id++, verts[j+2][imin+0], verts[j+1][imin+1] );
611 Segment s_lld(segment_id++, verts[j+0][imax+2], verts[j+1][imax+1] );
612 Segment s_lrd(segment_id++, verts[j+2][imax+2], verts[j+1][imax+1] );
616 mesh.
seg.push_back( s_clh );
617 mesh.
seg.push_back( s_crh );
618 mesh.
seg.push_back( s_uld );
619 mesh.
seg.push_back( s_urd );
620 mesh.
seg.push_back( s_lld );
621 mesh.
seg.push_back( s_lrd );
624 mesh.
tri.push_back(
Triangle( hseg[j+0][imin+0], vseg[j+0][imin+1], s_uld) );
625 mesh.
tri.push_back(
Triangle( s_clh, vseg[j+0][imin+0], s_uld) );
627 mesh.
tri.push_back(
Triangle( s_clh, s_urd, vseg[j+1][imin+0]) );
628 mesh.
tri.push_back(
Triangle( hseg[j+2][imin+0], s_urd, vseg[j+1][imin+1]) );
630 mesh.
tri.push_back(
Triangle( hseg[j+0][imax+1], s_lld, vseg[j+0][imax+1]) );
631 mesh.
tri.push_back(
Triangle( s_crh, s_lld, vseg[j+0][imax+2]) );
633 mesh.
tri.push_back(
Triangle( s_crh, vseg[j+1][imax+2], s_lrd) );
634 mesh.
tri.push_back(
Triangle( hseg[j+2][imax+1], vseg[j+1][imax+1], s_lrd) );
647 assert(split_param > 0);
653 splits = 2*(split_param-2)+1;
664 vector<Vertex> north_verts(splits+2);
665 vector<Vertex> south_verts(splits+2);
666 vector<Vertex> west_verts(splits+2);
667 vector<Vertex> east_verts(splits+2);
669 Vertex ll (vertex_id++, 0.0, 0.0);
670 Vertex ul (vertex_id++, 0.0, 1.0);
671 Vertex lr (vertex_id++, 1.0, 0.0);
672 Vertex ur (vertex_id++, 1.0, 1.0);
673 Vertex cc( vertex_id++, 0.5, 0.5);
675 mesh.
vert.push_back( ll );
676 mesh.
vert.push_back( ul );
677 mesh.
vert.push_back( ur );
678 mesh.
vert.push_back( lr );
679 mesh.
vert.push_back( cc );
686 north_verts[splits+1] = ur;
687 south_verts[splits+1] = lr;
688 west_verts[splits+1] = ul;
689 east_verts[splits+1] = ur;
691 for (j = 1; j <= splits; j++)
693 north_verts[j].init (vertex_id++, (
double)j * (1.0 / (splits+1)), 1.0);
694 south_verts[j].init (vertex_id++, (
double)j * (1.0 / (splits+1)), 0.0);
695 west_verts[j].init (vertex_id++, 0.0, (
double)j * (1.0 / (splits+1)));
696 east_verts[j].init (vertex_id++, 1.0, (
double)j * (1.0 / (splits+1)));
698 mesh.
vert.push_back( north_verts[j] );
699 mesh.
vert.push_back( south_verts[j] );
700 mesh.
vert.push_back( west_verts[j] );
701 mesh.
vert.push_back( east_verts[j] );
705 vector<Segment> north_seg(splits+1);
706 vector<Segment> south_seg(splits+1);
707 vector<Segment> west_seg(splits+1);
708 vector<Segment> east_seg(splits+1);
710 vector<Segment> north_diag_seg(splits+2);
711 vector<Segment> south_diag_seg(splits+2);
712 vector<Segment> west_diag_seg(splits+2);
713 vector<Segment> east_diag_seg(splits+2);
715 Segment s_ul (segment_id++, ul, cc);
716 Segment s_ll (segment_id++, ll, cc);
717 Segment s_lr (segment_id++, lr, cc);
718 Segment s_ur (segment_id++, ur, cc);
720 mesh.
seg.push_back( s_ul );
721 mesh.
seg.push_back( s_ll );
722 mesh.
seg.push_back( s_lr );
723 mesh.
seg.push_back( s_ur );
725 north_diag_seg[0] = s_ul;
726 south_diag_seg[0] = s_ll;
727 west_diag_seg[0] = s_ll;
728 east_diag_seg[0] = s_lr;
729 north_diag_seg[splits+1] = s_ur;
730 south_diag_seg[splits+1] = s_lr;
731 west_diag_seg[splits+1] = s_ul;
732 east_diag_seg[splits+1] = s_ur;
734 north_seg[0].init (segment_id++, north_verts[0], north_verts[1]);
735 south_seg[0].init (segment_id++, south_verts[0], south_verts[1]);
736 west_seg[0].init (segment_id++, west_verts[0], west_verts[1]);
737 east_seg[0].init (segment_id++, east_verts[0], east_verts[1]);
739 mesh.
seg.push_back( north_seg[0] );
740 mesh.
seg.push_back( south_seg[0] );
741 mesh.
seg.push_back( west_seg[0] );
742 mesh.
seg.push_back( east_seg[0] );
744 for (j = 1; j <= splits; j++)
746 north_seg[j].init (segment_id++, north_verts[j], north_verts[j+1]);
747 south_seg[j].init (segment_id++, south_verts[j], south_verts[j+1]);
748 west_seg[j].init (segment_id++, west_verts[j], west_verts[j+1]);
749 east_seg[j].init (segment_id++, east_verts[j], east_verts[j+1]);
751 north_diag_seg[j].init (segment_id++, north_verts[j], cc);
752 south_diag_seg[j].init (segment_id++, south_verts[j], cc);
753 west_diag_seg[j].init (segment_id++, west_verts[j], cc);
754 east_diag_seg[j].init (segment_id++, east_verts[j], cc);
756 mesh.
seg.push_back( north_seg[j] );
757 mesh.
seg.push_back( south_seg[j] );
758 mesh.
seg.push_back( west_seg[j] );
759 mesh.
seg.push_back( east_seg[j] );
761 mesh.
seg.push_back( north_diag_seg[j] );
762 mesh.
seg.push_back( south_diag_seg[j] );
763 mesh.
seg.push_back( west_diag_seg[j] );
764 mesh.
seg.push_back( east_diag_seg[j] );
767 mesh.
north.insert( mesh.
north.end(), north_seg.begin(), north_seg.end() );
768 mesh.
south.insert( mesh.
south.end(), south_seg.begin(), south_seg.end() );
769 mesh.
west.insert( mesh.
west.end(), west_seg.begin(), west_seg.end() );
770 mesh.
east.insert( mesh.
east.end(), east_seg.begin(), east_seg.end() );
773 for (j = 0; j < splits+1; j++)
775 mesh.
tri.push_back(
Triangle( north_diag_seg[j], north_diag_seg[j+1], north_seg[j] ) );
776 mesh.
tri.push_back(
Triangle( south_diag_seg[j], south_seg[j], south_diag_seg[j+1] ) );
777 mesh.
tri.push_back(
Triangle( west_diag_seg[j], west_diag_seg[j+1], west_seg[j] ) );
778 mesh.
tri.push_back(
Triangle( east_diag_seg[j], east_seg[j], east_diag_seg[j+1] ) );
784 int main(
int argc,
char *argv[])
786 vector<double> xc, yc;
805 type = boost::lexical_cast<
int>(argv[1]);
806 splits = boost::lexical_cast<
int>(argv[2]);
807 nx = boost::lexical_cast<
int>(argv[3]);
808 ny = boost::lexical_cast<
int>(argv[4]);
809 sx = boost::lexical_cast<
double>(argv[5]);
810 sy = boost::lexical_cast<
double>(argv[6]);
811 nummodes = boost::lexical_cast<
int>(argv[7]);
812 output_file = boost::lexical_cast<
string>(argv[8]);
813 output.open(output_file.c_str());
815 for (i = 0; i < nx; i++)
817 xc.push_back( (
double)i * (1.0 / (nx-1)) );
819 for (i = 0; i < ny; i++)
821 yc.push_back( (
double)i * (1.0 / (ny-1)) );
840 cout <<
"strange mesh type, aborting" << endl;
844 output<<
"<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
845 output<<
"<NEKTAR>" << endl;
846 output <<
"<EXPANSIONS>" << endl;
847 output <<
"<E COMPOSITE=\"C[0]\" NUMMODES=\"" << nummodes <<
"\" FIELDS=\"u\" TYPE=\"MODIFIED\" />" <<endl;
848 output <<
"</EXPANSIONS>\n" << endl;
852 output <<
"<GEOMETRY DIM=\"2\" SPACE=\"2\">" << endl;
858 output <<
"<VERTEX XSCALE=\"" << sx <<
"\" YSCALE=\"" << sy <<
"\">" << endl;
860 for (i = 0; i < mesh.
vert.size(); i++)
862 output <<
" <V ID=\"";
863 output << mesh.
vert[i].m_id <<
"\">\t";
864 output << mesh.
vert[i].m_x <<
" ";
865 output << mesh.
vert[i].m_y <<
" 0.0 </V>" << endl;
868 output <<
" </VERTEX>\n" << endl;
874 output <<
" <EDGE>" << endl;
875 for(i = 0; i < mesh.
seg.size(); i++)
877 output <<
" <E ID=\"";
878 output << mesh.
seg[i].m_id <<
"\">\t";
879 output << mesh.
seg[i].m_v1.m_id <<
" ";
880 output << mesh.
seg[i].m_v2.m_id <<
" </E>" << endl;
882 output <<
" </EDGE>\n" << endl;
888 output <<
" <ELEMENT>" << endl;
889 for(i = 0; i < mesh.
tri.size(); i++)
891 output <<
" <T ID=\"";
892 output << i <<
"\">\t";
893 output << mesh.
tri[i].m_s1.m_id <<
" ";
894 output << mesh.
tri[i].m_s2.m_id <<
" ";
895 output << mesh.
tri[i].m_s3.m_id <<
" </T>" << endl;
897 output <<
" </ELEMENT>\n" << endl;
903 output <<
"<COMPOSITE>" << endl;
904 output <<
"<C ID=\"0\"> T[0-" << mesh.
tri.size()-1 <<
"] </C>" << endl;
907 output <<
"<C ID=\"1\"> E[";
908 for(i = 0; i < mesh.
south.size(); ++i)
910 output << mesh.
south[i].m_id;
911 if(i != mesh.
south.size()-1)
916 output <<
"] </C> // south border" << endl;
918 output <<
"<C ID=\"2\"> E[";
920 for(i = 0; i < mesh.
west.size(); ++i)
922 output << mesh.
west[i].m_id;
923 if(i != mesh.
west.size()-1)
928 output <<
"] </C> // west border" << endl;
930 output <<
"<C ID=\"3\"> E[";
931 for(i = 0; i < mesh.
north.size(); ++i)
933 output << mesh.
north[i].m_id;
934 if(i != mesh.
north.size()-1)
939 output <<
"] </C> // north border" << endl;
941 output <<
"<C ID=\"4\"> E[";
942 for(i = 0; i < mesh.
east.size(); ++i)
944 output << mesh.
east[i].m_id;
945 if(i != mesh.
east.size()-1)
950 output <<
"] </C> // East border" << endl;
952 output <<
"</COMPOSITE>\n" << endl;
955 output <<
"<DOMAIN> C[0] </DOMAIN>\n" << endl;
956 output <<
"</GEOMETRY>\n" << endl;
957 output <<
"</NEKTAR>" << endl;
962 cerr <<
"Something went wrong. Caught an exception, stop." << endl;
974 output <<
"<CONDITIONS>" << endl;
976 output <<
"<SOLVERINFO>" << endl;
977 output <<
"<I PROPERTY=\"GlobalSysSoln\" VALUE=\"DirectFull\"/>" << endl;
978 output <<
"</SOLVERINFO>\n" << endl;
980 output <<
"<PARAMETERS>" << endl;
981 output <<
"<P> TimeStep = 0.002 </P>" << endl;
982 output <<
"<P> Lambda = 1 </P>" << endl;
983 output <<
"</PARAMETERS>\n" << endl;
985 output <<
"<VARIABLES>" << endl;
986 output <<
" <V ID=\"0\"> u </V>" << endl;
987 output <<
"</VARIABLES>\n" << endl;
989 output <<
"<BOUNDARYREGIONS>" << endl;
990 output <<
"<B ID=\"0\"> C[1] </B>" << endl;
991 output <<
"<B ID=\"1\"> C[2] </B>" << endl;
992 output <<
"<B ID=\"2\"> C[3] </B>" << endl;
993 output <<
"<B ID=\"3\"> C[4] </B>" << endl;
994 output <<
"</BOUNDARYREGIONS>\n" << endl;
996 output <<
"<BOUNDARYCONDITIONS>" << endl;
997 output <<
" <REGION REF=\"0\"> // South border " << endl;
998 output <<
" <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
999 output <<
" </REGION>" << endl;
1001 output <<
" <REGION REF=\"1\"> // West border " << endl;
1002 output <<
" <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
1003 output <<
" </REGION>" << endl;
1005 output <<
" <REGION REF=\"2\"> // North border " << endl;
1006 output <<
" <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
1007 output <<
" </REGION>" << endl;
1009 output <<
" <REGION REF=\"3\"> // East border " << endl;
1010 output <<
" <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
1011 output <<
" </REGION>" << endl;
1012 output <<
"</BOUNDARYCONDITIONS>" << endl;
1014 output <<
" <FUNCTION NAME=\"Forcing\">" << endl;
1015 output <<
" <E VAR=\"u\" VALUE=\"-(Lambda + 2*PI*PI/4)*sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
1016 output <<
" </FUNCTION>" << endl;
1018 output <<
" <FUNCTION NAME=\"ExactSolution\">" << endl;
1019 output <<
" <E VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
1020 output <<
" </FUNCTION>" << endl;
1022 output <<
"</CONDITIONS>" << endl;
Mesh generateSimilarDiamondsMesh(vector< double > &xc, vector< double > &yc, int splits)
Vertex(int id, double x, double y)
Mesh generateDiamondMeshWithHorizontalShifts(vector< double > &xc, vector< double > &yc, int splits)
Mesh generateDiamondMeshDifferentlySplit(vector< double > &xc, vector< double > &yc, int splits)
void init(int id, Vertex v1, Vertex v2)
int main(int argc, char *argv[])
void print_usage_info(char *binary_name)
Triangle(Segment s1, Segment s2, Segment s3)
Represents a vertex in the mesh.
void PrintConditions(ofstream &output)
void init(int id, double x, double y)
Segment(int id, Vertex v1, Vertex v2)
Mesh generateStarcutSingleQuadMesh(int split_param)