8 #include <boost/lexical_cast.hpp>
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";
33 int main(
int argc,
char *argv[])
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]);
61 map<string, LibUtilities::ShapeType> acceptable_shapes;
66 if (type == acceptable_shapes.end())
68 cerr <<
"Wrong mesh element type " << type_name << endl << endl;
73 output.open(output_file.c_str());
75 for (i = 0; i < nx; i++)
77 xc.push_back( (
double)i * (1.0 / (nx-1)) );
79 for (i = 0; i < ny; i++)
81 yc.push_back( (
double)i * (1.0 / (ny-1)) );
84 output<<
"<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
85 output<<
"<NEKTAR>" << endl;
88 output <<
"<EXPANSIONS>" << endl;
89 output <<
"<E COMPOSITE=\"C[0]\" NUMMODES=\"" << nummodes <<
"\" FIELDS=\"u\" TYPE=\"MODIFIED\" />" <<endl;
90 output <<
"</EXPANSIONS>\n" << endl;
95 output <<
"<GEOMETRY DIM=\"2\" SPACE=\"2\">" << endl;
100 output <<
" XSCALE=\"" << sx <<
"\"";
104 output <<
" YSCALE=\"" << sy <<
"\"";
106 output <<
">"<< endl;
110 for(j = 0; j < ny; ++j)
112 for(k = 0; k < nx; ++k)
114 output <<
" <V ID=\"" << j*nx+k <<
"\">\t";
115 output << xc[k] <<
" " << yc[j] <<
" 0.0";
116 output <<
" </V>" << endl;
120 output <<
" </VERTEX>\n" << endl;
124 output <<
" <EDGE>" << endl;
126 for(j = 0; j < ny-1; ++j)
128 for(i = 0; i < nx-1; ++i)
130 output <<
" <E ID=\"" << cnt++ <<
"\">\t";
131 output << j*nx+i <<
" " << j*nx + i+1 ;
132 output <<
" </E>" << endl;
135 for(i = 0; i < nx; ++i)
137 output <<
" <E ID=\"" << cnt++ <<
"\">\t";
138 output << j*nx+i <<
" " << (j+1)*nx+i;
139 output <<
" </E>" << endl;
143 for(i = 0; i < nx-1; ++i)
145 output <<
" <E ID=\"" << cnt++ <<
"\">\t";
146 output << j*nx+i <<
" " << j*nx + i+1 ;
147 output <<
" </E>" << endl;
152 int total_quad_edges = cnt-1;
159 for(j = 0; j < ny-1; ++j)
161 for(i = 0; i < nx-1; ++i)
163 output <<
" <E ID=\"" << cnt++ <<
"\">\t";
164 output << j*nx+i <<
" " << (j+1)*nx + i+1 ;
165 output <<
" </E>" << endl;
169 output <<
" </EDGE>\n" << endl;
171 output <<
" <ELEMENT>" << endl;
176 for(j = 0; j < ny-1; ++j)
178 for(i = 0; i < nx-1; ++i)
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;
188 for(j = 0; j < ny-1; ++j)
190 for(i = 0; i < nx-1; ++i)
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;
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;
205 cerr <<
"unknown element type\n";
208 output <<
" </ELEMENT>\n" << endl;
211 output <<
"<COMPOSITE>" << endl;
216 output <<
"<C ID=\"0\"> Q[0-" << (nx-1)*(ny-1)-1 <<
"] </C>" << endl;
219 output <<
"<C ID=\"0\"> T[0-" << 2*(nx-1)*(ny-1)-1 <<
"] </C>" << endl;
222 cerr <<
"unknown element type\n";
227 output <<
"<C ID=\"1\"> E[";
228 for(i = 0; i < nx-1; ++i)
236 output <<
"] </C> // south border" << endl;
238 output <<
"<C ID=\"2\"> E[";
239 for(i = 0; i < ny-1; ++i)
241 output << (nx-1)*(i+1) + nx*i;
247 output <<
"] </C> // west border" << endl;
249 output <<
"<C ID=\"3\"> E[";
250 for(i = 0; i < nx-1; ++i)
252 output << (nx-1)*(ny-1)+ nx*(ny-1)+ i;
258 output <<
"] </C> // north border" << endl;
260 output <<
"<C ID=\"4\"> E[";
261 for(i = 0; i < ny-1; ++i)
263 output << (nx-1)*(i+1) + nx*i + nx-1;
269 output <<
"] </C> // East border" << endl;
272 output <<
"</COMPOSITE>\n" << endl;
275 output <<
"<DOMAIN> C[0] </DOMAIN>\n" << endl;
276 output <<
"</GEOMETRY>\n" << endl;
278 output <<
"</NEKTAR>" << endl;
294 output <<
"<CONDITIONS>" << endl;
296 output <<
"<SOLVERINFO>" << endl;
297 output <<
"<I PROPERTY=\"GlobalSysSoln\" VALUE=\"DirectFull\"/>" << endl;
298 output <<
"</SOLVERINFO>\n" << endl;
300 output <<
"<PARAMETERS>" << endl;
301 output <<
"<P> TimeStep = 0.002 </P>" << endl;
302 output <<
"<P> Lambda = 1 </P>" << endl;
303 output <<
"</PARAMETERS>\n" << endl;
305 output <<
"<VARIABLES>" << endl;
306 output <<
" <V ID=\"0\"> u </V>" << endl;
307 output <<
"</VARIABLES>\n" << endl;
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;
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;
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;
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;
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;
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;
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;
342 output <<
"</CONDITIONS>" << endl;
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void PrintConditions(ofstream &output)
void print_usage_info(char *binary_name)
int main(int argc, char *argv[])