#!/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: 2 february 2018

##
#@commandfile mkgeo_contraction mesh of an abrupt contraction
#@addindex command `mkgeo_contraction`
#@addindex command `geo`
#@addindex command `bamg`
#@addindex mesh
#@addindex file format `.geo` mesh
#
#Synopsis
#--------
#
#    mkgeo_contraction [options] [nx [ny]]
#
#Example
#--------
#The following command build a triangle-based 2d unstructured mesh
#of the [-10,0]x[0,4] u [0,10]x[0,1] mesh
#
#        mkgeo_contraction
#        geo contraction.geo
#
#Description
#-----------
#This command is convenient for building a mesh for
#the abrupt contraction, as it is a very classical benchmark in complex
#fluid flow problems.
#It calls `bamg` unstructured mesh generator
#with anisotropy feature.
#The mesh files goes on name.geo where *name*
#is the basename for the output (see option `-name` below).
#The three auxiliary files required for automatic mesh generation
#with Rheolef combined with `bamg` are also provided:
#
#            name.bamg
#            name.bamgcad
#            name.dmn
#
#The geometry
#------------
#The geometry is [-Lu,0]x[0,d] u [0,Ld]x[0,1].
#By default c=4 is the contraction ratio and Lu=Ld=10
#are the upstream and downstream pipe lengths.
#
#`-c`  *float* \n
#`-Lu` *float* \n
#`-Ld` *float*
#>  	These options control the geometry parameters.
#
#`-cartesian`  \n
#`-zr` \n
#`-rz`
#> 	These options control the geometry coordinate system.
#>	Default is Cartesian.
#
#The discretization
#------------------
#The optional *nx* and *ny* arguments are integers
#that specifies the subdivision in each direction.
#By default nx=1 and ny=nx.
#Changing the density applies the factor 1/nx to all the mesh edges lengths.
#
#`-hmin` *float*
#>	Controls the edge length at the re-entrant corner of the contraction.
#>	Default is `hmin`=0.1.
#>	Changing the density by the previous options 
#>	applies the factor 1/n to `hmin` also.
#
#Boundary domains
#----------------
#The boundary sides are represented by domains: `axis`, `wall`, `upstream`
#and `downstream`.
#
#Others options
#--------------
#
#`-name` *string*
#>       Set the basename for the output files.
#>       By default, the basename is `contraction`.
#
#`-[no]split`
#>       Split each triangle in three subtriangles by inserting the barycenter as
#>	 an additional node. This is useful for using the Scott and Vogelius 
#>	 incompressible mixed finite element P2-P1d for velocity-pressure approximation.
#>       Default is no split.
#
#`-[no]clean`
#>	Clear temporary files (this is the default).
# 
#`-[no]verbose`
#>       In verbose mode, print to stderr all commands and logs.
#
#Implementation
#--------------
#@showfromfile

nx=""
nx_default=1
c=4
Lu=20
Ld=20
hmin=0.1
split=false
clean=true
verbose=false
name="contraction"
usage="usage: mkgeo_contraction 
	[nx=$nx_default [ny=nx]]
	[-c float=$c]
	[-Lu float=$Lu]
	[-Ld float=$Ld]
	[-rz|-zr]
        [-hmin float=$hmin]
	[-name string=$name]
	[-[no]split]
	[-[no]clean]
	[-[no]verbose]
"
pkgbindir="`rheolef-config --pkglibdir`"

while test $# -ne 0; do
  case $1 in
  [0-9]*) if test x$nx = x""; then nx=$1; else ny=$1; fi;;
  -c)   c="$2"; shift;;
  -Lu) Lu="$2"; shift;;
  -Ld) Ld="$2"; shift;;
  -zr|-rz|-cartesian) sys_coord_opt="$1";;
  -name)      name=$2; shift;;
  -hmin)      hmin=$2; shift;;
  -split)     split=true;;
  -nosplit)   split=false;;
  -clean)     clean=true;;
  -noclean)   clean=false;;
  -verbose)   verbose=true;;
  -noverbose) verbose=false;;
  -h) /bin/echo -E ${usage} >&2; exit 0;;
  *)  /bin/echo -E ${usage} >&2; exit 1;;
  esac
  shift
done
if test x"$nx" = x""; then
  nx=$nx_default
fi
if test x"$ny" = x""; then
  ny=$nx
fi

to_clean=""


