Loading...
Searching...
No Matches
mkgeo_obstacle.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: 2 february 2018
nx=""
nx_default=1
N=100
c=2
L="unset"
hmin=0.05
split=false
quarter=1
clean=true
verbose=true
name="obstacle"
usage="usage: mkgeo_obstacle
[nx=$nx_default [ny=nx]]
[-c float=$c]
[-L float=$L]
[-N float=$N]
[-quarter|-half]
[-rz]
[-hmin float=$hmin]
[-name string=$name]
[-[no]split]
[-[no]clean]
[-[no]verbose]
"
pkgbindir="`rheolef-config --pkglibdir`"
while test $# -ne 0; do
case $1 in
[0-9]*) if test x$nx = x""; then nx=$1; else ny=$1; fi;;
-c) c="$2"; shift;;
-L) L="$2"; shift;;
-N) N="$2"; shift;;
-quarter) quarter=1;;
-half) quarter=0;;
-zr|-cartesian) sys_coord_opt="$1";;
-name) name=$2; shift;;
-hmin) hmin=$2; shift;;
-split) split=true;;
-nosplit) split=false;;
-clean) clean=true;;
-noclean) clean=false;;
-verbose) verbose=true;;
-noverbose) verbose=false;;
-h) /bin/echo -E ${usage} >&2; exit 0;;
*) /bin/echo -E ${usage} >&2; exit 1;;
esac
shift
done
if test x"$nx" = x""; then
nx=$nx_default
fi
if test x"$ny" = x""; then
ny=$nx
fi
if test $L = "unset"; then
L=`echo | awk -v c=$c '{print c < 2 ? 2 : 2*c}'`
fi
c_error=`echo | awk -v c=$c '{print (c <= 1 ? "true" : "false")}'`
if $c_error; then
echo "$0: invalid c=$c" 1>&2 ; exit 1
fi
to_clean=""
#
# background mesh for bamg mesh generator:
# with anisotropic (1/hx^2, 1/hy^2) metric
#
# half geometry:
# +-------------------+-------------------+ y=c
# |6+N 5+N 4+N|
# | |
# | --- | y=1
# | / \ |
# |1 2| |2+N 3+N|
# +----------------+ + +----------------+ y=0
# x=-L x=-1 x=0 x=1 x=L
#
# quarter geometry:
# +-------------------+ y=c
# 4+N 3+N|
# | |
# 1+- | y=1
# \ |
# |1+N 2+N|
# + +----------------+ y=0
# x=0 x=1 x=L
#
awk -v L=$L -v c=$c -v N=$N -v quarter=$quarter '
BEGIN {
pi = 3.14159265358979323846
print "MeshVersionFormatted"
print " 0"
print "Dimension"
print " 2"
print "Vertices"
print " ", (!quarter ? 6+N : 4+N)
if (! quarter) {
printf (" %.16g %.16g %d\n", -L, 0, 1)
printf (" %.16g %.16g %d\n", -1, 0, 2)
for (i = 1; i <= N; ++i) {
theta = pi - pi*(1.0*i/N)
printf (" %.16g %.16g %d\n", cos(theta), sin(theta), 2+i)
}
printf (" %.16g %.16g %d\n", L, 0, 3+N)
printf (" %.16g %.16g %d\n", L, c, 4+N)
printf (" %.16g %.16g %d\n", 0, c, 5+N)
printf (" %.16g %.16g %d\n", -L, c, 6+N)
} else {
printf (" %.16g %.16g %d\n", 0, 1, 1)
for (i = 1; i <= N; ++i) {
theta = pi/2 - pi/2*(1.0*i/N)
printf (" %.16g %.16g %d\n", cos(theta), sin(theta), 1+i)
}
printf (" %.16g %.16g %d\n", L, 0, 2+N)
printf (" %.16g %.16g %d\n", L, c, 3+N)
printf (" %.16g %.16g %d\n", 0, c, 4+N)
}
print "Edges"
print " ", (!quarter ? 6+N : 4+N)
if (! quarter) {
print " 1 2 101"
for (i = 0; i < N; ++i) {
print " ", 2+i, 2+i+1, 102
}
print " ", 2+N, 3+N, 101
print " ", 3+N, 4+N, 103
print " ", 4+N, 5+N, 104
print " ", 5+N, 6+N, 104
print " ", 6+N, 1, 105
} else {
for (i = 0; i < N; ++i) {
print " ", 1+i, 1+i+1, 101
}
print " ", 1+N, 2+N, 102
print " ", 2+N, 3+N, 103
print " ", 3+N, 4+N, 104
print " ", 4+N, 1, 105
}
}
' \
> $name.bamgcad
$verbose && echo "! $name.bamgcad created" 1>&2
if test $quarter -eq 0; then
cat > $name.dmn << EOF1
EdgeDomainNames
5
axis
obstacle
downstream
wall
upstream
EOF1
else
cat > $name.dmn << EOF2
EdgeDomainNames
5
obstacle
axis
downstream
wall
vaxis
EOF2
fi
$verbose && echo "! $name.dmn created" 1>&2
#echo "! nx=$nx ny=$ny hmin=$hmin" 1>&2
# -------------------------------------
# metric for adaptive mesh
# -------------------------------------
awk -v L=$L -v c=$c -v N=$N -v hmin=$hmin -v nx=$nx -v ny=$ny -v quarter=$quarter '
BEGIN {
hm = 0.1
hx = hm*L/nx
hy = hm*c/ny
h0 = hmin # (hmin >= 0.5*(c-1)) ? 0.5*(c-1) : hmin
hf = (c-1)*h0
h1 = (c-1 < 1) ? hf : h0
m0 = h0**(-2)
m1 = h1**(-2)
mf = hf**(-2)
mx = hx**(-2)
my = hy**(-2)
print (!quarter ? 6+N : 4+N), 3
if (!quarter) {
printf (" %.16g %.16g %.16g\n", mx, 0, my)
for (i = 0; i <= N; ++i) {
# TODO: N -> N/2
theta = (i < N/2) ? 1.0*i/(0.5*N) : 1.0*(N-i)/(0.5*N);
hi = h1*theta + h0*(1-theta)
mi = hi**(-2)
printf (" %.16g %.16g %.16g\n", mi, 0, mi)
}
printf (" %.16g %.16g %.16g\n", mx, 0, my)
printf (" %.16g %.16g %.16g\n", mx, 0, my)
printf (" %.16g %.16g %.16g\n", mf, 0, mf)
printf (" %.16g %.16g %.16g\n", mx, 0, my)
} else {
for (i = 0; i <= N; ++i) {
theta = 1.0*(N-i)/N;
hi = h1*theta + h0*(1-theta)
mi = hi**(-2)
printf (" %.16g %.16g %.16g\n", mi, 0, mi)
}
printf (" %.16g %.16g %.16g\n", mx, 0, my)
printf (" %.16g %.16g %.16g\n", mx, 0, my)
printf (" %.16g %.16g %.16g\n", mf, 0, mf)
}
}
' \
> $name.mtr
$verbose && echo "! $name.mtr created" 1>&2
to_clean="$to_clean $name.mtr"
command="bamg -g $name.bamgcad -M $name.mtr -o $name.bamg"
if $verbose; then
command="$command 1>&2"
else
command="$command 1> $name.bamglog"
to_clean="$to_clean $name.bamglog"
fi
$verbose && echo "! $command" 1>&2
eval $command
status=$?
if test $status -ne 0; then
if $verbose; then true; else cat $name.bamglog 1>&2; fi
echo "$0: command failed" 1>&2
exit $status
fi
echo "! $name.bamg created" 1>&2
if $split; then
filter="| ${pkgbindir}/geo_split | geo -upgrade -geo -"
else
filter=""
fi
command="bamg2geo $name.bamg $name.dmn $sys_coord_opt $filter > $name.geo"
$verbose && echo "! $command" 1>&2
eval $command
status=$?
if test $status -ne 0; then
echo "$0: command failed" 1>&2
exit $status
fi
echo "! $name.geo created" 1>&2
if $clean; then
command="rm -f $to_clean"
$verbose && echo "! $command" 1>&2
eval $command
fi