rheolef 6.7-1+b4 / usr / bin / mkgeo_ball

#!/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
#
# -------------------------------------------------------------------------

#Prog:mkgeo_ball
#NAME: @code{mkgeo_ball} -- build an unstructured mesh of an ellipsoid, in 2d or 3d
#@pindex mkgeo_ball
#@pindex geo
#@cindex mesh
#@fiindex @file{.geo} mesh
#SYNOPSIS:
#@example
#	mkgeo_ball @var{options} [@var{n}]
#@end example
#EXAMPLE: 
#@noindent
#  The following command build a triangle based 2d unstructured mesh
#  of the unit ball
#@example
#	mkgeo_ball -t 10 > ball-10.geo
#	geo -mayavi ball-10.geo
#@end example
#@noindent
#  or in one comand line:
#@example
#	mkgeo_ball -t 10 | geo -mayavi -
#@end example
#DESCRIPTION:       
#@pindex gmsh
#  @noindent
#  This command is usefull when testing programs on simple geometries.
#  Invocation is similar to @code{mkgeo_grid}
#  (see @ref{mkgeo_grid command}, @ref{mkgeo_ugrid command}).
#  It calls @code{gmsh} as unstructured mesh generator.
#  It avoid the preparation of an input file for a mesh generator.
#  The optional @var{n} argument is an integer
#  that specifies the subdivision in each direction. By default 
#  @var{n}=10.
#  The mesh files goes on standard output.
#ELEMENT TYPE OPTIONS:
#@table @code
#@item -t
#	2d mesh using triangles.
#@item -q
#	2d mesh using quadrangles.
#@item -tq
#	2d mesh using both triangles and quadrangles.
#@end table
#THE MESH ORDER:
#@pindex gmsh
#@cindex mesh order
#@table @code
#@item -order int
#	The polynomial pproximation mesh order,
#       for the curved boundary, as defined by @code{gmsh}.
#       Defualt is order=1.
#@end table
#THE GEOMETRY:
#  @noindent
#  The geometry can be ellipse/ellipsoid inscribted in the [a,b]x[c,d] rectangle
#  or the [a,b]x[c,d]x[f,g] parallelotope. By default a=c=f=-1 and b=d=g=1, thus,
#  the unit balls are considered.
#@table @code
#@item -s
#	Only a surface mesh is generated: a curved line in 2d or a curved surface mesh in 3d.
#	In 3d, supports both '-t' and '-q' options.
#@item -a @var{float}
#@itemx -b @var{float}
#@itemx -c @var{float}
#@itemx -d @var{float}
#@itemx -f @var{float}
#@itemx -g @var{float}
#@end table
#BOUNDARY DOMAINS:
#  The meshes defines a domain named @code{boundary} that groups all boundary sides.
#OTHERS OPTIONS:
#@table @code
#  @item -fix
#  @itemx -nofix
#        By default, internal face and volume node from gmsh are recomputed,
#        since they have incorrect coordinate that destroy the convergence 
#        properties of isoparametric high order elements (order >= 3).
#        The @code{-nofix} option skip this correction: this option is
#        available for test purpose.
#  @item -clean
#        clear temporary files (this is the default).
#  @item -noclean
#        does not clear temporary files.
#  @item -verbose
#  @item -noverbose
#        print intermediate commands and information messages.
#@end table
#
#DATE:
#    1 nov 2011
#END:

# ------------------------------------------
# utility
# ------------------------------------------
verbose=false

