rheolef 6.5-1+b1 / usr / bin / mkgeo_ugrid

#!/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_ugrid
#NAME: @code{mkgeo_ugrid} -- build an unstructured mesh of a parallelotope, in 1d, 2d or 3d
#@pindex mkgeo_ugrid
#@pindex geo
#@cindex mesh
#@fiindex @file{.geo} mesh
#SYNOPSIS:
#@example
#	mkgeo_ugrid @var{options} [@var{n}]
#@end example
#EXAMPLE: 
#@noindent
#  The following command build a triangle based 2d unstructured mesh
#  of the unit square:
#@example
#	mkgeo_ugrid -t 10 > square-10.geo
#	geo -mayavi square-10.geo
#@end example
#@noindent
#  or in one comand line:
#@example
#	mkgeo_ugrid -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}).
#  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.
#
#  @noindent
#  The command supports all the possible element types: edges, triangles,
#  rectangles, tetraedra, prisms and hexahedra. It supports also mixed 2D 
#  with triangles and quadrangles:  
#@example
#	mkgeo_ugrid -tq 10 | geo -mayavi -
#@end example
#  and mixed 3D with tetraedra, prisms and/or hjexaedra:
#@example
#	mkgeo_ugrid -TP  10 | geo -mayavi -
#	mkgeo_ugrid -PH  10 | geo -mayavi -
#	mkgeo_ugrid -TPH 10 | geo -mayavi -
#@end example
#ELEMENT TYPE OPTIONS:
#@table @code
#@item -e
#	1d mesh using edges.
#@item -t
#	2d mesh using triangles.
#@item -q
#	2d mesh using quadrangles.
#@item -tq
#	2d mesh using both triangles and quadrangles.
#@item -T
#	3d mesh using tetraedra.
#@item -P
#	3d mesh using prisms.
#@item -H
#	3d mesh using hexahedra.
#@item -TP
#@itemx -PH
#@itemx -TPH
#	3d mesh using a mixt between tetraedra, prisms and/or hexahedra.
#@end table
#THE GEOMETRY:
#  @noindent
#  The geometry can be any [a,b] segment, [a,b]x[c,d] rectangle
#  or [a,b]x[c,d]x[f,g] parallelotope. By default a=c=f=0 and b=d=g=1, thus,
#  the unit boxes are considered. For instance, the following
#  command meshes the [-2,2]x[-1.5, 1.5] rectangle:
#@example
#	mkgeo_ugrid -t 10 -a -2 -b 2 -c -1.5 -d 1.5 | geo -
#@end example
#@table @code
#@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:
#@table @code
#@item -sides
#@itemx -nosides
#  The boundary sides are representd by domains: @code{left}, @code{right},
#  @code{top}, @code{bottom},@code{front} and @code{back}.
#@item -boundary
#@itemx -noboundary
#  This option defines a domain named @code{boundary} that groups all sides.
#@end table
#By default, both sides and the whole boundary are defined as domains:
#@example
#	mkgeo_ugrid -t 10 > square.geo 
#       geo square.geo
#	mkgeo_ugrid -t 10 -nosides > square.geo 
#       geo square.geo
#	mkgeo_ugrid -t 10 -noboundary > square.geo 
#       geo square.geo
#	mkgeo_ugrid -t 10 -noboundary -nosides > square.geo 
#       geo square.geo
#@end example
#REGIONS:
#@table @code
#@item -region
#@itemx -noregion
#  The whole domain is splitted into two subdomains: @code{east} and @code{west},
#  This option is used for testing computations with subdomains (e.g. transmission
#  problem; see the user manual).
#@end table
#@example
#	mkgeo_ugrid -t 10 -region | geo -
#@end example
#CORNERS:
#@table @code
#@item -corner
#@itemx -nocorner
#  The corners (four in 2D and eight in 3D) are defined as OD-domains.
#  This could be usefull for some special boundary conditions.
#@end table
#@example
#	mkgeo_ugrid -t 10 -corner | geo -
#	mkgeo_ugrid -T  5 -corner | geo -
#@end example
#
#THE MESH ORDER:
#@pindex gmsh
#@cindex mesh order
#@table @code
#@item -order int
#       The polynomial approximation mesh order, as defined by @code{gmsh}.
#       This option enable a possible curved boundary, when applying
#	a suitable nonlinear transformation to the mesh.
#       Defualt is order=1.
#@end table
#OTHERS OPTIONS:
#@table @code
#  @item -clean
#        clear temporary files (this is the default).
#  @item -noclean
#        does not clear temporary files.
#@end table
#
#DATE:
#    1 nov 2011
#END:

