Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Regular2DMeshGenerator.cpp
Go to the documentation of this file.
1 #include <cstdlib>
2 #include <string>
3 
4 #include <iostream>
5 #include <vector>
6 #include <map>
7 
8 #include <boost/lexical_cast.hpp>
9 
10 #include<SpatialDomains/Geometry.h>
11 
12 using namespace std;
13 using namespace Nektar;
14 
15 
16 void PrintConditions(ofstream& output);
17 
18 void print_usage_info(char* binary_name)
19 {
20  // argv: 1 2 3 4 5 6 7
21  cerr << "Usage "<<binary_name<<" type nx ny sx sy nummodes output_file.xml\n";
22  cerr << "where type is either 'quads' or 'triangles' with no quotation marks,\n";
23  cerr << " nx is the number of points in x direction,\n";
24  cerr << " ny is the number of points in y direction,\n";
25  cerr << " sx is the coordinate scaling in x direction,\n";
26  cerr << " sy is the coordinate scaling in y direction,\n";
27  cerr << " nummodes is the number of boundary modes.\n";
28  cerr << "It generates regular mesh with triangles or quadrilaterals (but not both)\n";
29  cerr << "All vertex coordinates are evenly distributed within unit square but\n";
30  cerr << "then scaled by means of XSCALE and YSCALE attributes of VERTEX section.\n";
31 }
32 
33 int main(int argc, char *argv[])
34 {
35  vector<double> xc,yc;
36  int nx = 0;
37  int ny = 0;
38  double sx = 1.0;
39  double sy = 1.0;
40  int nummodes = 7;
41  string type_name;
42  int i,j,k;
43  string output_file;
44  ofstream output;
45 
46  if(argc != 8)
47  {
48  print_usage_info(argv[0]);
49  exit(1);
50  }
51 
52  try{
53  type_name = boost::lexical_cast<string>(argv[1]);
54  nx = boost::lexical_cast<int>(argv[2]);
55  ny = boost::lexical_cast<int>(argv[3]);
56  sx = boost::lexical_cast<double>(argv[4]);
57  sy = boost::lexical_cast<double>(argv[5]);
58  nummodes = boost::lexical_cast<int>(argv[6]);
59  output_file = boost::lexical_cast<string>(argv[7]);
60 
61  map<string, LibUtilities::ShapeType> acceptable_shapes;
62  acceptable_shapes["triangles"] = LibUtilities::eTriangle;
63  acceptable_shapes["quads"] = LibUtilities::eQuadrilateral;
64 
65  map<std::string, LibUtilities::ShapeType>::iterator type = acceptable_shapes.find(type_name);
66  if (type == acceptable_shapes.end())
67  {
68  cerr << "Wrong mesh element type " << type_name << endl << endl;
69  print_usage_info(argv[0]);
70  throw 1;
71  }
72 
73  output.open(output_file.c_str());
74 
75  for (i = 0; i < nx; i++)
76  {
77  xc.push_back( (double)i * (1.0 / (nx-1)) );
78  }
79  for (i = 0; i < ny; i++)
80  {
81  yc.push_back( (double)i * (1.0 / (ny-1)) );
82  }
83 
84  output<< "<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
85  output<< "<NEKTAR>" << endl;
86 
87 
88  output << "<EXPANSIONS>" << endl;
89  output << "<E COMPOSITE=\"C[0]\" NUMMODES=\"" << nummodes << "\" FIELDS=\"u\" TYPE=\"MODIFIED\" />" <<endl;
90  output << "</EXPANSIONS>\n" << endl;
91 
92  PrintConditions(output);
93 
94  //Vertices
95  output << "<GEOMETRY DIM=\"2\" SPACE=\"2\">" << endl;
96 
97  output << " <VERTEX";
98  if (sx != 1.0)
99  {
100  output << " XSCALE=\"" << sx << "\"";
101  }
102  if (sy != 1.0)
103  {
104  output << " YSCALE=\"" << sy << "\"";
105  }
106  output << ">"<< endl;
107 
108  nx = xc.size();
109  ny = yc.size();
110  for(j = 0; j < ny; ++j)
111  {
112  for(k = 0; k < nx; ++k)
113  {
114  output << " <V ID=\"" << j*nx+k << "\">\t";
115  output << xc[k] << " " << yc[j] << " 0.0";
116  output << " </V>" << endl;
117  }
118  }
119 
120  output << " </VERTEX>\n" << endl;
121 
122  // Edges. By default it generates edges for
123  // quadrilateral mesh.
124  output << " <EDGE>" << endl;
125  int cnt = 0;
126  for(j = 0; j < ny-1; ++j)
127  {
128  for(i = 0; i < nx-1; ++i)
129  {
130  output << " <E ID=\"" << cnt++ << "\">\t";
131  output << j*nx+i <<" " << j*nx + i+1 ;
132  output << " </E>" << endl;
133  }
134 
135  for(i = 0; i < nx; ++i)
136  {
137  output << " <E ID=\"" << cnt++ << "\">\t";
138  output << j*nx+i <<" " << (j+1)*nx+i;
139  output << " </E>" << endl;
140  }
141  }
142 
143  for(i = 0; i < nx-1; ++i)
144  {
145  output << " <E ID=\"" << cnt++ << "\">\t";
146  output << j*nx+i <<" " << j*nx + i+1 ;
147  output << " </E>" << endl;
148  }
149 
150  // total number of quad edges. Useful in renumbering
151  // diagonal edges later on.
152  int total_quad_edges = cnt-1;
153 
154  // Triangular mesh is made by adding diagonal segments
155  if (type->second == LibUtilities::eTriangle)
156  {
157  // generating diagonal edges
158 
159  for(j = 0; j < ny-1; ++j)
160  {
161  for(i = 0; i < nx-1; ++i)
162  {
163  output << " <E ID=\"" << cnt++ << "\">\t";
164  output << j*nx+i <<" " << (j+1)*nx + i+1 ;
165  output << " </E>" << endl;
166  }
167  }
168  }
169  output << " </EDGE>\n" << endl;
170 
171  output << " <ELEMENT>" << endl;
172  cnt = 0;
173  switch(type->second)
174  {
175  case LibUtilities::eQuadrilateral:
176  for(j = 0; j < ny-1; ++j)
177  {
178  for(i = 0; i < nx-1; ++i)
179  {
180  output << " <Q ID=\"" << cnt++ << "\">\t";
181  output << j*((nx-1)+nx)+i <<" " << j*((nx-1)+nx) + (nx-1)+i+1 << " " ;
182  output << (j+1)*((nx-1)+nx)+i <<" " << j*((nx-1)+nx) + (nx-1)+i ;
183  output << " </Q>" << endl;
184  }
185  }
186  break;
187  case LibUtilities::eTriangle:
188  for(j = 0; j < ny-1; ++j)
189  {
190  for(i = 0; i < nx-1; ++i)
191  {
192  output << " <T ID=\"" << cnt++ << "\">\t";
193  output << total_quad_edges + 1 + (j*(nx-1)+i) << " " ;
194  output << (j+1)*((nx-1)+nx)+i <<" " << j*((nx-1)+nx) + (nx-1)+i ;
195  output << " </T>" << endl;
196 
197  output << " <T ID=\"" << cnt++ << "\">\t";
198  output << j*((nx-1)+nx)+i <<" " << j*((nx-1)+nx) + (nx-1)+i+1 << " " ;
199  output << total_quad_edges + 1 + (j*(nx-1)+i);
200  output << " </T>" << endl;
201  }
202  }
203  break;
204  default:
205  cerr << "unknown element type\n";
206  throw 1;
207  }
208  output << " </ELEMENT>\n" << endl;
209 
210 
211  output << "<COMPOSITE>" << endl;
212 
213  switch(type->second)
214  {
215  case LibUtilities::eQuadrilateral:
216  output << "<C ID=\"0\"> Q[0-" << (nx-1)*(ny-1)-1 << "] </C>" << endl;
217  break;
218  case LibUtilities::eTriangle:
219  output << "<C ID=\"0\"> T[0-" << 2*(nx-1)*(ny-1)-1 << "] </C>" << endl;
220  break;
221  default:
222  cerr << "unknown element type\n";
223  throw 1;
224  }
225 
226  // boundary composites coincide for both mesh element types
227  output << "<C ID=\"1\"> E[";
228  for(i = 0; i < nx-1; ++i)
229  {
230  output << i;
231  if(i != nx-2)
232  {
233  output << ",";
234  }
235  }
236  output << "] </C> // south border" << endl;
237 
238  output << "<C ID=\"2\"> E[";
239  for(i = 0; i < ny-1; ++i)
240  {
241  output << (nx-1)*(i+1) + nx*i;
242  if(i != ny-2)
243  {
244  output << ",";
245  }
246  }
247  output << "] </C> // west border" << endl;
248 
249  output << "<C ID=\"3\"> E[";
250  for(i = 0; i < nx-1; ++i)
251  {
252  output << (nx-1)*(ny-1)+ nx*(ny-1)+ i;
253  if(i != nx-2)
254  {
255  output << ",";
256  }
257  }
258  output << "] </C> // north border" << endl;
259 
260  output << "<C ID=\"4\"> E[";
261  for(i = 0; i < ny-1; ++i)
262  {
263  output << (nx-1)*(i+1) + nx*i + nx-1;
264  if(i != ny-2)
265  {
266  output << ",";
267  }
268  }
269  output << "] </C> // East border" << endl;
270 
271 
272  output << "</COMPOSITE>\n" << endl;
273 
274 
275  output << "<DOMAIN> C[0] </DOMAIN>\n" << endl;
276  output << "</GEOMETRY>\n" << endl;
277 
278  output << "</NEKTAR>" << endl;
279 
280  }
281  catch(...)
282  {
283  return 1;
284  }
285 
286  return 0;
287 
288 }
289 
290 
291 
292 void PrintConditions(ofstream& output)
293 {
294  output << "<CONDITIONS>" << endl;
295 
296  output << "<SOLVERINFO>" << endl;
297  output << "<I PROPERTY=\"GlobalSysSoln\" VALUE=\"DirectFull\"/>" << endl;
298  output << "</SOLVERINFO>\n" << endl;
299 
300  output << "<PARAMETERS>" << endl;
301  output << "<P> TimeStep = 0.002 </P>" << endl;
302  output << "<P> Lambda = 1 </P>" << endl;
303  output << "</PARAMETERS>\n" << endl;
304 
305  output << "<VARIABLES>" << endl;
306  output << " <V ID=\"0\"> u </V>" << endl;
307  output << "</VARIABLES>\n" << endl;
308 
309  output << "<BOUNDARYREGIONS>" << endl;
310  output << "<B ID=\"0\"> C[1] </B>" << endl;
311  output << "<B ID=\"1\"> C[2] </B>" << endl;
312  output << "<B ID=\"2\"> C[3] </B>" << endl;
313  output << "<B ID=\"3\"> C[4] </B>" << endl;
314  output << "</BOUNDARYREGIONS>\n" << endl;
315 
316  output << "<BOUNDARYCONDITIONS>" << endl;
317  output << " <REGION REF=\"0\"> // South border " << endl;
318  output << " <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
319  output << " </REGION>" << endl;
320 
321  output << " <REGION REF=\"1\"> // West border " << endl;
322  output << " <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
323  output << " </REGION>" << endl;
324 
325  output << " <REGION REF=\"2\"> // North border " << endl;
326  output << " <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
327  output << " </REGION>" << endl;
328 
329  output << " <REGION REF=\"3\"> // East border " << endl;
330  output << " <N VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
331  output << " </REGION>" << endl;
332  output << "</BOUNDARYCONDITIONS>" << endl;
333 
334  output << " <FUNCTION NAME=\"Forcing\">" << endl;
335  output << " <E VAR=\"u\" VALUE=\"-(Lambda + 2*PI*PI/4)*sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
336  output << " </FUNCTION>" << endl;
337 
338  output << " <FUNCTION NAME=\"ExactSolution\">" << endl;
339  output << " <E VAR=\"u\" VALUE=\"sin(PI/2*x)*sin(PI/2*y)\" />" << endl;
340  output << " </FUNCTION>" << endl;
341 
342  output << "</CONDITIONS>" << endl;
343 }
344 
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2