Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
mkgeo_obstacle.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: 2 february 2018
24
25
128
129nx=""
130nx_default=1
131N=100
132c=2
133L="unset"
134hmin=0.05
135split=false
136quarter=1
137clean=true
138verbose=true
139name="obstacle"
140usage="usage: mkgeo_obstacle
141 [nx=$nx_default [ny=nx]]
142 [-c float=$c]
143 [-L float=$L]
144 [-N float=$N]
145 [-quarter|-half]
146 [-rz]
147 [-hmin float=$hmin]
148 [-name string=$name]
149 [-[no]split]
150 [-[no]clean]
151 [-[no]verbose]
152"
153pkgbindir="`rheolef-config --pkglibdir`"
154
155while test $# -ne 0; do
156 case $1 in
157 [0-9]*) if test x$nx = x""; then nx=$1; else ny=$1; fi;;
158 -c) c="$2"; shift;;
159 -L) L="$2"; shift;;
160 -N) N="$2"; shift;;
161 -quarter) quarter=1;;
162 -half) quarter=0;;
163 -zr|-cartesian) sys_coord_opt="$1";;
164 -name) name=$2; shift;;
165 -hmin) hmin=$2; shift;;
166 -split) split=true;;
167 -nosplit) split=false;;
168 -clean) clean=true;;
169 -noclean) clean=false;;
170 -verbose) verbose=true;;
171 -noverbose) verbose=false;;
172 -h) /bin/echo -E ${usage} >&2; exit 0;;
173 *) /bin/echo -E ${usage} >&2; exit 1;;
174 esac
175 shift
176done
177if test x"$nx" = x""; then
178 nx=$nx_default
179fi
180if test x"$ny" = x""; then
181 ny=$nx
182fi
183if test $L = "unset"; then
184 L=`echo | awk -v c=$c '{print c < 2 ? 2 : 2*c}'`
185fi
186c_error=`echo | awk -v c=$c '{print (c <= 1 ? "true" : "false")}'`
187if $c_error; then
188 echo "$0: invalid c=$c" 1>&2 ; exit 1
189fi
190to_clean=""
191
192#
193# background mesh for bamg mesh generator:
194# with anisotropic (1/hx^2, 1/hy^2) metric
195#
196# half geometry:
197# +-------------------+-------------------+ y=c
198# |6+N 5+N 4+N|
199# | |
200# | --- | y=1
201# | / \ |
202# |1 2| |2+N 3+N|
203# +----------------+ + +----------------+ y=0
204# x=-L x=-1 x=0 x=1 x=L
205#
206# quarter geometry:
207# +-------------------+ y=c
208# 4+N 3+N|
209# | |
210# 1+- | y=1
211# \ |
212# |1+N 2+N|
213# + +----------------+ y=0
214# x=0 x=1 x=L
215#
216awk -v L=$L -v c=$c -v N=$N -v quarter=$quarter '
217BEGIN {
218 pi = 3.14159265358979323846
219 print "MeshVersionFormatted"
220 print " 0"
221 print "Dimension"
222 print " 2"
223 print "Vertices"
224 print " ", (!quarter ? 6+N : 4+N)
225 if (! quarter) {
226 printf (" %.16g %.16g %d\n", -L, 0, 1)
227 printf (" %.16g %.16g %d\n", -1, 0, 2)
228 for (i = 1; i <= N; ++i) {
229 theta = pi - pi*(1.0*i/N)
230 printf (" %.16g %.16g %d\n", cos(theta), sin(theta), 2+i)
231 }
232 printf (" %.16g %.16g %d\n", L, 0, 3+N)
233 printf (" %.16g %.16g %d\n", L, c, 4+N)
234 printf (" %.16g %.16g %d\n", 0, c, 5+N)
235 printf (" %.16g %.16g %d\n", -L, c, 6+N)
236 } else {
237 printf (" %.16g %.16g %d\n", 0, 1, 1)
238 for (i = 1; i <= N; ++i) {
239 theta = pi/2 - pi/2*(1.0*i/N)
240 printf (" %.16g %.16g %d\n", cos(theta), sin(theta), 1+i)
241 }
242 printf (" %.16g %.16g %d\n", L, 0, 2+N)
243 printf (" %.16g %.16g %d\n", L, c, 3+N)
244 printf (" %.16g %.16g %d\n", 0, c, 4+N)
245 }
246 print "Edges"
247 print " ", (!quarter ? 6+N : 4+N)
248 if (! quarter) {
249 print " 1 2 101"
250 for (i = 0; i < N; ++i) {
251 print " ", 2+i, 2+i+1, 102
252 }
253 print " ", 2+N, 3+N, 101
254 print " ", 3+N, 4+N, 103
255 print " ", 4+N, 5+N, 104
256 print " ", 5+N, 6+N, 104
257 print " ", 6+N, 1, 105
258 } else {
259 for (i = 0; i < N; ++i) {
260 print " ", 1+i, 1+i+1, 101
261 }
262 print " ", 1+N, 2+N, 102
263 print " ", 2+N, 3+N, 103
264 print " ", 3+N, 4+N, 104
265 print " ", 4+N, 1, 105
266 }
267}
268' \
269> $name.bamgcad
270$verbose && echo "! $name.bamgcad created" 1>&2
271
272if test $quarter -eq 0; then
273cat > $name.dmn << EOF1
274EdgeDomainNames
275 5
276 axis
277 obstacle
278 downstream
279 wall
280 upstream
281EOF1
282else
283cat > $name.dmn << EOF2
284EdgeDomainNames
285 5
286 obstacle
287 axis
288 downstream
289 wall
290 vaxis
291EOF2
292fi
293$verbose && echo "! $name.dmn created" 1>&2
294
295#echo "! nx=$nx ny=$ny hmin=$hmin" 1>&2
296
297# -------------------------------------
298# metric for adaptive mesh
299# -------------------------------------
300awk -v L=$L -v c=$c -v N=$N -v hmin=$hmin -v nx=$nx -v ny=$ny -v quarter=$quarter '
301BEGIN {
302 hm = 0.1
303 hx = hm*L/nx
304 hy = hm*c/ny
305 h0 = hmin # (hmin >= 0.5*(c-1)) ? 0.5*(c-1) : hmin
306 hf = (c-1)*h0
307 h1 = (c-1 < 1) ? hf : h0
308 m0 = h0**(-2)
309 m1 = h1**(-2)
310 mf = hf**(-2)
311 mx = hx**(-2)
312 my = hy**(-2)
313 print (!quarter ? 6+N : 4+N), 3
314 if (!quarter) {
315 printf (" %.16g %.16g %.16g\n", mx, 0, my)
316 for (i = 0; i <= N; ++i) {
317 # TODO: N -> N/2
318 theta = (i < N/2) ? 1.0*i/(0.5*N) : 1.0*(N-i)/(0.5*N);
319 hi = h1*theta + h0*(1-theta)
320 mi = hi**(-2)
321 printf (" %.16g %.16g %.16g\n", mi, 0, mi)
322 }
323 printf (" %.16g %.16g %.16g\n", mx, 0, my)
324 printf (" %.16g %.16g %.16g\n", mx, 0, my)
325 printf (" %.16g %.16g %.16g\n", mf, 0, mf)
326 printf (" %.16g %.16g %.16g\n", mx, 0, my)
327 } else {
328 for (i = 0; i <= N; ++i) {
329 theta = 1.0*(N-i)/N;
330 hi = h1*theta + h0*(1-theta)
331 mi = hi**(-2)
332 printf (" %.16g %.16g %.16g\n", mi, 0, mi)
333 }
334 printf (" %.16g %.16g %.16g\n", mx, 0, my)
335 printf (" %.16g %.16g %.16g\n", mx, 0, my)
336 printf (" %.16g %.16g %.16g\n", mf, 0, mf)
337 }
338}
339' \
340> $name.mtr
341$verbose && echo "! $name.mtr created" 1>&2
342
343to_clean="$to_clean $name.mtr"
344
345command="bamg -g $name.bamgcad -M $name.mtr -o $name.bamg"
346if $verbose; then
347 command="$command 1>&2"
348else
349 command="$command 1> $name.bamglog"
350 to_clean="$to_clean $name.bamglog"
351fi
352$verbose && echo "! $command" 1>&2
353eval $command
354status=$?
355if test $status -ne 0; then
356 if $verbose; then true; else cat $name.bamglog 1>&2; fi
357 echo "$0: command failed" 1>&2
358 exit $status
359fi
360echo "! $name.bamg created" 1>&2
361
362if $split; then
363 filter="| ${pkgbindir}/geo_split | geo -upgrade -geo -"
364else
365 filter=""
366fi
367command="bamg2geo $name.bamg $name.dmn $sys_coord_opt $filter > $name.geo"
368$verbose && echo "! $command" 1>&2
369eval $command
370status=$?
371if test $status -ne 0; then
372 echo "$0: command failed" 1>&2
373 exit $status
374fi
375echo "! $name.geo created" 1>&2
376if $clean; then
377 command="rm -f $to_clean"
378 $verbose && echo "! $command" 1>&2
379 eval $command
380fi