# ------------------------------------------
# utility
# ------------------------------------------
verbose=false
my_eval () {
  command="$*"
  if $verbose; then echo "! $command" 1>&2; fi
  eval $command
  if test $? -ne 0; then
    echo "$0: error on command: $command" >&2
    exit 1
  fi
}
# ------------------------------------------
# 1d case
# ------------------------------------------
mkgmsh_1d () {
  variant=$1; shift
  n=$1; shift
  sides=$1; shift;
  boundary=$1; shift;
  corner=$1; shift;
  a=$1
  b=$2

cat << EOF_1D
n = $n;
a = $a;
b = $b;
h = 1.0*(b-a)/n;
Point(1) = {a, 0, 0, h};
Point(2) = {b, 0, 0, h};
Line(3) = {1,2};
EOF_1D
if $sides || $corner; then
  echo "Physical Point(\"left\")  = {1};"
  echo "Physical Point(\"right\") = {2};"
fi
if $boundary; then
  echo "Physical Point(\"boundary\") = {1,2};"
fi
echo "Physical Line(\"interior\") = {3};"
}
# ------------------------------------------
# 2d case: t,q,tq variants
# ------------------------------------------
mkgmsh_2d () {
  variant=$1; shift
  n=$1; shift
  sides=$1; shift;
  boundary=$1; shift;
  corner=$1; shift;
  a=$1
  b=$2
  c=$3
  d=$4

echo "n = $n;"
echo "a = $a; b = $b;"
echo "c = $c; d = $d;"
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_1
Point(1) = {a, c, 0, h};
Point(2) = {b, c, 0, h};
Point(3) = {b, d, 0, h};
Point(4) = {a, d, 0, h};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line Loop(5) = {1,2,3,4};
Mesh.Algorithm = 1;
Plane Surface(6) = {5} ;
EOF_2D_1
if test $variant = q; then
  echo "Mesh.RecombinationAlgorithm = 1;"
  echo "Mesh.SubdivisionAlgorithm = 1;"
  echo "Recombine Surface {6};"
fi
if test $variant = tq; then
  echo "Mesh.RecombinationAlgorithm = 0;"
  echo "angle = 90.0;"
  echo "Recombine Surface {6} = angle;"
fi
if $corner; then
  echo "Physical Point(\"left_bottom\")  = {1};"
  echo "Physical Point(\"right_bottom\") = {2};"
  echo "Physical Point(\"right_top\")    = {3};"
  echo "Physical Point(\"left_top\")     = {4};"
fi
if $boundary; then
  echo "Physical Line(\"boundary\") = {1,2,3,4};"
fi
if $sides; then
  echo "Physical Line(\"bottom\") = {1};"
  echo "Physical Line(\"right\")  = {2};"
  echo "Physical Line(\"top\")    = {3};"
  echo "Physical Line(\"left\")   = {4};"
fi
echo "Physical Surface(\"interior\") = {6};"
}
# ------------------------------------------
# 3d case: T,TP,TPH variants
# ------------------------------------------
mkgmsh_3d_T_TP_TPH () {
  variant=$1; shift
  n=$1; shift
  sides=$1; shift;
  boundary=$1; shift;
  corner=$1; shift;
  a=$1
  b=$2
  c=$3
  d=$4
  f=$5
  g=$6

cat << EOF_3D_T_1
n = $n;
a = $a; c = $c; f = $f;
b = $b; d = $d; g = $g;
EOF_3D_T_1
if test $variant = T; then
  d=" d"
  g=" g"
elif test $variant = TP; then
  echo "fg = (f+g)/2.0;"
  d=" d"
  g="fg"
elif test $variant = TPH; then
  echo "cd = (c+d)/2.0;"
  echo "fg = (f+g)/2.0;"
  d="cd"
  g="fg"
fi
cat << EOF_3D_T_2
h = (b-a)/n;
Point(1) = {a,  c,  f, h};
Point(2) = {b,  c,  f, h};
Point(3) = {b, $d,  f, h};
Point(4) = {a, $d,  f, h};
Point(5) = {a,  c, $g, h};
Point(6) = {b,  c, $g, h};
Point(7) = {b, $d, $g, h};
Point(8) = {a, $d, $g, h};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};
Line(9)  = {1,5};
Line(10) = {2,6};
Line(11) = {3,7};
Line(12) = {4,8};
Line Loop(21) = {-1,-4,-3,-2};
Line Loop(22) = {5,6,7,8};
Line Loop(23) = {1,10,-5,-9};
Line Loop(25) = {2,11,-6,-10};
Line Loop(24) = {12,-7,-11,3};
Line Loop(26) = {9,-8,-12,4};
Plane Surface(31) = {21} ;
Plane Surface(32) = {22} ;
Plane Surface(33) = {23} ;
Plane Surface(34) = {24} ;
Plane Surface(35) = {25} ;
Plane Surface(36) = {26} ;
Surface Loop(41) = {31,32,33,34,35,36};
Volume(51) = {41};
EOF_3D_T_2
if test $variant = TP; then
   echo "extr[] = Extrude{0,0,(g-f)/2}{ Surface{32}; Layers{n/2}; Recombine;};"
