Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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# ------------------------------------------
164verbose=false
165my_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# ------------------------------------------
177mkgmsh_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
186cat << EOF_1D
187n = $n;
188a = $a;
189b = $b;
190h = 1.0*(b-a)/n;
191Point(1) = {a, 0, 0, h};
192Point(2) = {b, 0, 0, h};
193Line(3) = {1,2};
194EOF_1D
195if $sides || $corner; then
196 echo "Physical Point(\"left\") = {1};"
197 echo "Physical Point(\"right\") = {2};"
198fi
199if $boundary; then
200 echo "Physical Point(\"boundary\") = {1,2};"
201fi
202echo "Physical Line(\"interior\") = {3};"
203}
204# ------------------------------------------
205# 2d case: t,q,tq variants
206# ------------------------------------------
207mkgmsh_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
218echo "n = $n;"
219echo "a = $a; b = $b;"
220echo "c = $c; d = $d;"
221if test $variant = t -o $variant = tq; then
222 echo "h = 1.0*(b-a)/n;"
223else
224 echo "h = 2.0*(b-a)/n;"
225fi
226cat << EOF_2D_1
227Point(1) = {a, c, 0, h};
228Point(2) = {b, c, 0, h};
229Point(3) = {b, d, 0, h};
230Point(4) = {a, d, 0, h};
231Line(1) = {1,2};
232Line(2) = {2,3};
233Line(3) = {3,4};
234Line(4) = {4,1};
235Line Loop(5) = {1,2,3,4};
236Mesh.Algorithm = 1;
237Plane Surface(6) = {5} ;
238EOF_2D_1
239if test $variant = q; then
240 echo "Mesh.RecombinationAlgorithm = 1;"
241 echo "Mesh.SubdivisionAlgorithm = 1;"
242 echo "Recombine Surface {6};"
243fi
244if test $variant = tq; then
245 echo "Mesh.RecombinationAlgorithm = 0;"
246 echo "angle = 90.0;"
247 echo "Recombine Surface {6} = angle;"
248fi
249if $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};"
254fi
255if $boundary; then
256 echo "Physical Line(\"boundary\") = {1,2,3,4};"
257fi
258if $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};"
263fi
264echo "Physical Surface(\"interior\") = {6};"
265}
266# ------------------------------------------
267# 3d case: T,TP,TPH variants
268# ------------------------------------------
269mkgmsh_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
282cat << EOF_3D_T_1
283n = $n;
284a = $a; c = $c; f = $f;
285b = $b; d = $d; g = $g;
286EOF_3D_T_1
287if test $variant = T; then
288 d=" d"
289 g=" g"
290elif test $variant = TP; then
291 echo "fg = (f+g)/2.0;"
292 d=" d"
293 g="fg"
294elif test $variant = TPH; then
295 echo "cd = (c+d)/2.0;"
296 echo "fg = (f+g)/2.0;"
297 d="cd"
298 g="fg"
299fi
300cat << EOF_3D_T_2
301h = (b-a)/n;
302Point(1) = {a, c, f, h};
303Point(2) = {b, c, f, h};
304Point(3) = {b, $d, f, h};
305Point(4) = {a, $d, f, h};
306Point(5) = {a, c, $g, h};
307Point(6) = {b, c, $g, h};
308Point(7) = {b, $d, $g, h};
309Point(8) = {a, $d, $g, h};
310Line(1) = {1,2};
311Line(2) = {2,3};
312Line(3) = {3,4};
313Line(4) = {4,1};
314Line(5) = {5,6};
315Line(6) = {6,7};
316Line(7) = {7,8};
317Line(8) = {8,5};
318Line(9) = {1,5};
319Line(10) = {2,6};
320Line(11) = {3,7};
321Line(12) = {4,8};
322Line Loop(21) = {-1,-4,-3,-2};
323Line Loop(22) = {5,6,7,8};
324Line Loop(23) = {1,10,-5,-9};
325Line Loop(25) = {2,11,-6,-10};
326Line Loop(24) = {12,-7,-11,3};
327Line Loop(26) = {9,-8,-12,4};
328Plane Surface(31) = {21} ;
329Plane Surface(32) = {22} ;
330Plane Surface(33) = {23} ;
331Plane Surface(34) = {24} ;
332Plane Surface(35) = {25} ;
333Plane Surface(36) = {26} ;
334Surface Loop(41) = {31,32,33,34,35,36};
335Volume(51) = {41};
336EOF_3D_T_2
337if test $variant = TP; then
338 echo "extr[] = Extrude{0,0,(g-f)/2}{ Surface{32}; Layers{n/2}; Recombine;};"
339elif 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;};"
342fi
343if 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};"
356elif 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]};"
369elif 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]};"
384fi
385}
386# ------------------------------------------
387# 3d case: P,H,PH variants
388# ------------------------------------------
389mkgmsh_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
402if test $variant = H -a $n -ne 1; then
403 n="$n/2";
404fi
405cat << EOF_3D_PH
406n = $n;
407a = $a; c = $c; f = $f;
408b = $b; d = $d; g = $g;
409h = 1.0*(b-a)/n;
410Point(1) = {a, c, f, h};
411Point(2) = {b, c, f, h};
412Point(6) = {b, c, g, h};
413Point(5) = {a, c, g, h};
414Line(1) = {1,2};
415Line(6) = {2,6};
416Line(9) = {5,6};
417Line(5) = {1,5};
418Line Loop(300) = {1,6,-9,-5};
419Mesh.Algorithm = 1;
420EOF_3D_PH
421if test $variant = P; then
422 echo "Plane Surface(3) = {300};"
423elif 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;"
429elif 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;"
434fi
435echo "extr[] = Extrude{0,d-c,0}{ Surface{3}; Layers{n}; Recombine;};"
436if $boundary; then
437 echo "Physical Surface(\"boundary\") = {-extr[2], -extr[4], 3, -extr[3], -extr[0], -extr[5]};"
438fi
439if $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]};"
446fi
447echo "Physical Volume(\"interior\") = {extr[1]};"
448}
449#----------------------------------------------
450# multi-dim switch
451#----------------------------------------------
452mkgmsh () {
453dim=$1; shift
454case $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
462esac
463}
464#----------------------------------------------
465# main
466#----------------------------------------------
467usage="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
479if test $# -eq 0; then
480 echo ${usage} >&2
481 exit 0
482fi
483
484clean=true
485bindir=""
486dim=2
487variant=t
488n=10
489a=0; b=1
490c=0; d=1
491f=0; g=1
492boundary=true
493sides=true
494corner=false
495order=1
496while 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
528done
529if $clean; then
530 tmp="/tmp/tmp$$"
531else
532 tmp="output"
533fi
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
539mkgmsh $dim $variant $n $sides $boundary $corner $a $b $c $d $f $g > $tmp.mshcad
540my_eval "gmsh -$dim -order $order -format msh2 $tmp.mshcad -o $tmp.msh > $tmp.log"
541if test ! -f $tmp.msh; then
542 echo "$0: gmsh failed"
543 exit 1
544fi
545MSH2GEO="${bindir}msh2geo"
546GEO="${bindir}geo"
547my_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"
551if $clean; then
552 my_eval "rm -f $tmp.mshcad $tmp.log $tmp.msh"
553fi
554