my_eval () {
  command="$*"
  if test "$verbose" = true; then echo "! $command" 1>&2; fi
  eval $command
  if test $? -ne 0; then
    echo "$0: error on command: $command" >&2
    exit 1
  fi
}
# ------------------------------------------
# 2d case: t,q,tq variants
# ------------------------------------------
mkgmsh_2d () {
  surface_only=$1; shift
  variant=$1; shift
  n=$1; shift
  a=$1
  b=$2
  c=$3
  d=$4

cat << EOF_2D_1
Mesh.SecondOrderIncomplete = 0;
/* Mesh.ElementOrder = $order; */
/*
Mesh.SecondOrderLinear = 0;
Mesh.SecondOrderExperimental = 0;
Mesh.SmoothInternalEdges = 100;
Mesh.Smoothing = 100;
*/
n = $n; // the density of discretisation
a = $a;  b = $b;
c = $c;  d = $d;
EOF_2D_1
if test $variant = t -o $variant = tq; then
  echo "h = 1.0*(b-a)/n;"
else
  echo "h = 2.0*(b-a)/n;"
fi
cat << EOF_2D_2
ab = (a+b)/2;
cd = (c+d)/2;
Point(40) = {ab, cd, 0, h};
Point(41) = { b, cd, 0, h};
Point(42) = {ab,  d, 0, h};
Point(43) = { a, cd, 0, h};
Point(44) = {ab,  c, 0, h};
Circle(21) = {41, 40, 42};
Circle(22) = {42, 40, 43};
Circle(23) = {43, 40, 44};
Circle(24) = {44, 40, 41};
Line Loop (51) = {21,22,23,24};
Physical Line("boundary") = {21,22,23,24};
EOF_2D_2
if test "${surface_only}" != true; then
  echo "Mesh.Algorithm = 1;"
  echo "Plane Surface(61) = {51};"
  echo "Physical Surface(\"interior\") = {61};"
  if test $variant = q; then
    echo "Mesh.RecombinationAlgorithm = 1;" ; 
    # Note:  Mesh.RecombinationAlgorithm = 0 : standard => triangle with 2 bdry edges
    #        Mesh.RecombinationAlgorithm = 1 : blossom => triangle with 2 bdry edges, too
    #  => mkgeo_ball_gmsh_fix fails
    #  => try a higher mesh density...  sometimes its ok...
    echo "Mesh.SubdivisionAlgorithm = 1;"
    echo "Recombine Surface {61};"
  fi
  if test $variant = tq; then
    echo "Mesh.RecombinationAlgorithm = 0;"
    echo "angle = 90.0;"
    echo "Recombine Surface {61} = angle;"
  fi
fi
}
# ------------------------------------------
# 3d case
# http://sites.google.com/site/auxcapucins/maillage-3d-en-gmsh---maillage-d-une-sphere
# ------------------------------------------
mkgmsh_3d () {
  surface_only=$1; shift
  variant=$1; shift
  n=$1; shift
  a=$1
  b=$2
  c=$3
  d=$4
  f=$5
  g=$6

cat << EOF_3D_1
Mesh.SecondOrderIncomplete = 0;
n = $n; // the density of discretisation
a = $a;  b = $b;
c = $c;  d = $d;
f = $f;  g = $g;
ab = (a+b)/2;
cd = (c+d)/2;
fg = (f+g)/2;
EOF_3D_1
if test $variant = t -o $variant = tq; then
  echo "h = 1.0*(b-a)/n;"
else
  echo "h = 2.0*(b-a)/n;"
fi
cat << EOF_3D_2
Point(1) = {ab, cd, fg, h};
Point(2) = { b, cd, fg, h};
Point(3) = {ab,  d, fg, h};
Point(4) = {ab, cd,  g, h};
Point(5) = {a,  cd, fg, h};
Point(6) = {ab, c,  fg, h};
Point(7) = {ab, cd, f,  h};
Circle(1) = {2,1,3}; // {start,center,end}
Circle(2) = {3,1,5};
Circle(3) = {5,1,6};
Circle(4) = {6,1,2};
Circle(5) = {2,1,7};
Circle(6) = {7,1,5};
Circle(7) = {5,1,4};
Circle(8) = {4,1,2};
Circle(9) = {6,1,7};
Circle(10) = {7,1,3};
Circle(11) = {3,1,4};
Circle(12) = {4,1,6};
Line Loop(1) = {1,11,8};
Line Loop(2) = {2,7,-11};
Line Loop(3) = {3,-12,-7};
Line Loop(4) = {4,-8,12};
Line Loop(5) = {5,10,-1};
Line Loop(6) = {-2,-10,6};
Line Loop(7) = {-3,-6,-9};
Line Loop(8) = {-4,9,-5};
Mesh.Algorithm = 1;
Ruled Surface(1) = {1};
Ruled Surface(2) = {2};
Ruled Surface(3) = {3};
Ruled Surface(4) = {4};
Ruled Surface(5) = {5};
Ruled Surface(6) = {6};
Ruled Surface(7) = {7};
Ruled Surface(8) = {8};
Surface Loop (1) = {1,2,3,4,5,6,7,8};
Physical Surface("boundary") = {1,2,3,4,5,6,7,8};
EOF_3D_2
if test "${surface_only}" = true; then
  if test $variant = q; then
    echo "Mesh.RecombinationAlgorithm = 1;"
    echo "Mesh.SubdivisionAlgorithm = 1;"
    echo "Recombine Surface {1};"
  fi
  if test $variant = tq; then
    echo "Mesh.RecombinationAlgorithm = 0;"
    echo "angle = 90.0;"
    echo "Recombine Surface {1} = angle;"
  fi
else
  echo "Mesh.Algorithm3D = 1; // attention: Tetgen is not available in debian packages"
  echo "Mesh.OptimizeNetgen = 1;"
  echo "Mesh.Optimize = 1;"
  echo "Mesh.Smoothing = 10;"
  echo "Volume (1) = {1};"
  if test $variant = H; then
    echo "Mesh.RecombinationAlgorithm = 1;"
    echo "Mesh.SubdivisionAlgorithm = 2;"
  fi
  echo "Physical Volume(\"internal\") = {1};"
fi
}
#----------------------------------------------
# multi-dim switch
#----------------------------------------------
mkgmsh () {
dim=$1; shift
case $dim in
 2) mkgmsh_2d $*;;
 *) mkgmsh_3d $*;;
