#!/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
#
# -------------------------------------------------------------------------
# author: Pierre.Saramito@imag.fr
# date: 1 nov 2011

##
#@commandfile mkgeo_ugrid unstructured mesh of a parallelepiped
#@addindex command `mkgeo_ugrid`
#@addindex command `geo`
#@addindex command `gmsh`
#@addindex mesh
#@addindex file format `.geo` mesh
#
#Synopsis
#---------
#
#    mkgeo_ugrid [options] n
#
#Example
#-------
#
#        mkgeo_ugrid -t 10 > square-10.geo
#        geo square-10.geo
#
#The previous command build a triangle based 2d unstructured mesh
#of the unit square:
#
#Description
#-----------
#This command is useful when testing programs on simple geometries.
#Invocation is similar to those of @ref mkgeo_grid_1 .
#It callsg `gmsh` as unstructured mesh generator.
#It avoid the preparation of an input file for a mesh generator.
#The optional *n* argument is an integer
#that specifies the subdivision in each direction. By default 
#n=10.
#The mesh files goes on standard output.
#
#The command supports all the possible element types: edges, triangles,
#rectangles, tetrahedra, prisms and hexahedra. It supports also mixed 2D 
#with triangles and quadrangles:  
#
#        mkgeo_ugrid -tq 10 | geo -
#
#and mixed 3D with tetraedra, prisms and/or hexahedra:
#
#        mkgeo_ugrid -TP  10 | geo -paraview -
#        mkgeo_ugrid -PH  10 | geo -paraview -
#        mkgeo_ugrid -TPH 10 | geo -paraview -
#
#Element type option
#-------------------
#
#`-e`
#>	1d mesh using edges.
#`-t`
#>	2d mesh using triangles.
#`-q`
#>	2d mesh using quadrangles (rectangles).
#`-T`
#>	3d mesh using tetrahedra.
#`-P`
#>	3d mesh using prisms.
#`-H`
#>	3d mesh using hexahedra.
#
#`-tq`
#>	2d mesh using both triangles and quadrangles.
#`-TP` \n
#`-PH` \n
#`-TPH`
#>	3d mixed mesh combining tetrahedra, prisms and/or hexahedra.
#
#The geometry
#------------
#`-a` *float* \n
#`-b` *float* \n
#`-c` *float* \n
#`-d` *float* \n
#`-f` *float* \n
#`-g` *float*
#>	The geometry can be any [a,b] segment, [a,b]x[c,d] rectangle
#>	or [a,b]x[c,d]x[f,g] parallelepiped. By default a=c=f=0 and b=d=g=1, thus,
#>	the unit boxes are considered.
#
#Boundary domains
#----------------
#`-[no]sides`
#>	The boundary sides are represented by domains: `left`, `right`,
#>	`top`, `bottom`, `front` and `back`.
#`-[no]boundary`
#>	This option defines a domain named `boundary` that groups all sides.
#
#Regions
#-------
#`-[no]region`
#>	The whole domain is split into two subdomains: `east` and `west`.
#>	Also, the separating domain is named `interface` in the mesh.
#>	This option is used for testing computations with subdomains (e.g. transmission
#>	problem; see @ref usersguide_page ).
#
#Corners
#-------
#`-[no]corner`
#>	The corners (four in 2D and eight in 3D) are defined as OD-domains.
#>	This could be useful for some special boundary conditions.
#
#Coordinate system options
#-------------------------
#@addindex axisymmetric coordinate system
#
#`-cartesian` \n
#`-rz` \n
#`-zr`
#>       Specifies the coordinate system.
#>       The default is `cartesian` while `-rz`
#>       and `-zr` denotes some axisymmetric coordinate systems.
#>       Recall that most of Rheolef codes are coordinate-system independent
#>       and the coordinate system is specified in the geometry file `.geo`.
#
#@addindex mesh order
#The mesh order
#--------------
#This option is related to curved boundaries.
#Since boundaries are here flat, this option has no practical effect
#and is provided for test purpose only.
#
#`-order` *int*
#>      The polynomial approximation mesh order, as defined by `gmsh`.
#>      This option enable a possible curved boundary, when applying
#>	a suitable nonlinear transformation to the mesh.
#>       Default is order=1.
#
#Implementation
#--------------
#@showfromfile

# ------------------------------------------
# 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]
	[-[no]clean]
	[-[no]verbose]
"

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;;
  -verbose)    verbose=true;;
  -noverbose)  verbose=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 -format msh2 $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"
# upgrade: no more need with msh2geo
#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"
fi

