mkgeo_ball.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 test "$verbose" = true; then echo "! $command" 1>&2; fi
eval $command
if test $? -ne 0; then
echo "$0: error on command: $command" >&2
exit 1
fi
}
# ------------------------------------------
# 2d case: t,q,tq variants
# ------------------------------------------
mkgmsh_2d () {
surface_only=$1; shift
variant=$1; shift
n=$1; shift
a=$1
b=$2
c=$3
d=$4
cat << EOF_2D_1
Mesh.SecondOrderIncomplete = 0;
/* Mesh.ElementOrder = $order; */
/*
Mesh.SecondOrderLinear = 0;
Mesh.SecondOrderExperimental = 0;
Mesh.SmoothInternalEdges = 100;
Mesh.Smoothing = 100;
*/
n = $n; // the density of discretisation
a = $a; b = $b;
c = $c; d = $d;
EOF_2D_1
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_2
ab = (a+b)/2;
cd = (c+d)/2;
Point(40) = {ab, cd, 0, h};
Point(41) = { b, cd, 0, h};
Point(42) = {ab, d, 0, h};
Point(43) = { a, cd, 0, h};
Point(44) = {ab, c, 0, h};
Circle(21) = {41, 40, 42};
Circle(22) = {42, 40, 43};
Circle(23) = {43, 40, 44};
Circle(24) = {44, 40, 41};
Line Loop (51) = {21,22,23,24};
Physical Line("boundary") = {21,22,23,24};
EOF_2D_2
if test "${surface_only}" != true; then
echo "Mesh.Algorithm = 1;"
echo "Plane Surface(61) = {51};"
echo "Physical Surface(\"interior\") = {61};"
if test $variant = q; then
echo "Mesh.RecombinationAlgorithm = 1;" ;
# Note: Mesh.RecombinationAlgorithm = 0 : standard => triangle with 2 bdry edges
# Mesh.RecombinationAlgorithm = 1 : blossom => triangle with 2 bdry edges, too
# => mkgeo_ball_gmsh_fix fails
# => try a higher mesh density... sometimes its ok...
echo "Mesh.SubdivisionAlgorithm = 1;"
echo "Recombine Surface {61};"
fi
if test $variant = tq; then
echo "Mesh.RecombinationAlgorithm = 0;"
echo "angle = 90.0;"
echo "Recombine Surface {61} = angle;"
fi
fi
}
# ------------------------------------------
# 3d case
# http://sites.google.com/site/auxcapucins/maillage-3d-en-gmsh---maillage-d-une-sphere
# ------------------------------------------
mkgmsh_3d () {
surface_only=$1; shift
variant=$1; shift
n=$1; shift
a=$1
b=$2
c=$3
d=$4
f=$5
g=$6
cat << EOF_3D_1
Mesh.SecondOrderIncomplete = 0;
n = $n; // the density of discretisation
a = $a; b = $b;
c = $c; d = $d;
f = $f; g = $g;
ab = (a+b)/2;
cd = (c+d)/2;
fg = (f+g)/2;
EOF_3D_1
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_3D_2
Point(1) = {ab, cd, fg, h};
Point(2) = { b, cd, fg, h};
Point(3) = {ab, d, fg, h};
Point(4) = {ab, cd, g, h};
Point(5) = {a, cd, fg, h};
Point(6) = {ab, c, fg, h};
Point(7) = {ab, cd, f, h};
Circle(1) = {2,1,3}; // {start,center,end}
Circle(2) = {3,1,5};
Circle(3) = {5,1,6};
Circle(4) = {6,1,2};
Circle(5) = {2,1,7};
Circle(6) = {7,1,5};
Circle(7) = {5,1,4};
Circle(8) = {4,1,2};
Circle(9) = {6,1,7};
Circle(10) = {7,1,3};
Circle(11) = {3,1,4};
Circle(12) = {4,1,6};
Line Loop(1) = {1,11,8};
Line Loop(2) = {2,7,-11};
Line Loop(3) = {3,-12,-7};
Line Loop(4) = {4,-8,12};
Line Loop(5) = {5,10,-1};
Line Loop(6) = {-2,-10,6};
Line Loop(7) = {-3,-6,-9};
Line Loop(8) = {-4,9,-5};
Mesh.Algorithm = 1;
Ruled Surface(1) = {1};
Ruled Surface(2) = {2};
Ruled Surface(3) = {3};
Ruled Surface(4) = {4};
Ruled Surface(5) = {5};
Ruled Surface(6) = {6};
Ruled Surface(7) = {7};
Ruled Surface(8) = {8};
Surface Loop (1) = {1,2,3,4,5,6,7,8};
Physical Surface("boundary") = {1,2,3,4,5,6,7,8};
EOF_3D_2
if test "${surface_only}" = true; then
if test $variant = q; then
echo "Mesh.RecombinationAlgorithm = 1;"
echo "Mesh.SubdivisionAlgorithm = 1;"
echo "Recombine Surface {1};"
fi
if test $variant = tq; then
echo "Mesh.RecombinationAlgorithm = 0;"
echo "angle = 90.0;"
echo "Recombine Surface {1} = angle;"
fi
else
echo "Mesh.Algorithm3D = 1; // attention: Tetgen is not available in debian packages"
echo "Mesh.OptimizeNetgen = 1;"
echo "Mesh.Optimize = 1;"
echo "Mesh.Smoothing = 10;"
echo "Volume (1) = {1};"
if test $variant = H; then
echo "Mesh.RecombinationAlgorithm = 1;"
echo "Mesh.SubdivisionAlgorithm = 2;"
fi
echo "Physical Volume(\"internal\") = {1};"
fi
}
#----------------------------------------------
# multi-dim switch
#----------------------------------------------
mkgmsh () {
dim=$1; shift
case $dim in
2) mkgmsh_2d $*;;
*) mkgmsh_3d $*;;
esac
}
#----------------------------------------------
# main
#----------------------------------------------
usage="mkgeo_ball
[-{tq}|-tq]
[n]
[-s]
[-order k]
[-{abcdfg} float]
[-[no]fix]
[-[no]clean]
[-[no]verbose]
"
if test $# -eq 0; then
echo ${usage} >&2
exit 0
fi
GMSH=${GMSH-"gmsh"}
pkgbindir=`rheolef-config --pkglibdir`
verbose=false
clean=true
bindir=""
map_dim=2
variant=t
n=10
a=-1; b=1
c=-1; d=1
f=-1; g=1
order=1
surface_only="false"
fix=true
while test $# -ne 0; do
case $1 in
-h) echo ${usage} >&2; exit 0;;
-clean) clean=true;;
-noclean) clean=false; verbose=true;;
-verbose) verbose=true;;
-noverbose) verbose=false;;
-fix) fix=true;;
-nofix) fix=false;;
-s) surface_only=true;;
-e) map_dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
-[tq]|-tq) map_dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
-[TPH]|-TP|-PH|-TPH) map_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 "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
expr="$var=$2"
eval $expr
shift
;;
-order) if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
order=$2
shift
;;
[0-9]*) n=$1;;
-bindir)
if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
bindir="$2/"
shift
;;
*) echo "$0: invalid option: $1" >&2; echo ${usage} >&2; exit 1;;
esac
shift
done
if $clean; then
tmp="/tmp/tmp$$"
else
tmp="output"
fi
if test "${surface_only}" = true; then
dim=`expr ${map_dim} + 1`
else
dim=${map_dim}
fi
#echo "dim=${dim}" 1>&2
#echo "map_dim=${map_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} ${surface_only} $variant $n $a $b $c $d $f $g > $tmp.mshcad
if $fix; then
# mkgeo_ball_gmsh_fix will set mesh order:
my_eval "$GMSH -${map_dim} $tmp.mshcad -format msh2 -o $tmp.msh > $tmp.log 2>&1"
else
my_eval "$GMSH -${map_dim} -order $order $tmp.mshcad -format msh2 -o $tmp.msh > $tmp.log 2>&1"
fi
if test ! -f $tmp.msh; then
echo "$0: gmsh failed"
exit 1
fi
MSH2GEO="${bindir}msh2geo"
GEO="${bindir}geo"
my_eval "$MSH2GEO $tmp.msh > ${tmp}-v2.geo"
# upgrade: no more need with msh2geo
# my_eval "$GEO -upgrade -geo -check - < ${tmp}-v1.geo > ${tmp}-v2.geo"
if $fix; then
fix_filter="$pkgbindir/mkgeo_ball_gmsh_fix -order $order"
else
fix_filter="cat"
fi
my_eval "$fix_filter < ${tmp}-v2.geo"
if $clean; then
my_eval "rm -f $tmp.log ${tmp}-v1.geo ${tmp}-v2.geo"
my_eval "rm -f $tmp.mshcad $tmp.msh"
fi