Rheolef  7.1
an efficient C++ finite element environment
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 
129 nx=""
130 nx_default=1
131 N=100
132 c=2
133 L="unset"
134 hmin=0.05
135 split=false
136 quarter=1
137 clean=true
138 verbose=true
139 name="obstacle"
140 usage="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 "
153 pkgbindir="`rheolef-config --pkglibdir`"
154 
155 while 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
176 done
177 if test x"$nx" = x""; then
178  nx=$nx_default
179 fi
180 if test x"$ny" = x""; then
181  ny=$nx
182 fi
183 if test $L = "unset"; then
184  L=`echo | awk -v c=$c '{print c < 2 ? 2 : 2*c}'`
185 fi
186 c_error=`echo | awk -v c=$c '{print (c <= 1 ? "true" : "false")}'`
187 if $c_error; then
188  echo "$0: invalid c=$c" 1>&2 ; exit 1
189 fi
190 to_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 #
216 awk -v L=$L -v c=$c -v N=$N -v quarter=$quarter '
217 BEGIN {
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 
272 if test $quarter -eq 0; then
273 cat > $name.dmn << EOF1
274 EdgeDomainNames
275  5
276  axis
277  obstacle
278  downstream
279  wall
280  upstream
281 EOF1
282 else
283 cat > $name.dmn << EOF2
284 EdgeDomainNames
285  5
286  obstacle
287  axis
288  downstream
289  wall
290  vaxis
291 EOF2
292 fi
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 # -------------------------------------
300 awk -v L=$L -v c=$c -v N=$N -v hmin=$hmin -v nx=$nx -v ny=$ny -v quarter=$quarter '
301 BEGIN {
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 
343 to_clean="$to_clean $name.mtr"
344 
345 command="bamg -g $name.bamgcad -M $name.mtr -o $name.bamg"
346 if $verbose; then
347  command="$command 1>&2"
348 else
349  command="$command 1> $name.bamglog"
350  to_clean="$to_clean $name.bamglog"
351 fi
352 $verbose && echo "! $command" 1>&2
353 eval $command
354 status=$?
355 if 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
359 fi
360 echo "! $name.bamg created" 1>&2
361 
362 if $split; then
363  filter="| ${pkgbindir}/geo_split | geo -upgrade -geo -"
364 else
365  filter=""
366 fi
367 command="bamg2geo $name.bamg $name.dmn $sys_coord_opt $filter > $name.geo"
368 $verbose && echo "! $command" 1>&2
369 eval $command
370 status=$?
371 if test $status -ne 0; then
372  echo "$0: command failed" 1>&2
373  exit $status
374 fi
375 echo "! $name.geo created" 1>&2
376 if $clean; then
377  command="rm -f $to_clean"
378  $verbose && echo "! $command" 1>&2
379  eval $command
380 fi