Rheolef  7.2
an efficient C++ finite element environment
mkgeo_ugrid.sh
Go to the documentation of this file.
1 #!/bin/sh
2 #
3 # This file is part of Rheolef.
4 #
5 # Copyright (C) 2000-2009 Pierre Saramito
6 #
7 # Rheolef is free software; you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation; either version 2 of the License, or
10 # (at your option) any later version.
11 #
12 # Rheolef is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with Rheolef; if not, write to the Free Software
19 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 #
21 # -------------------------------------------------------------------------
22 # author: Pierre.Saramito@imag.fr
23 # date: 1 nov 2011
24 
25 
160 
161 # ------------------------------------------
162 # utility
163 # ------------------------------------------
164 verbose=false
165 my_eval () {
166  command="$*"
167  if $verbose; then echo "! $command" 1>&2; fi
168  eval $command
169  if test $? -ne 0; then
170  echo "$0: error on command: $command" >&2
171  exit 1
172  fi
173 }
174 # ------------------------------------------
175 # 1d case
176 # ------------------------------------------
177 mkgmsh_1d () {
178  variant=$1; shift
179  n=$1; shift
180  sides=$1; shift;
181  boundary=$1; shift;
182  corner=$1; shift;
183  a=$1
184  b=$2
185 
186 cat << EOF_1D
187 n = $n;
188 a = $a;
189 b = $b;
190 h = 1.0*(b-a)/n;
191 Point(1) = {a, 0, 0, h};
192 Point(2) = {b, 0, 0, h};
193 Line(3) = {1,2};
194 EOF_1D
195 if $sides || $corner; then
196  echo "Physical Point(\"left\") = {1};"
197  echo "Physical Point(\"right\") = {2};"
198 fi
199 if $boundary; then
200  echo "Physical Point(\"boundary\") = {1,2};"
201 fi
202 echo "Physical Line(\"interior\") = {3};"
203 }
204 # ------------------------------------------
205 # 2d case: t,q,tq variants
206 # ------------------------------------------
207 mkgmsh_2d () {
208  variant=$1; shift
209  n=$1; shift
210  sides=$1; shift;
211  boundary=$1; shift;
212  corner=$1; shift;
213  a=$1
214  b=$2
215  c=$3
216  d=$4
217 
218 echo "n = $n;"
219 echo "a = $a; b = $b;"
220 echo "c = $c; d = $d;"
221 if test $variant = t -o $variant = tq; then
222  echo "h = 1.0*(b-a)/n;"
223 else
224  echo "h = 2.0*(b-a)/n;"
225 fi
226 cat << EOF_2D_1
227 Point(1) = {a, c, 0, h};
228 Point(2) = {b, c, 0, h};
229 Point(3) = {b, d, 0, h};
230 Point(4) = {a, d, 0, h};
231 Line(1) = {1,2};
232 Line(2) = {2,3};
233 Line(3) = {3,4};
234 Line(4) = {4,1};
235 Line Loop(5) = {1,2,3,4};
236 Mesh.Algorithm = 1;
237 Plane Surface(6) = {5} ;
238 EOF_2D_1
239 if test $variant = q; then
240  echo "Mesh.RecombinationAlgorithm = 1;"
241  echo "Mesh.SubdivisionAlgorithm = 1;"
242  echo "Recombine Surface {6};"
243 fi
244 if test $variant = tq; then
245  echo "Mesh.RecombinationAlgorithm = 0;"
246  echo "angle = 90.0;"
247  echo "Recombine Surface {6} = angle;"
248 fi
249 if $corner; then
250  echo "Physical Point(\"left_bottom\") = {1};"
251  echo "Physical Point(\"right_bottom\") = {2};"
252  echo "Physical Point(\"right_top\") = {3};"
253  echo "Physical Point(\"left_top\") = {4};"
254 fi
255 if $boundary; then
256  echo "Physical Line(\"boundary\") = {1,2,3,4};"
257 fi
258 if $sides; then
259  echo "Physical Line(\"bottom\") = {1};"
260  echo "Physical Line(\"right\") = {2};"
261  echo "Physical Line(\"top\") = {3};"
262  echo "Physical Line(\"left\") = {4};"
263 fi
264 echo "Physical Surface(\"interior\") = {6};"
265 }
266 # ------------------------------------------
267 # 3d case: T,TP,TPH variants
268 # ------------------------------------------
269 mkgmsh_3d_T_TP_TPH () {
270  variant=$1; shift
271  n=$1; shift
272  sides=$1; shift;
273  boundary=$1; shift;
274  corner=$1; shift;
275  a=$1
276  b=$2
277  c=$3
278  d=$4
279  f=$5
280  g=$6
281 
282 cat << EOF_3D_T_1
283 n = $n;
284 a = $a; c = $c; f = $f;
285 b = $b; d = $d; g = $g;
286 EOF_3D_T_1
287 if test $variant = T; then
288  d=" d"
289  g=" g"
290 elif test $variant = TP; then
291  echo "fg = (f+g)/2.0;"
292  d=" d"
293  g="fg"
294 elif test $variant = TPH; then
295  echo "cd = (c+d)/2.0;"
296  echo "fg = (f+g)/2.0;"
297  d="cd"
298  g="fg"
299 fi
300 cat << EOF_3D_T_2
301 h = (b-a)/n;
302 Point(1) = {a, c, f, h};
303 Point(2) = {b, c, f, h};
304 Point(3) = {b, $d, f, h};
305 Point(4) = {a, $d, f, h};
306 Point(5) = {a, c, $g, h};
307 Point(6) = {b, c, $g, h};
308 Point(7) = {b, $d, $g, h};
309 Point(8) = {a, $d, $g, h};
310 Line(1) = {1,2};
311 Line(2) = {2,3};
312 Line(3) = {3,4};
313 Line(4) = {4,1};
314 Line(5) = {5,6};
315 Line(6) = {6,7};
316 Line(7) = {7,8};
317 Line(8) = {8,5};
318 Line(9) = {1,5};
319 Line(10) = {2,6};
320 Line(11) = {3,7};
321 Line(12) = {4,8};
322 Line Loop(21) = {-1,-4,-3,-2};
323 Line Loop(22) = {5,6,7,8};
324 Line Loop(23) = {1,10,-5,-9};
325 Line Loop(25) = {2,11,-6,-10};
326 Line Loop(24) = {12,-7,-11,3};
327 Line Loop(26) = {9,-8,-12,4};
328 Plane Surface(31) = {21} ;
329 Plane Surface(32) = {22} ;
330 Plane Surface(33) = {23} ;
331 Plane Surface(34) = {24} ;
332 Plane Surface(35) = {25} ;
333 Plane Surface(36) = {26} ;
334 Surface Loop(41) = {31,32,33,34,35,36};
335 Volume(51) = {41};
336 EOF_3D_T_2
337 if test $variant = TP; then
338  echo "extr[] = Extrude{0,0,(g-f)/2}{ Surface{32}; Layers{n/2}; Recombine;};"
339 elif test $variant = TPH; then
340  echo "extr_y[] = Extrude{0,(d-c)/2,0}{ Surface{34}; Layers{n/2}; Recombine;};"
341  echo "extr_z[] = Extrude{0,0,(g-f)/2}{ Surface{32,extr_y[3]}; Layers{n/2}; Recombine;};"
342 fi
343 if test $variant = T; then
344  if $boundary; then
345  echo "Physical Surface(\"boundary\") = {31, 32, 33, 34, 35, 36};"
346  fi
347  if $sides; then
348  echo "Physical Surface(\"bottom\") = {31};"
349  echo "Physical Surface(\"top\") = {32};"
350  echo "Physical Surface(\"left\") = {33};"
351  echo "Physical Surface(\"front\") = {35};"
352  echo "Physical Surface(\"right\") = {34};"
353  echo "Physical Surface(\"back\") = {36};"
354  fi
355  echo "Physical Volume(\"internal\") = {51};"
356 elif test $variant = TP; then
357  if $boundary; then
358  echo "Physical Surface(\"boundary\") = {31, extr[0], 33, extr[2], 35, extr[3], 34, extr[4], 36, extr[5]};"
359  fi
360  if $sides; then
361  echo "Physical Surface(\"bottom\") = {31};"
362  echo "Physical Surface(\"top\") = {extr[0]};"
363  echo "Physical Surface(\"left\") = {33, extr[2]};"
364  echo "Physical Surface(\"front\") = {35, extr[3]};"
365  echo "Physical Surface(\"right\") = {34, extr[4]};"
366  echo "Physical Surface(\"back\") = {36, extr[5]};"
367  fi
368  echo "Physical Volume(\"internal\") = {51,extr[1]};"
369 elif test $variant = TPH; then
370  if $boundary; then
371  echo "Physical Surface(\"boundary\") = {31, extr_y[5], extr_z[0], extr_z[6],"
372  echo " 33,extr_z[2], 35, extr_y[4], extr_z[3], extr_z[9],"
373  echo " extr_y[0],extr_z[10], 36, extr_y[2], extr_z[5], extr_z[11]};"
374  fi
375  if $sides; then
376  echo "Physical Surface(\"bottom\") = {31, extr_y[5]};"
377  echo "Physical Surface(\"top\") = {extr_z[0], extr_z[6]};"
378  echo "Physical Surface(\"left\") = {33,extr_z[2]};"
379  echo "Physical Surface(\"front\") = {35, extr_y[4], extr_z[3], extr_z[9]};"
380  echo "Physical Surface(\"right\") = {extr_y[0],extr_z[10]};"
381  echo "Physical Surface(\"back\") = {36, extr_y[2], extr_z[5], extr_z[11]};"
382  fi
383  echo "Physical Volume(\"internal\") = {51,extr_y[1],extr_z[1],extr_z[7]};"
384 fi
385 }
386 # ------------------------------------------
387 # 3d case: P,H,PH variants
388 # ------------------------------------------
389 mkgmsh_3d_P_H_PH () {
390  variant=$1; shift
391  n=$1; shift
392  sides=$1; shift;
393  boundary=$1; shift;
394  corner=$1; shift;
395  a=$1
396  b=$2
397  c=$3
398  d=$4
399  f=$5
400  g=$6
401 
402 if test $variant = H -a $n -ne 1; then
403  n="$n/2";
404 fi
405 cat << EOF_3D_PH
406 n = $n;
407 a = $a; c = $c; f = $f;
408 b = $b; d = $d; g = $g;
409 h = 1.0*(b-a)/n;
410 Point(1) = {a, c, f, h};
411 Point(2) = {b, c, f, h};
412 Point(6) = {b, c, g, h};
413 Point(5) = {a, c, g, h};
414 Line(1) = {1,2};
415 Line(6) = {2,6};
416 Line(9) = {5,6};
417 Line(5) = {1,5};
418 Line Loop(300) = {1,6,-9,-5};
419 Mesh.Algorithm = 1;
420 EOF_3D_PH
421 if test $variant = P; then
422  echo "Plane Surface(3) = {300};"
423 elif test $variant = H; then
424  echo "Plane Surface(3) = {300};"
425  echo "Mesh.SubdivisionAlgorithm = 1;"
426  echo "Mesh.RecombinationAlgorithm = 1;"
427  echo "Recombine Surface {3};"
428  echo "Mesh.SubdivisionAlgorithm = 2;"
429 elif test $variant = PH; then
430  echo "Plane Surface(3) = {300};"
431  echo "Mesh.RecombinationAlgorithm = 0;"
432  echo "angle = 45.0;"
433  echo "Recombine Surface {3} = angle;"
434 fi
435 echo "extr[] = Extrude{0,d-c,0}{ Surface{3}; Layers{n}; Recombine;};"
436 if $boundary; then
437  echo "Physical Surface(\"boundary\") = {-extr[2], -extr[4], 3, -extr[3], -extr[0], -extr[5]};"
438 fi
439 if $sides; then
440  echo "Physical Surface(\"bottom\") = {-extr[2]};"
441  echo "Physical Surface(\"top\") = {-extr[4]};"
442  echo "Physical Surface(\"left\") = {3};"
443  echo "Physical Surface(\"front\") = {-extr[3]};"
444  echo "Physical Surface(\"right\") = {-extr[0]};"
445  echo "Physical Surface(\"back\") = {-extr[5]};"
446 fi
447 echo "Physical Volume(\"interior\") = {extr[1]};"
448 }
449 #----------------------------------------------
450 # multi-dim switch
451 #----------------------------------------------
452 mkgmsh () {
453 dim=$1; shift
454 case $dim in
455  1) mkgmsh_1d $*;;
456  2) mkgmsh_2d $*;;
457  *) variant=$1
458  case $variant in
459  T*) mkgmsh_3d_T_TP_TPH $*;;
460  *) mkgmsh_3d_P_H_PH $*;;
461  esac
462 esac
463 }
464 #----------------------------------------------
465 # main
466 #----------------------------------------------
467 usage="mkgeo_ugrid
468  [-{eptqTPH}|-tq|-TP|-PH|-TPH]
469  [n]
470  [-{abcdfg} float]
471  [-[no]sides]
472  [-[no]boundary]
473  [-[no]corner]
474  [-order int]
475  [-[no]clean]
476  [-[no]verbose]
477 "
478 
479 if test $# -eq 0; then
480  echo ${usage} >&2
481  exit 0
482 fi
483 
484 clean=true
485 bindir=""
486 dim=2
487 variant=t
488 n=10
489 a=0; b=1
490 c=0; d=1
491 f=0; g=1
492 boundary=true
493 sides=true
494 corner=false
495 order=1
496 while test $# -ne 0; do
497  case $1 in
498  -h) echo ${usage} >&2; exit 0;;
499  -e) dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
500  -[tq]|-tq) dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
501  -[TPH]|-TP|-PH|-TPH) dim=3; variant=`echo x$1 | sed -e 's/x-//'`;;
502  -[abcdfg]) var=`echo x$1 | sed -e 's/x-//'`
503  if test x"$2" = x""; then echo ${usage} >&2; exit 1; fi
504  expr="$var=$2"
505  eval $expr
506  shift
507  ;;
508  [0-9]*) n=$1;;
509  -order) order=$2; shift;;
510  -sides) sides=true;;
511  -nosides) sides=false;;
512  -boundary) boundary=true;;
513  -noboundary) boundary=false;;
514  -corner) corner=true;;
515  -nocorner) corner=false;;
516  -verbose) verbose=true;;
517  -noverbose) verbose=false;;
518  -clean) clean=true;;
519  -noclean) clean=false;;
520  -bindir)
521  if test x"$2" = x""; then echo ${usage} >&2; exit 1; fi
522  bindir="$2/"
523  shift
524  ;;
525  *) echo ${usage} >&2; exit 1;;
526  esac
527  shift
528 done
529 if $clean; then
530  tmp="/tmp/tmp$$"
531 else
532  tmp="output"
533 fi
534 #echo "dim=$dim" 1>&2
535 #echo "variant=$variant" 1>&2
536 #echo "n=$n" 1>&2
537 #echo "a=$a" 1>&2
538 #echo "b=$b" 1>&2
539 mkgmsh $dim $variant $n $sides $boundary $corner $a $b $c $d $f $g > $tmp.mshcad
540 my_eval "gmsh -$dim -order $order -format msh2 $tmp.mshcad -o $tmp.msh > $tmp.log"
541 if test ! -f $tmp.msh; then
542  echo "$0: gmsh failed"
543  exit 1
544 fi
545 MSH2GEO="${bindir}msh2geo"
546 GEO="${bindir}geo"
547 my_eval "$MSH2GEO $tmp.msh"
548 # upgrade: no more need with msh2geo
549 #my_eval "$MSH2GEO $tmp.msh > ${tmp}-v1.geo"
550 #my_eval "$GEO -upgrade -geo - < ${tmp}-v1.geo"
551 if $clean; then
552  my_eval "rm -f $tmp.mshcad $tmp.log $tmp.msh"
553 fi
554