elif test $variant = TPH; then
   echo "extr_y[] = Extrude{0,(d-c)/2,0}{ Surface{34}; Layers{n/2}; Recombine;};"
   echo "extr_z[] = Extrude{0,0,(g-f)/2}{ Surface{32,extr_y[3]}; Layers{n/2}; Recombine;};"
fi
if test $variant = T; then
  if $boundary; then
    echo "Physical Surface(\"boundary\") = {31, 32, 33, 34, 35, 36};"
  fi
  if $sides; then
    echo "Physical Surface(\"bottom\") = {31};"
    echo "Physical Surface(\"top\")    = {32};"
    echo "Physical Surface(\"left\")   = {33};"
    echo "Physical Surface(\"front\")  = {35};"
    echo "Physical Surface(\"right\")  = {34};"
    echo "Physical Surface(\"back\")   = {36};"
 fi
 echo "Physical Volume(\"internal\") = {51};"
elif test $variant = TP; then
  if $boundary; then
    echo "Physical Surface(\"boundary\") = {31, extr[0], 33, extr[2], 35, extr[3], 34, extr[4], 36, extr[5]};"
  fi
  if $sides; then
    echo "Physical Surface(\"bottom\") = {31};"
    echo "Physical Surface(\"top\")   = {extr[0]};"
    echo "Physical Surface(\"left\")  = {33, extr[2]};"
    echo "Physical Surface(\"front\") = {35, extr[3]};"
    echo "Physical Surface(\"right\") = {34, extr[4]};"
    echo "Physical Surface(\"back\")  = {36, extr[5]};"
  fi
  echo "Physical Volume(\"internal\") = {51,extr[1]};"
elif test $variant = TPH; then
  if $boundary; then
    echo "Physical Surface(\"boundary\") = {31, extr_y[5], extr_z[0], extr_z[6],"
    echo "				33,extr_z[2], 35, extr_y[4], extr_z[3], extr_z[9],"
    echo "				extr_y[0],extr_z[10], 36, extr_y[2], extr_z[5], extr_z[11]};"
  fi
  if $sides; then
    echo "Physical Surface(\"bottom\") = {31, extr_y[5]};"
    echo "Physical Surface(\"top\")   = {extr_z[0], extr_z[6]};"
    echo "Physical Surface(\"left\")  = {33,extr_z[2]};"
    echo "Physical Surface(\"front\") = {35, extr_y[4], extr_z[3], extr_z[9]};"
    echo "Physical Surface(\"right\")  = {extr_y[0],extr_z[10]};"
    echo "Physical Surface(\"back\")  = {36, extr_y[2], extr_z[5], extr_z[11]};"
  fi
  echo "Physical Volume(\"internal\") = {51,extr_y[1],extr_z[1],extr_z[7]};"
