mkgeo_ugrid.sh
Go to the documentation of this file.
#!/bin/sh
#
# This file is part of Rheolef.
#
# Copyright (C) 2000-2009 Pierre Saramito
#
# Rheolef is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# Rheolef is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Rheolef; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# -------------------------------------------------------------------------
# author: Pierre.Saramito@imag.fr
# date: 1 nov 2011
# ------------------------------------------
# utility
# ------------------------------------------
verbose=false
my_eval () {
command="$*"
if $verbose; then echo "! $command" 1>&2; fi
eval $command
if test $? -ne 0; then
echo "$0: error on command: $command" >&2
exit 1
fi
}
# ------------------------------------------
# 1d case
# ------------------------------------------
mkgmsh_1d () {
variant=$1; shift
n=$1; shift
sides=$1; shift;
boundary=$1; shift;
corner=$1; shift;
a=$1
b=$2
cat << EOF_1D
n = $n;
a = $a;
b = $b;
h = 1.0*(b-a)/n;
Point(1) = {a, 0, 0, h};
Point(2) = {b, 0, 0, h};
Line(3) = {1,2};
EOF_1D
if $sides || $corner; then
echo "Physical Point(\"left\") = {1};"
echo "Physical Point(\"right\") = {2};"
fi
if $boundary; then
echo "Physical Point(\"boundary\") = {1,2};"
fi
echo "Physical Line(\"interior\") = {3};"
}
# ------------------------------------------
# 2d case: t,q,tq variants
# ------------------------------------------
mkgmsh_2d () {
variant=$1; shift
n=$1; shift
sides=$1; shift;
boundary=$1; shift;
corner=$1; shift;
a=$1
b=$2
c=$3
d=$4
echo "n = $n;"
echo "a = $a; b = $b;"
echo "c = $c; d = $d;"
if test $variant = t -o $variant = tq; then
echo "h = 1.0*(b-a)/n;"
else
echo "h = 2.0*(b-a)/n;"
fi
cat << EOF_2D_1
Point(1) = {a, c, 0, h};
Point(2) = {b, c, 0, h};
Point(3) = {b, d, 0, h};
Point(4) = {a, d, 0, h};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line Loop(5) = {1,2,3,4};
Mesh.Algorithm = 1;
Plane Surface(6) = {5} ;
EOF_2D_1
if test $variant = q; then
echo "Mesh.RecombinationAlgorithm = 1;"
echo "Mesh.SubdivisionAlgorithm = 1;"
echo "Recombine Surface {6};"
fi
if test $variant = tq; then
echo "Mesh.RecombinationAlgorithm = 0;"
echo "angle = 90.0;"
echo "Recombine Surface {6} = angle;"
fi
if $corner; then
echo "Physical Point(\"left_bottom\") = {1};"
echo "Physical Point(\"right_bottom\") = {2};"
echo "Physical Point(\"right_top\") = {3};"
echo "Physical Point(\"left_top\") = {4};"
fi
if $boundary; then
echo "Physical Line(\"boundary\") = {1,2,3,4};"
fi
if $sides; then
echo "Physical Line(\"bottom\") = {1};"
echo "Physical Line(\"right\") = {2};"
echo "Physical Line(\"top\") = {3};"
echo "Physical Line(\"left\") = {4};"
fi
echo "Physical Surface(\"interior\") = {6};"
}
# ------------------------------------------
# 3d case: T,TP,TPH variants
# ------------------------------------------
mkgmsh_3d_T_TP_TPH () {
variant=$1; shift
n=$1; shift
sides=$1; shift;
boundary=$1; shift;
corner=$1; shift;
a=$1
b=$2
c=$3
d=$4
f=$5
g=$6
cat << EOF_3D_T_1
n = $n;
a = $a; c = $c; f = $f;
b = $b; d = $d; g = $g;
EOF_3D_T_1
if test $variant = T; then
d=" d"
g=" g"
elif test $variant = TP; then
echo "fg = (f+g)/2.0;"
d=" d"
g="fg"
elif test $variant = TPH; then
echo "cd = (c+d)/2.0;"
echo "fg = (f+g)/2.0;"
d="cd"
g="fg"
fi
cat << EOF_3D_T_2
h = (b-a)/n;
Point(1) = {a, c, f, h};
Point(2) = {b, c, f, h};
Point(3) = {b, $d, f, h};
Point(4) = {a, $d, f, h};
Point(5) = {a, c, $g, h};
Point(6) = {b, c, $g, h};
Point(7) = {b, $d, $g, h};
Point(8) = {a, $d, $g, h};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};
Line(9) = {1,5};
Line(10) = {2,6};
Line(11) = {3,7};
Line(12) = {4,8};
Line Loop(21) = {-1,-4,-3,-2};
Line Loop(22) = {5,6,7,8};
Line Loop(23) = {1,10,-5,-9};
Line Loop(25) = {2,11,-6,-10};
Line Loop(24) = {12,-7,-11,3};
Line Loop(26) = {9,-8,-12,4};
Plane Surface(31) = {21} ;
Plane Surface(32) = {22} ;
Plane Surface(33) = {23} ;
Plane Surface(34) = {24} ;
Plane Surface(35) = {25} ;
Plane Surface(36) = {26} ;
Surface Loop(41) = {31,32,33,34,35,36};
Volume(51) = {41};
EOF_3D_T_2
if test $variant = TP; then
echo "extr[] = Extrude{0,0,(g-f)/2}{ Surface{32}; Layers{n/2}; Recombine;};"
elif test $variant = TPH; then
echo "extr_y[] = Extrude{0,(d-c)/2,0}{ Surface{34}; Layers{n/2}; Recombine;};"
echo "extr_z[] = Extrude{0,0,(g-f)/2}{ Surface{32,extr_y[3]}; Layers{n/2}; Recombine;};"
fi
if test $variant = T; then
if $boundary; then
echo "Physical Surface(\"boundary\") = {31, 32, 33, 34, 35, 36};"
fi
if $sides; then
echo "Physical Surface(\"bottom\") = {31};"
echo "Physical Surface(\"top\") = {32};"
echo "Physical Surface(\"left\") = {33};"
echo "Physical Surface(\"front\") = {35};"
echo "Physical Surface(\"right\") = {34};"
echo "Physical Surface(\"back\") = {36};"
fi
echo "Physical Volume(\"internal\") = {51};"
elif test $variant = TP; then
if $boundary; then
echo "Physical Surface(\"boundary\") = {31, extr[0], 33, extr[2], 35, extr[3], 34, extr[4], 36, extr[5]};"
fi
if $sides; then
echo "Physical Surface(\"bottom\") = {31};"
echo "Physical Surface(\"top\") = {extr[0]};"
echo "Physical Surface(\"left\") = {33, extr[2]};"
echo "Physical Surface(\"front\") = {35, extr[3]};"
echo "Physical Surface(\"right\") = {34, extr[4]};"
echo "Physical Surface(\"back\") = {36, extr[5]};"
fi
echo "Physical Volume(\"internal\") = {51,extr[1]};"
elif test $variant = TPH; then
if $boundary; then
echo "Physical Surface(\"boundary\") = {31, extr_y[5], extr_z[0], extr_z[6],"
echo " 33,extr_z[2], 35, extr_y[4], extr_z[3], extr_z[9],"
echo " extr_y[0],extr_z[10], 36, extr_y[2], extr_z[5], extr_z[11]};"
fi
if $sides; then
echo "Physical Surface(\"bottom\") = {31, extr_y[5]};"
echo "Physical Surface(\"top\") = {extr_z[0], extr_z[6]};"
echo "Physical Surface(\"left\") = {33,extr_z[2]};"
echo "Physical Surface(\"front\") = {35, extr_y[4], extr_z[3], extr_z[9]};"
echo "Physical Surface(\"right\") = {extr_y[0],extr_z[10]};"
echo "Physical Surface(\"back\") = {36, extr_y[2], extr_z[5], extr_z[11]};"
fi
echo "Physical Volume(\"internal\") = {51,extr_y[1],extr_z[1],extr_z[7]};"
fi
}
# ------------------------------------------
# 3d case: P,H,PH variants
# ------------------------------------------
mkgmsh_3d_P_H_PH () {
variant=$1; shift
n=$1; shift
sides=$1; shift;
boundary=$1; shift;
corner=$1; shift;
a=$1
b=$2
c=$3
d=$4
f=$5
g=$6
if test $variant = H -a $n -ne 1; then
n="$n/2";
fi
cat << EOF_3D_PH
n = $n;
a = $a; c = $c; f = $f;
b = $b; d = $d; g = $g;
h = 1.0*(b-a)/n;
Point(1) = {a, c, f, h};
Point(2) = {b, c, f, h};
Point(6) = {b, c, g, h};
Point(5) = {a, c, g, h};
Line(1) = {1,2};
Line(6) = {2,6};
Line(9) = {5,6};
Line(5) = {1,5};
Line Loop(300) = {1,6,-9,-5};
Mesh.Algorithm = 1;
EOF_3D_PH
if test $variant = P; then
echo "Plane Surface(3) = {300};"
elif test $variant = H; then
echo "Plane Surface(3) = {300};"
echo "Mesh.SubdivisionAlgorithm = 1;"
echo "Mesh.RecombinationAlgorithm = 1;"
echo "Recombine Surface {3};"
echo "Mesh.SubdivisionAlgorithm = 2;"
elif test $variant = PH; then
echo "Plane Surface(3) = {300};"
echo "Mesh.RecombinationAlgorithm = 0;"
echo "angle = 45.0;"
echo "Recombine Surface {3} = angle;"
fi
echo "extr[] = Extrude{0,d-c,0}{ Surface{3}; Layers{n}; Recombine;};"
if $boundary; then
echo "Physical Surface(\"boundary\") = {-extr[2], -extr[4], 3, -extr[3], -extr[0], -extr[5]};"
fi
if $sides; then
echo "Physical Surface(\"bottom\") = {-extr[2]};"
echo "Physical Surface(\"top\") = {-extr[4]};"
echo "Physical Surface(\"left\") = {3};"
echo "Physical Surface(\"front\") = {-extr[3]};"
echo "Physical Surface(\"right\") = {-extr[0]};"
echo "Physical Surface(\"back\") = {-extr[5]};"
fi
echo "Physical Volume(\"interior\") = {extr[1]};"
}
#----------------------------------------------
# multi-dim switch
#----------------------------------------------
mkgmsh () {
dim=$1; shift
case $dim in
1) mkgmsh_1d $*;;
2) mkgmsh_2d $*;;
*) variant=$1
case $variant in
T*) mkgmsh_3d_T_TP_TPH $*;;
*) mkgmsh_3d_P_H_PH $*;;
esac
esac
}
#----------------------------------------------
# main
#----------------------------------------------
usage="mkgeo_ugrid
[-{eptqTPH}|-tq|-TP|-PH|-TPH]
[n]
[-{abcdfg} float]
[-[no]sides]
[-[no]boundary]
[-[no]corner]
[-order int]
[-[no]clean]
[-[no]verbose]
"
if test $# -eq 0; then
echo ${usage} >&2
exit 0
fi
clean=true
bindir=""
dim=2
variant=t
n=10
a=0; b=1
c=0; d=1
f=0; g=1
boundary=true
sides=true
corner=false
order=1
while test $# -ne 0; do
case $1 in
-h) echo ${usage} >&2; exit 0;;
-e) dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
-[tq]|-tq) dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
-[TPH]|-TP|-PH|-TPH) dim=3; variant=`echo x$1 | sed -e 's/x-//'`;;
-[abcdfg]) var=`echo x$1 | sed -e 's/x-//'`
if test x"$2" = x""; then echo ${usage} >&2; exit 1; fi
expr="$var=$2"
eval $expr
shift
;;
[0-9]*) n=$1;;
-order) order=$2; shift;;
-sides) sides=true;;
-nosides) sides=false;;
-boundary) boundary=true;;
-noboundary) boundary=false;;
-corner) corner=true;;
-nocorner) corner=false;;
-verbose) verbose=true;;
-noverbose) verbose=false;;
-clean) clean=true;;
-noclean) clean=false;;
-bindir)
if test x"$2" = x""; then echo ${usage} >&2; exit 1; fi
bindir="$2/"
shift
;;
*) echo ${usage} >&2; exit 1;;
esac
shift
done
if $clean; then
tmp="/tmp/tmp$$"
else
tmp="output"
fi
#echo "dim=$dim" 1>&2
#echo "variant=$variant" 1>&2
#echo "n=$n" 1>&2
#echo "a=$a" 1>&2
#echo "b=$b" 1>&2
mkgmsh $dim $variant $n $sides $boundary $corner $a $b $c $d $f $g > $tmp.mshcad
my_eval "gmsh -$dim -order $order -format msh2 $tmp.mshcad -o $tmp.msh > $tmp.log"
if test ! -f $tmp.msh; then
echo "$0: gmsh failed"
exit 1
fi
MSH2GEO="${bindir}msh2geo"
GEO="${bindir}geo"
my_eval "$MSH2GEO $tmp.msh"
# upgrade: no more need with msh2geo
#my_eval "$MSH2GEO $tmp.msh > ${tmp}-v1.geo"
#my_eval "$GEO -upgrade -geo - < ${tmp}-v1.geo"
if $clean; then
my_eval "rm -f $tmp.mshcad $tmp.log $tmp.msh"
fi