esac
}
#----------------------------------------------
# main
#----------------------------------------------
usage="mkgeo_ball
	[-{tq}|-tq]
	[n]
	[-s]
	[-order k]
	[-{abcdfg} float]
	[-[no]fix]
	[-[no]clean]
	[-[no]verbose]
"

if test $# -eq 0; then
  echo ${usage} >&2
  exit 0
fi

GMSH=${GMSH-"gmsh"}
pkgbindir=`rheolef-config --pkglibdir`
verbose=false
clean=true
bindir=""
map_dim=2
variant=t
n=10
a=-1; b=1
c=-1; d=1
f=-1; g=1
order=1
surface_only="false"
fix=true
while test $# -ne 0; do
  case $1 in
  -h) echo ${usage} >&2; exit 0;;
  -clean)      clean=true;;
  -noclean)    clean=false; verbose=true;;
  -verbose)    verbose=true;;
  -noverbose)  verbose=false;;
  -fix)        fix=true;;
  -nofix)      fix=false;;
  -s)	       surface_only=true;;
  -e)                  map_dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
  -[tq]|-tq)           map_dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
  -[TPH]|-TP|-PH|-TPH) map_dim=3; variant=`echo x$1 | sed -e 's/x-//'`;;
  -[abcdfg]) var=`echo x$1 | sed -e 's/x-//'`
        if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
	expr="$var=$2"
	eval $expr
	shift
	;;
  -order) if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
	order=$2
	shift
	;;
  [0-9]*)      n=$1;;
  -bindir)    
        if test x"$2" = x""; then echo "$0: $1: missing arg" >&2; echo ${usage} >&2; exit 1; fi
	bindir="$2/"
	shift
	;;
  *) echo "$0: invalid option: $1" >&2; echo ${usage} >&2; exit 1;;
  esac
  shift
done
if $clean; then
  tmp="/tmp/tmp$$"
else
  tmp="output"
fi
if test "${surface_only}" = true; then
  dim=`expr ${map_dim} + 1`
else
  dim=${map_dim}
fi
#echo "dim=${dim}" 1>&2
#echo "map_dim=${map_dim}" 1>&2
#echo "variant=$variant" 1>&2
#echo "n=$n" 1>&2
#echo "a=$a" 1>&2
#echo "b=$b" 1>&2
mkgmsh ${dim} ${surface_only} $variant $n $a $b $c $d $f $g > $tmp.mshcad
if $fix; then
  # mkgeo_ball_gmsh_fix will set mesh order:
  my_eval "$GMSH -${map_dim}               $tmp.mshcad -o $tmp.msh > $tmp.log 2>&1"
else
  my_eval "$GMSH -${map_dim} -order $order $tmp.mshcad -o $tmp.msh > $tmp.log 2>&1"
fi
if test ! -f $tmp.msh; then
  echo "$0: gmsh failed"
  exit 1
fi
MSH2GEO="${bindir}msh2geo"
GEO="${bindir}geo"
my_eval "$MSH2GEO $tmp.msh > ${tmp}-v1.geo"
my_eval "$GEO -upgrade -geo -check - < ${tmp}-v1.geo > ${tmp}-v2.geo"
if $fix; then
  fix_filter="$pkgbindir/mkgeo_ball_gmsh_fix -order $order"
else
  fix_filter="cat"
fi
my_eval "$fix_filter < ${tmp}-v2.geo"
if $clean; then
  my_eval "rm -f $tmp.log ${tmp}-v1.geo ${tmp}-v2.geo"
  my_eval "rm -f $tmp.mshcad $tmp.msh"
fi