fi
}
# ------------------------------------------
# 3d case: P,H,PH variants
# ------------------------------------------
mkgmsh_3d_P_H_PH () {
  variant=$1; shift
  n=$1; shift
  sides=$1; shift;
  boundary=$1; shift;
  corner=$1; shift;
  a=$1
  b=$2
  c=$3
  d=$4
  f=$5
  g=$6

if test $variant = H -a $n -ne 1; then
  n="$n/2";
fi
cat << EOF_3D_PH
n = $n;
a = $a; c = $c; f = $f;
b = $b; d = $d; g = $g;
h = 1.0*(b-a)/n;
Point(1) = {a, c, f, h};
Point(2) = {b, c, f, h};
Point(6) = {b, c, g, h};
Point(5) = {a, c, g, h};
Line(1) = {1,2};
Line(6) = {2,6};
Line(9) = {5,6};
Line(5) = {1,5};
Line Loop(300) = {1,6,-9,-5};
Mesh.Algorithm = 1;
EOF_3D_PH
if test $variant = P; then
  echo "Plane Surface(3) = {300};"
elif test $variant = H; then
  echo "Plane Surface(3) = {300};"
  echo "Mesh.SubdivisionAlgorithm = 1;"
  echo "Mesh.RecombinationAlgorithm = 1;"
  echo "Recombine Surface {3};"
  echo "Mesh.SubdivisionAlgorithm = 2;"
elif test $variant = PH; then
  echo "Plane Surface(3) = {300};"
  echo "Mesh.RecombinationAlgorithm = 0;"
  echo "angle = 45.0;"
  echo "Recombine Surface {3} = angle;"
fi
echo "extr[] = Extrude{0,d-c,0}{ Surface{3}; Layers{n}; Recombine;};"
if $boundary; then
  echo "Physical Surface(\"boundary\") = {-extr[2], -extr[4], 3, -extr[3], -extr[0], -extr[5]};"
fi
if $sides; then
  echo "Physical Surface(\"bottom\")  = {-extr[2]};"
  echo "Physical Surface(\"top\")     = {-extr[4]};"
  echo "Physical Surface(\"left\")    = {3};"
  echo "Physical Surface(\"front\")   = {-extr[3]};"
  echo "Physical Surface(\"right\")   = {-extr[0]};"
  echo "Physical Surface(\"back\")    = {-extr[5]};"
fi
echo "Physical Volume(\"interior\") = {extr[1]};"
}
#----------------------------------------------
# multi-dim switch
#----------------------------------------------
mkgmsh () {
dim=$1; shift
case $dim in
 1) mkgmsh_1d $*;;
 2) mkgmsh_2d $*;;
 *) variant=$1
    case $variant in
    T*) mkgmsh_3d_T_TP_TPH  $*;;
    *)  mkgmsh_3d_P_H_PH    $*;;
    esac
esac
}
#----------------------------------------------
# main
#----------------------------------------------
usage="mkgeo_ugrid
	[-{eptqTPH}|-tq|-TP|-PH|-TPH]
	[n]
	[-{abcdfg} float]
	[-[no]sides]
	[-[no]boundary]
	[-[no]corner]
	[-order int]
"

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

clean=true
bindir=""
dim=2
variant=t
n=10
a=0; b=1
c=0; d=1
f=0; g=1
boundary=true
sides=true
corner=false
order=1
while test $# -ne 0; do
  case $1 in
  -h) echo ${usage} >&2; exit 0;;
  -e)                  dim=1; variant=`echo x$1 | sed -e 's/x-//'`;;
  -[tq]|-tq)           dim=2; variant=`echo x$1 | sed -e 's/x-//'`;;
  -[TPH]|-TP|-PH|-TPH) 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 ${usage} >&2; exit 1; fi
	expr="$var=$2"
	eval $expr
	shift
	;;
  [0-9]*) n=$1;;
  -order)      order=$2; shift;;
  -sides)      sides=true;;
  -nosides)    sides=false;;
  -boundary)   boundary=true;;
  -noboundary) boundary=false;;
  -corner)     corner=true;;
  -nocorner)   corner=false;;
  -clean)      clean=true;;
  -noclean)    clean=false;;
  -bindir)    
        if test x"$2" = x""; then echo ${usage} >&2; exit 1; fi
	bindir="$2/"
	shift
	;;
  *) echo ${usage} >&2; exit 1;;
  esac
  shift
done
if $clean; then
  tmp="/tmp/tmp$$"
else
  tmp="output"
fi
#echo "dim=$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 $variant $n $sides $boundary $corner $a $b $c $d $f $g > $tmp.mshcad
my_eval "gmsh -$dim -order $order $tmp.mshcad -o $tmp.msh > $tmp.log"
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 - < ${tmp}-v1.geo"
if $clean; then
  my_eval "rm -f $tmp.mshcad $tmp.log $tmp.msh ${tmp}-v1.geo"
fi