Rheolef  7.2
an efficient C++ finite element environment
mkgeo_ball.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: # 1 nov 2011
24 
25 
129 
130 # ------------------------------------------
131 # utility
132 # ------------------------------------------
133 verbose=false
134 
135 my_eval () {
136  command="$*"
137  if test "$verbose" = true; then echo "! $command" 1>&2; fi
138  eval $command
139  if test $? -ne 0; then
140  echo "$0: error on command: $command" >&2
141  exit 1
142  fi
143 }
144 # ------------------------------------------
145 # 2d case: t,q,tq variants
146 # ------------------------------------------
147 mkgmsh_2d () {
148  surface_only=$1; shift
149  variant=$1; shift
150  n=$1; shift
151  a=$1
152  b=$2
153  c=$3
154  d=$4
155 
156 cat << EOF_2D_1
157 Mesh.SecondOrderIncomplete = 0;
158 /* Mesh.ElementOrder = $order; */
159 /*
160 Mesh.SecondOrderLinear = 0;
161 Mesh.SecondOrderExperimental = 0;
162 Mesh.SmoothInternalEdges = 100;
163 Mesh.Smoothing = 100;
164 */
165 n = $n; // the density of discretisation
166 a = $a; b = $b;
167 c = $c; d = $d;
168 EOF_2D_1
169 if test $variant = t -o $variant = tq; then
170  echo "h = 1.0*(b-a)/n;"
171 else
172  echo "h = 2.0*(b-a)/n;"
173 fi
174 cat << EOF_2D_2
175 ab = (a+b)/2;
176 cd = (c+d)/2;
177 Point(40) = {ab, cd, 0, h};
178 Point(41) = { b, cd, 0, h};
179 Point(42) = {ab, d, 0, h};
180 Point(43) = { a, cd, 0, h};
181 Point(44) = {ab, c, 0, h};
182 Circle(21) = {41, 40, 42};
183 Circle(22) = {42, 40, 43};
184 Circle(23) = {43, 40, 44};
185 Circle(24) = {44, 40, 41};
186 Line Loop (51) = {21,22,23,24};
187 Physical Line("boundary") = {21,22,23,24};
188 EOF_2D_2
189 if test "${surface_only}" != true; then
190  echo "Mesh.Algorithm = 1;"
191  echo "Plane Surface(61) = {51};"
192  echo "Physical Surface(\"interior\") = {61};"
193  if test $variant = q; then
194  echo "Mesh.RecombinationAlgorithm = 1;" ;
195  # Note: Mesh.RecombinationAlgorithm = 0 : standard => triangle with 2 bdry edges
196  # Mesh.RecombinationAlgorithm = 1 : blossom => triangle with 2 bdry edges, too
197  # => mkgeo_ball_gmsh_fix fails
198  # => try a higher mesh density... sometimes its ok...
199  echo "Mesh.SubdivisionAlgorithm = 1;"
200  echo "Recombine Surface {61};"
201  fi
202  if test $variant = tq; then
203  echo "Mesh.RecombinationAlgorithm = 0;"
204  echo "angle = 90.0;"
205  echo "Recombine Surface {61} = angle;"
206  fi
207 fi
208 }
209 # ------------------------------------------
210 # 3d case
211 # http://sites.google.com/site/auxcapucins/maillage-3d-en-gmsh---maillage-d-une-sphere
212 # ------------------------------------------
213 mkgmsh_3d () {
214  surface_only=$1; shift
215  variant=$1; shift
216  n=$1; shift
217  a=$1
218  b=$2
219  c=$3
220  d=$4
221  f=$5
222  g=$6
223 
224 cat << EOF_3D_1
225 Mesh.SecondOrderIncomplete = 0;
226 n = $n; // the density of discretisation
227 a = $a; b = $b;
228 c = $c; d = $d;
229 f = $f; g = $g;
230 ab = (a+b)/2;
231 cd = (c+d)/2;
232 fg = (f+g)/2;
233 EOF_3D_1
234 if test $variant = t -o $variant = tq; then
235  echo "h = 1.0*(b-a)/n;"
236 else
237  echo "h = 2.0*(b-a)/n;"
238 fi
239 cat << EOF_3D_2
240 Point(1) = {ab, cd, fg, h};
241 Point(2) = { b, cd, fg, h};
242 Point(3) = {ab, d, fg, h};
243 Point(4) = {ab, cd, g, h};
244 Point(5) = {a, cd, fg, h};
245 Point(6) = {ab, c, fg, h};
246 Point(7) = {ab, cd, f, h};
247 Circle(1) = {2,1,3}; // {start,center,end}
248 Circle(2) = {3,1,5};
249 Circle(3) = {5,1,6};
250 Circle(4) = {6,1,2};
251 Circle(5) = {2,1,7};
252 Circle(6) = {7,1,5};
253 Circle(7) = {5,1,4};
254 Circle(8) = {4,1,2};
255 Circle(9) = {6,1,7};
256 Circle(10) = {7,1,3};
257 Circle(11) = {3,1,4};
258 Circle(12) = {4,1,6};
259 Line Loop(1) = {1,11,8};
260 Line Loop(2) = {2,7,-11};
261 Line Loop(3) = {3,-12,-7};
262 Line Loop(4) = {4,-8,12};
263 Line Loop(5) = {5,10,-1};
264 Line Loop(6) = {-2,-10,6};
265 Line Loop(7) = {-3,-6,-9};
266 Line Loop(8) = {-4,9,-5};
267 Mesh.Algorithm = 1;
268 Ruled Surface(1) = {1};
269 Ruled Surface(2) = {2};
270 Ruled Surface(3) = {3};
271 Ruled Surface(4) = {4};
272 Ruled Surface(5) = {5};
273 Ruled Surface(6) = {6};
274 Ruled Surface(7) = {7};
275 Ruled Surface(8) = {8};
276 Surface Loop (1) = {1,2,3,4,5,6,7,8};
277 Physical Surface("boundary") = {1,2,3,4,5,6,7,8};
278 EOF_3D_2
279 if test "${surface_only}" = true; then
280  if test $variant = q; then
281  echo "Mesh.RecombinationAlgorithm = 1;"
282  echo "Mesh.SubdivisionAlgorithm = 1;"
283  echo "Recombine Surface {1};"
284  fi
285  if test $variant = tq; then
286  echo "Mesh.RecombinationAlgorithm = 0;"
287  echo "angle = 90.0;"
288  echo "Recombine Surface {1} = angle;"
289  fi
290 else
291  echo "Mesh.Algorithm3D = 1; // attention: Tetgen is not available in debian packages"
292  echo "Mesh.OptimizeNetgen = 1;"
293  echo "Mesh.Optimize = 1;"
294  echo "Mesh.Smoothing = 10;"
295  echo "Volume (1) = {1};"
296  if test $variant = H; then
297  echo "Mesh.RecombinationAlgorithm = 1;"
298  echo "Mesh.SubdivisionAlgorithm = 2;"
299  fi
300  echo "Physical Volume(\"internal\") = {1};"
301 fi
302 }
303 #----------------------------------------------
304 # multi-dim switch
305 #----------------------------------------------
306 mkgmsh () {
307 dim=$1; shift
308 case $dim in
309  2) mkgmsh_2d $*;;
310  *) mkgmsh_3d $*;;
311 esac
312 }
313 #----------------------------------------------
314 # main
315 #----------------------------------------------
316 usage="mkgeo_ball
317  [-{tq}|-tq]
318  [n]
319  [-s]
320  [-order k]
321  [-{abcdfg} float]
322  [-[no]fix]
323  [-[no]clean]
324  [-[no]verbose]
325 "
326 
327 if test $# -eq 0; then
328  echo ${usage} >&2
329  exit 0
330 fi
331 
332 GMSH=${GMSH-"gmsh"}
333 pkgbindir=`rheolef-config --pkglibdir`
334 verbose=false
335 clean=true
336 bindir=""
337 map_dim=2
338 variant=t
339 n=10
340 a=-1; b=1
341 c=-1; d=1
342 f=-1; g=1
343 order=1
344 surface_only="false"
345 fix=true
346 while test $# -ne 0; do
347  case $1 in
348  -h) echo ${usage} >&2; exit 0;;
349  -clean) clean=true;;
350  -noclean) clean=false; verbose=true;;
351  -verbose) verbose=true;;
352  -noverbose) verbose=false;;
353  -fix) fix=true;;
354  -nofix) fix=false;;
355  -s) surface_only=true;;
356  -e) map_dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
357  -[tq]|-tq) map_dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
358  -[TPH]|-TP|-PH|-TPH) map_dim=3; variant=`echo x$1 | sed -e 's/x-//'`;;
359  -[abcdfg]) var=`echo x$1 | sed -e 's/x-//'`
360  if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
361  expr="$var=$2"
362  eval $expr
363  shift
364  ;;
365  -order) if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
366  order=$2
367  shift
368  ;;
369  [0-9]*) n=$1;;
370  -bindir)
371  if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
372  bindir="$2/"
373  shift
374  ;;
375  *) echo "$0: invalid option: $1" >&2; echo ${usage} >&2; exit 1;;
376  esac
377  shift
378 done
379 if $clean; then
380  tmp="/tmp/tmp$$"
381 else
382  tmp="output"
383 fi
384 if test "${surface_only}" = true; then
385  dim=`expr ${map_dim} + 1`
386 else
387  dim=${map_dim}
388 fi
389 #echo "dim=${dim}" 1>&2
390 #echo "map_dim=${map_dim}" 1>&2
391 #echo "variant=$variant" 1>&2
392 #echo "n=$n" 1>&2
393 #echo "a=$a" 1>&2
394 #echo "b=$b" 1>&2
395 mkgmsh ${dim} ${surface_only} $variant $n $a $b $c $d $f $g > $tmp.mshcad
396 if $fix; then
397  # mkgeo_ball_gmsh_fix will set mesh order:
398  my_eval "$GMSH -${map_dim} $tmp.mshcad -format msh2 -o $tmp.msh > $tmp.log 2>&1"
399 else
400  my_eval "$GMSH -${map_dim} -order $order $tmp.mshcad -format msh2 -o $tmp.msh > $tmp.log 2>&1"
401 fi
402 if test ! -f $tmp.msh; then
403  echo "$0: gmsh failed"
404  exit 1
405 fi
406 MSH2GEO="${bindir}msh2geo"
407 GEO="${bindir}geo"
408 my_eval "$MSH2GEO $tmp.msh > ${tmp}-v2.geo"
409 # upgrade: no more need with msh2geo
410 # my_eval "$GEO -upgrade -geo -check - < ${tmp}-v1.geo > ${tmp}-v2.geo"
411 if $fix; then
412  fix_filter="$pkgbindir/mkgeo_ball_gmsh_fix -order $order"
413 else
414  fix_filter="cat"
415 fi
416 my_eval "$fix_filter < ${tmp}-v2.geo"
417 if $clean; then
418  my_eval "rm -f $tmp.log ${tmp}-v1.geo ${tmp}-v2.geo"
419  my_eval "rm -f $tmp.mshcad $tmp.msh"
420 fi
421