Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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# ------------------------------------------
133verbose=false
134
135my_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# ------------------------------------------
147mkgmsh_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
156cat << EOF_2D_1
157Mesh.SecondOrderIncomplete = 0;
158/* Mesh.ElementOrder = $order; */
159/*
160Mesh.SecondOrderLinear = 0;
161Mesh.SecondOrderExperimental = 0;
162Mesh.SmoothInternalEdges = 100;
163Mesh.Smoothing = 100;
164*/
165n = $n; // the density of discretisation
166a = $a; b = $b;
167c = $c; d = $d;
168EOF_2D_1
169if test $variant = t -o $variant = tq; then
170 echo "h = 1.0*(b-a)/n;"
171else
172 echo "h = 2.0*(b-a)/n;"
173fi
174cat << EOF_2D_2
175ab = (a+b)/2;
176cd = (c+d)/2;
177Point(40) = {ab, cd, 0, h};
178Point(41) = { b, cd, 0, h};
179Point(42) = {ab, d, 0, h};
180Point(43) = { a, cd, 0, h};
181Point(44) = {ab, c, 0, h};
182Circle(21) = {41, 40, 42};
183Circle(22) = {42, 40, 43};
184Circle(23) = {43, 40, 44};
185Circle(24) = {44, 40, 41};
186Line Loop (51) = {21,22,23,24};
187Physical Line("boundary") = {21,22,23,24};
188EOF_2D_2
189if 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
207fi
208}
209# ------------------------------------------
210# 3d case
211# http://sites.google.com/site/auxcapucins/maillage-3d-en-gmsh---maillage-d-une-sphere
212# ------------------------------------------
213mkgmsh_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
224cat << EOF_3D_1
225Mesh.SecondOrderIncomplete = 0;
226n = $n; // the density of discretisation
227a = $a; b = $b;
228c = $c; d = $d;
229f = $f; g = $g;
230ab = (a+b)/2;
231cd = (c+d)/2;
232fg = (f+g)/2;
233EOF_3D_1
234if test $variant = t -o $variant = tq; then
235 echo "h = 1.0*(b-a)/n;"
236else
237 echo "h = 2.0*(b-a)/n;"
238fi
239cat << EOF_3D_2
240Point(1) = {ab, cd, fg, h};
241Point(2) = { b, cd, fg, h};
242Point(3) = {ab, d, fg, h};
243Point(4) = {ab, cd, g, h};
244Point(5) = {a, cd, fg, h};
245Point(6) = {ab, c, fg, h};
246Point(7) = {ab, cd, f, h};
247Circle(1) = {2,1,3}; // {start,center,end}
248Circle(2) = {3,1,5};
249Circle(3) = {5,1,6};
250Circle(4) = {6,1,2};
251Circle(5) = {2,1,7};
252Circle(6) = {7,1,5};
253Circle(7) = {5,1,4};
254Circle(8) = {4,1,2};
255Circle(9) = {6,1,7};
256Circle(10) = {7,1,3};
257Circle(11) = {3,1,4};
258Circle(12) = {4,1,6};
259Line Loop(1) = {1,11,8};
260Line Loop(2) = {2,7,-11};
261Line Loop(3) = {3,-12,-7};
262Line Loop(4) = {4,-8,12};
263Line Loop(5) = {5,10,-1};
264Line Loop(6) = {-2,-10,6};
265Line Loop(7) = {-3,-6,-9};
266Line Loop(8) = {-4,9,-5};
267Mesh.Algorithm = 1;
268Ruled Surface(1) = {1};
269Ruled Surface(2) = {2};
270Ruled Surface(3) = {3};
271Ruled Surface(4) = {4};
272Ruled Surface(5) = {5};
273Ruled Surface(6) = {6};
274Ruled Surface(7) = {7};
275Ruled Surface(8) = {8};
276Surface Loop (1) = {1,2,3,4,5,6,7,8};
277Physical Surface("boundary") = {1,2,3,4,5,6,7,8};
278EOF_3D_2
279if 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
290else
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};"
301fi
302}
303#----------------------------------------------
304# multi-dim switch
305#----------------------------------------------
306mkgmsh () {
307dim=$1; shift
308case $dim in
309 2) mkgmsh_2d $*;;
310 *) mkgmsh_3d $*;;
311esac
312}
313#----------------------------------------------
314# main
315#----------------------------------------------
316usage="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
327if test $# -eq 0; then
328 echo ${usage} >&2
329 exit 0
330fi
331
332GMSH=${GMSH-"gmsh"}
333pkgbindir=`rheolef-config --pkglibdir`
334verbose=false
335clean=true
336bindir=""
337map_dim=2
338variant=t
339n=10
340a=-1; b=1
341c=-1; d=1
342f=-1; g=1
343order=1
344surface_only="false"
345fix=true
346while 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
378done
379if $clean; then
380 tmp="/tmp/tmp$$"
381else
382 tmp="output"
383fi
384if test "${surface_only}" = true; then
385 dim=`expr ${map_dim} + 1`
386else
387 dim=${map_dim}
388fi
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
395mkgmsh ${dim} ${surface_only} $variant $n $a $b $c $d $f $g > $tmp.mshcad
396if $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"
399else
400 my_eval "$GMSH -${map_dim} -order $order $tmp.mshcad -format msh2 -o $tmp.msh > $tmp.log 2>&1"
401fi
402if test ! -f $tmp.msh; then
403 echo "$0: gmsh failed"
404 exit 1
405fi
406MSH2GEO="${bindir}msh2geo"
407GEO="${bindir}geo"
408my_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"
411if $fix; then
412 fix_filter="$pkgbindir/mkgeo_ball_gmsh_fix -order $order"
413else
414 fix_filter="cat"
415fi
416my_eval "$fix_filter < ${tmp}-v2.geo"
417if $clean; then
418 my_eval "rm -f $tmp.log ${tmp}-v1.geo ${tmp}-v2.geo"
419 my_eval "rm -f $tmp.mshcad $tmp.msh"
420fi
421