Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
VariableValence2DMeshGenerator.cpp
Go to the documentation of this file.
1 #include <cstdlib>
2 #include <string>
3 #include <cassert>
4 
5 #include <iostream>
6 #include <fstream>
7 #include <vector>
8 
9 #include <boost/lexical_cast.hpp>
10 
11 using namespace std;
12 
13 
14 void print_usage_info(char* binary_name)
15 {
16  // argv: 1 2 3 4 5 6 7 8
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";
37 }
38 
39 
40 
41 
42 void PrintConditions(ofstream& output);
43 
44 
45 struct Vertex
46 {
47  Vertex() {}
48  Vertex(int id, double x, double y):
49  m_id(id),
50  m_x(x),
51  m_y(y)
52  {}
53  void init(int id, double x, double y)
54  {
55  m_id = id;
56  m_x = x;
57  m_y = y;
58  }
59 
60  int m_id;
61  double m_x;
62  double m_y;
63 };
64 
65 struct Segment
66 {
67  Segment() {}
68  Segment(int id, Vertex v1, Vertex v2):
69  m_id(id),
70  m_v1(v1),
71  m_v2(v2)
72  {}
73  void init(int id, Vertex v1, Vertex v2)
74  {
75  m_id = id;
76  m_v1 = v1;
77  m_v2 = v2;
78  }
79 
80  int m_id;
81  Vertex m_v1, m_v2;
82 };
83 
84 struct Triangle
85 {
86  Triangle() {}
87  Triangle(Segment s1, Segment s2, Segment s3):
88  m_s1(s1),
89  m_s2(s2),
90  m_s3(s3)
91  {}
92 
93  Segment m_s1, m_s2, m_s3;
94 };
95 
96 struct Mesh
97 {
98  vector<Vertex> vert;
99  vector<Segment> seg;
100  vector<Triangle> tri;
101 
102  // boundary
103  vector<Segment> south, north, east, west;
104 };
105 
106 enum MeshType
107 {
108  eRegularGridOfSimilarDiamonds,
109  eRegularGridOfDiamondsDifferentlySplit,
110  eRegularGridofDiamondsWithHorizontalShifts,
111  eStarcutSingleQuadrilateral,
112  eDummy
113 };
114 
115 // Calculates geometrical info for four adjacent quadrilaterals
116 // subdivided by diagonal segments to form triangular submesh.
117 Mesh generateSimilarDiamondsMesh(vector<double>& xc, vector<double>& yc, int splits)
118 {
119  // at least 3 points
120  assert(xc.size() > 2);
121  assert(yc.size() > 2);
122  // number of points should be even
123  assert(xc.size() % 2 == 1);
124  assert(yc.size() % 2 == 1);
125 
126  double x_split_inc = (xc[1]-xc[0])/(splits+1);
127 
128  Mesh mesh;
129 
130  int vertex_id = 0;
131  int segment_id = 0;
132  int i = 0;
133  int j = 0;
134  int k = 0;
135 
136  // prepare regular grid of vertices
137  vector<vector<Vertex> > verts(yc.size());
138  for (j = 0; j < yc.size(); j++)
139  {
140  verts[j].resize(xc.size());
141  for (i = 0; i < xc.size(); i++)
142  {
143  verts[j][i].init (vertex_id++, xc[i], yc[j]);
144  mesh.vert.push_back( verts[j][i] );
145  }
146  }
147 
148  // prepare carcass of edges
149  vector<vector<Segment> > hseg(yc.size());
150  vector<vector<Segment> > vseg(yc.size()-1);
151 
152  for (j = 0; j < verts.size(); j++)
153  {
154  hseg[j].resize(xc.size()-1);
155  }
156  // prepare only odd horizontal lines.
157  // even lines are split below
158  for (j = 0; j < verts.size(); j+=2)
159  {
160  for (i = 0; i < verts[j].size()-1; i++)
161  {
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] );
164 
165  if (j == 0)
166  {
167  mesh.south.push_back( hseg[j][i] );
168  }
169  if (j == verts[j].size()-1)
170  {
171  mesh.north.push_back( hseg[j][i] );
172  }
173  }
174  }
175 
176  for (j = 0; j < verts.size()-1; j++)
177  {
178  vseg[j].resize(xc.size());
179  for (i = 0; i < verts[j].size(); i++)
180  {
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] );
183 
184  if (i == 0)
185  {
186  mesh.west.push_back( vseg[j][i] );
187  }
188  if (i == verts[j].size()-1)
189  {
190  mesh.east.push_back( vseg[j][i] );
191  }
192  }
193  }
194 
195  // loop through groups of 4 adjacent quadrilaterals
196  // and define internal triangulation of these groups
197  for (j = 0; j < yc.size()-1; j+=2)
198  {
199  for (i = 0; i < xc.size()-1; i+=2)
200  {
201  // vertices splitting even horizontal lines, including border vertices
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++)
206  {
207  // these number from outside of the group towards the center
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]) );
210  }
211  v_cli.push_back( verts[j+1][i+1] );
212  v_cri.push_back( verts[j+1][i+1] );
213 
214  // saving vertices
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 );
217 
218  // define diagonal segments: upper left diagonal, ..
219  vector<Segment> s_uld, s_urd, s_lld, s_lrd;
220  // as well as horizontal segments splitting even lines
221  vector<Segment> s_clh, s_crh;
222 
223  for (k = 0; k < splits+1; k++)
224  {
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]) );
229  }
230  //s_clh.push_back( hseg[j+1][i+0] );
231  //s_crh.push_back( hseg[j+1][i+1] );
232  for (k = 0; k < splits+1; k++)
233  {
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]) );
236  }
237 
238  // saving segments
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() );
245 
246  // corner triangles first
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]) );
251 
252  // internal triangles
253  for (k = 0; k < splits; k++)
254  {
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]) );
259  }
260 
261  // inner central triangles
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]) );
266  }
267  }
268  return mesh;
269 }
270 
271 
272 // Calculates geometrical info for four adjacent quadrilaterals
273 // subdivided by diagonal segments to form triangular submesh.
274 //
275 // Generate regular grid of diamonds where within each group of
276 // 4 quadrilaterals their left quads have 'splits' number of
277 // horizontal splits (same as for regular grid above) while each
278 // right hand side pair of quadrilaterals is not split at all.
279 Mesh generateDiamondMeshDifferentlySplit(vector<double>& xc, vector<double>& yc, int splits)
280 {
281  // at least 3 points
282  assert(xc.size() > 2);
283  assert(yc.size() > 2);
284  // number of points should be even
285  assert(xc.size() % 2 == 1);
286  assert(yc.size() % 2 == 1);
287 
288  double x_split_inc = (xc[1]-xc[0])/(splits+1);
289 
290  Mesh mesh;
291 
292  int vertex_id = 0;
293  int segment_id = 0;
294  int i = 0;
295  int j = 0;
296  int k = 0;
297 
298  // prepare regular grid of vertices
299  vector<vector<Vertex> > verts(yc.size());
300  for (j = 0; j < yc.size(); j++)
301  {
302  verts[j].resize(xc.size());
303  for (i = 0; i < xc.size(); i++)
304  {
305  verts[j][i].init (vertex_id++, xc[i], yc[j]);
306  mesh.vert.push_back( verts[j][i] );
307  }
308  }
309 
310  // prepare carcass of edges
311  vector<vector<Segment> > hseg(yc.size());
312  vector<vector<Segment> > vseg(yc.size()-1);
313 
314  for (j = 0; j < verts.size(); j++)
315  {
316  hseg[j].resize(xc.size()-1);
317  }
318  // prepare only odd horizontal lines.
319  // even lines are split below
320  for (j = 0; j < verts.size(); j+=2)
321  {
322  for (i = 0; i < verts[j].size()-1; i++)
323  {
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] );
326 
327  if (j == 0)
328  {
329  mesh.south.push_back( hseg[j][i] );
330  }
331  if (j == verts[j].size()-1)
332  {
333  mesh.north.push_back( hseg[j][i] );
334  }
335  }
336  }
337 
338  for (j = 0; j < verts.size()-1; j++)
339  {
340  vseg[j].resize(xc.size());
341  for (i = 0; i < verts[j].size(); i++)
342  {
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] );
345 
346  if (i == 0)
347  {
348  mesh.west.push_back( vseg[j][i] );
349  }
350  if (i == verts[j].size()-1)
351  {
352  mesh.east.push_back( vseg[j][i] );
353  }
354  }
355  }
356 
357  // loop through groups of 4 adjacent quadrilaterals
358  // and define internal triangulation of these groups
359  for (j = 0; j < yc.size()-1; j+=2)
360  {
361  for (i = 0; i < xc.size()-1; i+=2)
362  {
363  // make every even diamond in a horizontal direction
364  // be simple diamond with only one diagonal edge splitting
365  // each participating quad
366  int this_diamond_splits = splits;
367  if ((i % 4) >= 2)
368  {
369  this_diamond_splits = 0;
370  }
371 
372  // vertices splitting even horizontal lines, including border vertices
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++)
377  {
378  // these number from outside of the group towards the center
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]) );
381  }
382  v_cli.push_back( verts[j+1][i+1] );
383  v_cri.push_back( verts[j+1][i+1] );
384 
385  // saving vertices
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 );
388 
389  // define diagonal segments: upper left diagonal, ..
390  vector<Segment> s_uld, s_urd, s_lld, s_lrd;
391  // as well as horizontal segments splitting even lines
392  vector<Segment> s_clh, s_crh;
393 
394  for (k = 0; k < this_diamond_splits+1; k++)
395  {
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]) );
400  }
401  //s_clh.push_back( hseg[j+1][i+0] );
402  //s_crh.push_back( hseg[j+1][i+1] );
403  for (k = 0; k < this_diamond_splits+1; k++)
404  {
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]) );
407  }
408 
409  // saving segments
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() );
416 
417  // corner triangles first
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]) );
422 
423  // internal triangles
424  for (k = 0; k < this_diamond_splits; k++)
425  {
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]) );
430  }
431 
432  // inner central triangles
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]) );
437  }
438  }
439  return mesh;
440 }
441 
442 
443 // Calculates geometrical info for four adjacent quadrilaterals
444 // subdivided by diagonal segments to form triangular submesh.
445 Mesh generateDiamondMeshWithHorizontalShifts(vector<double>& xc, vector<double>& yc, int splits)
446 {
447  // at least 3 points
448  assert(xc.size() > 2);
449  assert(yc.size() > 2);
450  // number of points should be even
451  assert(xc.size() % 2 == 1);
452  assert(yc.size() % 2 == 1);
453 
454  double x_split_inc = (xc[1]-xc[0])/(splits+1);
455 
456  Mesh mesh;
457 
458  int vertex_id = 0;
459  int segment_id = 0;
460  int i = 0;
461  int j = 0;
462  int k = 0;
463 
464  // prepare regular grid of vertices
465  vector<vector<Vertex> > verts(yc.size());
466  for (j = 0; j < yc.size(); j++)
467  {
468  verts[j].resize(xc.size());
469  for (i = 0; i < xc.size(); i++)
470  {
471  verts[j][i].init (vertex_id++, xc[i], yc[j]);
472  mesh.vert.push_back( verts[j][i] );
473  }
474  }
475 
476  // prepare carcass of edges
477  vector<vector<Segment> > hseg(yc.size());
478  vector<vector<Segment> > vseg(yc.size()-1);
479 
480  for (j = 0; j < verts.size(); j++)
481  {
482  hseg[j].resize(xc.size()-1);
483  }
484  // prepare only odd horizontal lines.
485  // even lines are split below
486  for (j = 0; j < verts.size(); j+=2)
487  {
488  for (i = 0; i < verts[j].size()-1; i++)
489  {
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] );
492 
493  if (j == 0)
494  {
495  mesh.south.push_back( hseg[j][i] );
496  }
497  if (j == verts.size()-1)
498  {
499  mesh.north.push_back( hseg[j][i] );
500  }
501  }
502  }
503 
504  for (j = 0; j < verts.size()-1; j++)
505  {
506  vseg[j].resize(xc.size());
507  for (i = 0; i < verts[j].size(); i++)
508  {
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] );
511 
512  if (i == 0)
513  {
514  mesh.west.push_back( vseg[j][i] );
515  }
516  if (i == verts[j].size()-1)
517  {
518  mesh.east.push_back( vseg[j][i] );
519  }
520  }
521  }
522 
523  // loop through groups of 4 adjacent quadrilaterals
524  // and define internal triangulation of these groups
525  for (j = 0; j < yc.size()-1; j+=2)
526  {
527  for (i = ((j % 4) >= 2); i < xc.size()-2; i+=2)
528  {
529  // vertices splitting even horizontal lines, including border vertices
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++)
534  {
535  // these number from outside of the group towards the center
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]) );
538  }
539  v_cli.push_back( verts[j+1][i+1] );
540  v_cri.push_back( verts[j+1][i+1] );
541 
542  // saving vertices
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 );
545 
546  // define diagonal segments: upper left diagonal, ..
547  vector<Segment> s_uld, s_urd, s_lld, s_lrd;
548  // as well as horizontal segments splitting even lines
549  vector<Segment> s_clh, s_crh;
550 
551  for (k = 0; k < splits+1; k++)
552  {
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]) );
557  }
558  //s_clh.push_back( hseg[j+1][i+0] );
559  //s_crh.push_back( hseg[j+1][i+1] );
560  for (k = 0; k < splits+1; k++)
561  {
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]) );
564  }
565 
566  // saving segments
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() );
573 
574  // corner triangles first
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]) );
579 
580  // internal triangles
581  for (k = 0; k < splits; k++)
582  {
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]) );
587  }
588 
589  // inner central triangles
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]) );
594  }
595 
596  // simplified triangulation for boundary groups of horizontally
597  // shifted stripes
598  if ((j % 4) >= 2)
599  {
600 
601  int imin = 0;
602  int imax = xc.size()-3;
603 
604  // horizontal segments
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]);
607 
608  // diagonal segments
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] );
613 
614 
615  // saving segments
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 );
622 
623  // triangles
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) );
626 
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]) );
629 
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]) );
632 
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) );
635  }
636  }
637  return mesh;
638 }
639 
640 
641 // Calculates geometrical info for a single quadrilateral
642 // split into triangles in star-like manner
643 // to form triangular mesh.
644 Mesh generateStarcutSingleQuadMesh(int split_param)
645 {
646  // at least 4 triangles
647  assert(split_param > 0);
648 
649  // number of internal splits for each boundary edge of quad
650  int splits = 0;
651  if (split_param > 1)
652  {
653  splits = 2*(split_param-2)+1;
654  }
655 
656  Mesh mesh;
657 
658  int vertex_id = 0;
659  int segment_id = 0;
660  int j = 0;
661 
662 
663  // prepare boundary vertices
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);
668 
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);
674 
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 );
680 
681  // prepare vertices inner to every boundary edge
682  north_verts[0] = ul;
683  south_verts[0] = ll;
684  west_verts[0] = ll;
685  east_verts[0] = lr;
686  north_verts[splits+1] = ur;
687  south_verts[splits+1] = lr;
688  west_verts[splits+1] = ul;
689  east_verts[splits+1] = ur;
690 
691  for (j = 1; j <= splits; j++)
692  {
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)));
697 
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] );
702  }
703 
704  // prepare carcass of edges
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);
709 
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);
714 
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);
719 
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 );
724 
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;
733 
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]);
738 
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] );
743 
744  for (j = 1; j <= splits; j++)
745  {
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]);
750 
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);
755 
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] );
760 
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] );
765  }
766 
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() );
771 
772  // define internal triangulation of a quad
773  for (j = 0; j < splits+1; j++)
774  {
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] ) );
779  }
780  return mesh;
781 }
782 
783 
784 int main(int argc, char *argv[])
785 {
786  vector<double> xc, yc;
787  int nx = 0;
788  int ny = 0;
789  double sx = 1.0;
790  double sy = 1.0;
791  int nummodes = 7;
792  int splits = 0;
793  int i;
794  int type;
795  string output_file;
796  ofstream output;
797 
798  if(argc != 9)
799  {
800  print_usage_info(argv[0]);
801  exit(1);
802  }
803 
804  try{
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());
814 
815  for (i = 0; i < nx; i++)
816  {
817  xc.push_back( (double)i * (1.0 / (nx-1)) );
818  }
819  for (i = 0; i < ny; i++)
820  {
821  yc.push_back( (double)i * (1.0 / (ny-1)) );
822  }
823 
824  Mesh mesh;
825  switch (type)
826  {
827  case eRegularGridOfSimilarDiamonds:
828  mesh = generateSimilarDiamondsMesh(xc, yc, splits);
829  break;
830  case eRegularGridOfDiamondsDifferentlySplit:
831  mesh = generateDiamondMeshDifferentlySplit(xc, yc, splits);
832  break;
833  case eRegularGridofDiamondsWithHorizontalShifts:
834  mesh = generateDiamondMeshWithHorizontalShifts(xc, yc, splits);
835  break;
836  case eStarcutSingleQuadrilateral:
837  mesh = generateStarcutSingleQuadMesh(splits);
838  break;
839  default:
840  cout << "strange mesh type, aborting" << endl;
841  throw 1;
842  }
843 
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;
849 
850  PrintConditions(output);
851 
852  output << "<GEOMETRY DIM=\"2\" SPACE=\"2\">" << endl;
853 
854  // -----------------------------------
855  // Vertices
856  // -----------------------------------
857 
858  output << "<VERTEX XSCALE=\"" << sx << "\" YSCALE=\"" << sy << "\">" << endl;
859 
860  for (i = 0; i < mesh.vert.size(); i++)
861  {
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;
866  }
867 
868  output << " </VERTEX>\n" << endl;
869 
870  // -----------------------------------
871  // Edges
872  // -----------------------------------
873 
874  output << " <EDGE>" << endl;
875  for(i = 0; i < mesh.seg.size(); i++)
876  {
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;
881  }
882  output << " </EDGE>\n" << endl;
883 
884  // -----------------------------------
885  // Elements
886  // -----------------------------------
887 
888  output << " <ELEMENT>" << endl;
889  for(i = 0; i < mesh.tri.size(); i++)
890  {
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;
896  }
897  output << " </ELEMENT>\n" << endl;
898 
899  // -----------------------------------
900  // Composites
901  // -----------------------------------
902 
903  output << "<COMPOSITE>" << endl;
904  output << "<C ID=\"0\"> T[0-" << mesh.tri.size()-1 << "] </C>" << endl;
905 
906  // boundary composites coincide for both mesh element types
907  output << "<C ID=\"1\"> E[";
908  for(i = 0; i < mesh.south.size(); ++i)
909  {
910  output << mesh.south[i].m_id;
911  if(i != mesh.south.size()-1)
912  {
913  output << ",";
914  }
915  }
916  output << "] </C> // south border" << endl;
917 
918  output << "<C ID=\"2\"> E[";
919 
920  for(i = 0; i < mesh.west.size(); ++i)
921  {
922  output << mesh.west[i].m_id;
923  if(i != mesh.west.size()-1)
924  {
925  output << ",";
926  }
927  }
928  output << "] </C> // west border" << endl;
929 
930  output << "<C ID=\"3\"> E[";
931  for(i = 0; i < mesh.north.size(); ++i)
932  {
933  output << mesh.north[i].m_id;
934  if(i != mesh.north.size()-1)
935  {
936  output << ",";
937  }
938  }
939  output << "] </C> // north border" << endl;
940 
941  output << "<C ID=\"4\"> E[";
942  for(i = 0; i < mesh.east.size(); ++i)
943  {
944  output << mesh.east[i].m_id;
945  if(i != mesh.east.size()-1)
946  {
947  output << ",";
948  }
949  }
950  output << "] </C> // East border" << endl;
951 
952  output << "</COMPOSITE>\n" << endl;
953 
954 
955  output << "<DOMAIN> C[0] </DOMAIN>\n" << endl;
956  output << "</GEOMETRY>\n" << endl;
957  output << "</NEKTAR>" << endl;
958 
959  }
960  catch(...)
961  {
962  cerr << "Something went wrong. Caught an exception, stop." << endl;
963  return 1;
964  }
965 
966  return 0;
967 
968 }
969 
970 
971 
972 void PrintConditions(ofstream& output)
973 {
974  output << "<CONDITIONS>" << endl;
975 
976  output << "<SOLVERINFO>" << endl;
977  output << "<I PROPERTY=\"GlobalSysSoln\" VALUE=\"DirectFull\"/>" << endl;
978  output << "</SOLVERINFO>\n" << endl;
979 
980  output << "<PARAMETERS>" << endl;
981  output << "<P> TimeStep = 0.002 </P>" << endl;
982  output << "<P> Lambda = 1 </P>" << endl;
983  output << "</PARAMETERS>\n" << endl;
984 
985  output << "<VARIABLES>" << endl;
986  output << " <V ID=\"0\"> u </V>" << endl;
987  output << "</VARIABLES>\n" << endl;
988 
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;
995 
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;
1000 
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;
1004 
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;
1008 
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;
1013 
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;
1017 
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;
1021 
1022  output << "</CONDITIONS>" << endl;
1023 }
1024 
MeshType
MeshType
Definition: VariableValence2DMeshGenerator.cpp:106
generateSimilarDiamondsMesh
Mesh generateSimilarDiamondsMesh(vector< double > &xc, vector< double > &yc, int splits)
Definition: VariableValence2DMeshGenerator.cpp:117
eStarcutSingleQuadrilateral
Definition: VariableValence2DMeshGenerator.cpp:111
eDummy
Definition: VariableValence2DMeshGenerator.cpp:112
Mesh::tri
vector< Triangle > tri
Definition: VariableValence2DMeshGenerator.cpp:100
Vertex::Vertex
Vertex(int id, double x, double y)
Definition: VariableValence2DMeshGenerator.cpp:48
generateDiamondMeshWithHorizontalShifts
Mesh generateDiamondMeshWithHorizontalShifts(vector< double > &xc, vector< double > &yc, int splits)
Definition: VariableValence2DMeshGenerator.cpp:445
Mesh::north
vector< Segment > north
Definition: VariableValence2DMeshGenerator.cpp:103
std
STL namespace.
generateDiamondMeshDifferentlySplit
Mesh generateDiamondMeshDifferentlySplit(vector< double > &xc, vector< double > &yc, int splits)
Definition: VariableValence2DMeshGenerator.cpp:279
Segment::init
void init(int id, Vertex v1, Vertex v2)
Definition: VariableValence2DMeshGenerator.cpp:73
Mesh::east
vector< Segment > east
Definition: VariableValence2DMeshGenerator.cpp:103
eRegularGridOfDiamondsDifferentlySplit
Definition: VariableValence2DMeshGenerator.cpp:109
main
int main(int argc, char *argv[])
Definition: VariableValence2DMeshGenerator.cpp:784
print_usage_info
void print_usage_info(char *binary_name)
Definition: VariableValence2DMeshGenerator.cpp:14
eRegularGridofDiamondsWithHorizontalShifts
Definition: VariableValence2DMeshGenerator.cpp:110
Triangle::Triangle
Triangle(Segment s1, Segment s2, Segment s3)
Definition: VariableValence2DMeshGenerator.cpp:87
eRegularGridOfSimilarDiamonds
Definition: VariableValence2DMeshGenerator.cpp:108
Vertex
Represents a vertex in the mesh.
Definition: makeEllipsoidNektar.c:22
Mesh::west
vector< Segment > west
Definition: VariableValence2DMeshGenerator.cpp:103
Mesh::vert
vector< Vertex > vert
Definition: VariableValence2DMeshGenerator.cpp:98
Mesh::south
vector< Segment > south
Definition: VariableValence2DMeshGenerator.cpp:103
PrintConditions
void PrintConditions(ofstream &output)
Definition: VariableValence2DMeshGenerator.cpp:972
Vertex::init
void init(int id, double x, double y)
Definition: VariableValence2DMeshGenerator.cpp:53
Segment::Segment
Segment(int id, Vertex v1, Vertex v2)
Definition: VariableValence2DMeshGenerator.cpp:68
generateStarcutSingleQuadMesh
Mesh generateStarcutSingleQuadMesh(int split_param)
Definition: VariableValence2DMeshGenerator.cpp:644
Mesh::seg
vector< Segment > seg
Definition: VariableValence2DMeshGenerator.cpp:99
Generated on Wed May 4 2016 17:54:18 for Nektar++ by   doxygen 1.8.8