# edge lengths are scaled at upstream by c and in x-dir by Lu and Ld
h=0.3
h0=`echo  $hmin   $nx | awk '{print 1.0*$1/$2}'`
hd=`echo  $h      $nx | awk '{print 1.0*$1/$2}'`
hu=`echo  $h  $c  $nx | awk '{print 1.0*$1*$2/$3}'`
hux=`echo $h  $Lu $nx | awk '{print 1.0*$1*$2/$3}'`
huy=`echo $h  $c  $ny | awk '{print 1.0*$1*$2/$3}'`
hdx=`echo $h  $Ld $nx | awk '{print 1.0*$1*$2/$3}'`
hdy=`echo $h      $ny | awk '{print 1.0*$1/$2}'`

#echo "hdx=$hdx"
#echo "hdy=$hdy"

m0=` echo $h0  | awk '{print 1./($1*$1) }'`
md=` echo $hd  | awk '{print 1./($1*$1) }'`
mu=` echo $hu  | awk '{print 1./($1*$1) }'`
mux=`echo $hux | awk '{print 1./($1*$1) }'`
muy=`echo $huy | awk '{print 1./($1*$1) }'`
mdx=`echo $hdx | awk '{print 1./($1*$1) }'`
mdy=`echo $hdy | awk '{print 1./($1*$1) }'`

c2=`echo $c $Ld | awk '{ c2 = 2*$1; print (c2 < $2 ? c2 : $2) }'`

#
# background mesh for bamg mesh generator:
# with anisotropic (1/hx^2, 1/hy^2) metric
#
#  11                  10          9
# +-------------------+-----------+
# |                   |\          |
# |                   |  \  T3    |
# |                   |    \      |
# |                   |      \    |
# |      Q1           |        \  |
# |                   |          \|8           7                6
# |                   |   Q2      +-----------+---------------+
# |1                  |2          |3   Q4     |4    Q5        | 5
# +-------------------+-----------+-----------+---------------+
# -Lu                 -c          0           c2             Ld
#
cat > $name.bamgcad << EOF2
MeshVersionFormatted
  0
Dimension
  2
Vertices
  11
 -$Lu 0  1
 -$c  0  2
 -0   0  3
  $c2 0  4
  $Ld 0  5
  $Ld 1  6
  $c2 1  7
  0   1  8
  0   $c 9
 -$c  $c 10
 -$Lu $c 11
Edges
  11
  1 2    101
  2 3    101
  3 4    101
  4 5    101
  5 6    102
  6 7    103
  7 8    103
  8 9    103
  9 10   103
 10 11   103
 11 1    104
EOF2
echo "! $name.bamgcad created" 1>&2

cat > $name.mtr << EOF1b
11 3
$mux 0 $muy
$mu  0 $mu 
$md  0 $md
$md  0 $md
$mdx 0 $mdy
$mdx 0 $mdy
$md  0 $md
$m0  0 $m0
$mu  0 $mu
$mu  0 $mu
$mux 0 $muy
EOF1b
$verbose && echo "! $name.mtr created" 1>&2
to_clean="$to_clean $name.mtr"

cat > $name.dmn << EOF3
EdgeDomainNames
  4
  axis
  downstream
  wall
  upstream
EOF3
echo "! $name.dmn created" 1>&2

# bamg -coef 1 -err  0.01 -errg 0.1 -hmin 0.004 -hmax 0.3 -ratio 0 -anisomax 1e+06 -nbv 50000 -NbJacobi 1 -CutOff 1e-07 -splitpbedge  -RelError -b sector-10-P2-Bi-0.5-n-1-v2.bamg -Mbb sector-10-P2-Bi-0.5-n-1-v2-crit.bb -o sector-10-P2-Bi-0.5-n-1-v2-001.bamg 1>&2

command="bamg -g $name.bamgcad -M $name.mtr -o $name.bamg"
if $verbose; then
  command="$command 1>&2"
else
  command="$command 1> $name.bamglog"
  to_clean="$to_clean $name.bamglog"
fi
$verbose && echo "! $command" 1>&2
eval    $command
status=$?
if test $status -ne 0; then
  if $verbose; then true; else cat $name.bamglog 1>&2; fi
  echo "$0: command failed" 1>&2
  exit $status
fi
echo "! $name.bamg created" 1>&2

if $split; then
  filter="| ${pkgbindir}/geo_split | geo -upgrade -geo -"
else
  filter=""
fi
command="bamg2geo $name.bamg $name.dmn $sys_coord_opt $filter > $name.geo"
$verbose && echo "! $command" 1>&2
eval    $command
status=$?
if test $status -ne 0; then
  echo "$0: command failed" 1>&2
  exit $status
fi
echo "! $name.geo created" 1>&2
if $clean; then
  command="rm -f $to_clean"
  $verbose && echo "! $command" 1>&2
  eval    $command
fi
