#!/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_obstacle mesh of the obstacle-in-channel
#@addindex command `mkgeo_contraction`
#@addindex command `geo`
#@addindex command `bamg`
#@addindex mesh
#@addindex file format `.geo` mesh
#
#Synopsis
#--------
#
#  mkgeo_obstacle [options] [nx [ny]]
#
#Example
#-------
#The following command build a triangle-based 2d unstructured mesh
#of the [-L,L]x[0,c] with a unit circle hole at origin.
#
#        mkgeo_obstacle > obstacle.geo
#        geo obstacle.geo
#
#Description
#-----------
#This command is convenient for building a mesh for
#the obstacle-in-channel geometry, 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
#------------
#
#`-c` *float* \n
#`-L` *float*
#>	These options control the geometry parameters.
#>	Default is c=2 and L=2*c.
#
#`-cartesian` \n
#`-zr`
#>	These options control the geometry coordinate system.
#>	Default is Cartesian.
#
#`-quarter` \n
#`-half`
#>	When the `quarter` option is activated, the upstream-downstream
#>	symmetry is assumed and only one quarter of the geometry
#>	is considered.
#>	Otherwise, with the `half` option,
#>	the upstream-downstream geometry is generated
#>	Default is the `quarter` option.
#
#`-N` *int*
#>	The number of elements around the half of the circular obstacle.
#>	When the `quarter` option is activated, only a quarter  
#>	of the circular obstacle is considered.
#>	Default is `N=100`.
#
#The discretization
#------------------
#`-hmin` *float*
#>	Controls the edge length on the obstacle.
#>	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: `obstacle`, `axis`, `wall`,
#`upstream` and `downstream`.
#
#Others options
#--------------
#`-name` *string*
#>	Set the basename for the output files.
#>	By default, the basename is `obstacle`.
#
#`-[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
N=100
c=2
L="unset"
hmin=0.05
split=false
quarter=1
clean=true
verbose=true
name="obstacle"
usage="usage: mkgeo_obstacle 
	[nx=$nx_default [ny=nx]]
	[-c float=$c]
	[-L float=$L]
	[-N float=$N]
	[-quarter|-half]
	[-rz]
        [-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;;
  -L)   L="$2"; shift;;
  -N)   N="$2"; shift;;
  -quarter)   quarter=1;;
  -half)      quarter=0;;
  -zr|-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
if test $L = "unset"; then
  L=`echo | awk -v c=$c '{print c < 2 ? 2 : 2*c}'`
fi
c_error=`echo | awk -v c=$c '{print (c <= 1 ? "true" : "false")}'`
if $c_error; then
  echo "$0: invalid c=$c" 1>&2 ; exit 1
fi
to_clean=""

#
# background mesh for bamg mesh generator:
# with anisotropic (1/hx^2, 1/hy^2) metric
#
# half geometry:
# +-------------------+-------------------+ y=c
# |6+N               5+N               4+N|
# |                                       |
# |                  ---                  | y=1
# |                 /   \                 |
# |1              2|     |2+N          3+N|
# +----------------+  +  +----------------+ y=0
# x=-L          x=-1 x=0 x=1            x=L
#
# quarter geometry:
#                     +-------------------+ y=c
#                    4+N               3+N|
#                     |                   |
#                    1+-                  | y=1
#                       \                 |
#                        |1+N          2+N|
#                     +  +----------------+ y=0
#                    x=0 x=1            x=L
#
awk -v L=$L -v c=$c -v N=$N -v quarter=$quarter '
BEGIN {
  pi = 3.14159265358979323846
  print "MeshVersionFormatted"
  print "  0"
  print "Dimension"
  print "  2"
  print "Vertices"
  print " ", (!quarter ? 6+N : 4+N)
  if (! quarter) {
    printf ("  %.16g %.16g %d\n", -L, 0, 1)
    printf ("  %.16g %.16g %d\n", -1, 0, 2)
    for (i = 1; i <= N; ++i) {
      theta = pi - pi*(1.0*i/N)
      printf ("  %.16g %.16g %d\n", cos(theta), sin(theta), 2+i)
    }
    printf ("  %.16g %.16g %d\n",  L, 0, 3+N)
    printf ("  %.16g %.16g %d\n",  L, c, 4+N)
    printf ("  %.16g %.16g %d\n",  0, c, 5+N)
    printf ("  %.16g %.16g %d\n", -L, c, 6+N)
  } else {
    printf ("  %.16g %.16g %d\n", 0, 1, 1)
    for (i = 1; i <= N; ++i) {
      theta = pi/2 - pi/2*(1.0*i/N)
      printf ("  %.16g %.16g %d\n", cos(theta), sin(theta), 1+i)
    }
    printf ("  %.16g %.16g %d\n",  L, 0, 2+N)
    printf ("  %.16g %.16g %d\n",  L, c, 3+N)
    printf ("  %.16g %.16g %d\n",  0, c, 4+N)
  }
  print "Edges"
  print " ", (!quarter ? 6+N : 4+N)
  if (! quarter) {
    print "  1 2 101"
    for (i = 0; i < N; ++i) {
      print " ", 2+i, 2+i+1, 102
    }
    print " ", 2+N, 3+N, 101
    print " ", 3+N, 4+N, 103
    print " ", 4+N, 5+N, 104
    print " ", 5+N, 6+N, 104
    print " ", 6+N, 1,   105
  } else {
    for (i = 0; i < N; ++i) {
      print " ", 1+i, 1+i+1, 101
    }
    print " ", 1+N, 2+N, 102
    print " ", 2+N, 3+N, 103
    print " ", 3+N, 4+N, 104
    print " ", 4+N, 1,   105
  }
}
' \
> $name.bamgcad
$verbose && echo "! $name.bamgcad created" 1>&2

if test $quarter -eq 0; then
cat > $name.dmn << EOF1
EdgeDomainNames
  5
  axis
  obstacle
  downstream
  wall
  upstream
EOF1
else
cat > $name.dmn << EOF2
EdgeDomainNames
  5
  obstacle
  axis
  downstream
  wall
  vaxis
EOF2
fi
$verbose && echo "! $name.dmn created" 1>&2

#echo "! nx=$nx ny=$ny hmin=$hmin" 1>&2

# -------------------------------------
# metric for adaptive mesh
# -------------------------------------
awk -v L=$L -v c=$c -v N=$N -v hmin=$hmin -v nx=$nx -v ny=$ny -v quarter=$quarter '
BEGIN {
  hm = 0.1
  hx = hm*L/nx
  hy = hm*c/ny
  h0 = hmin # (hmin >= 0.5*(c-1)) ? 0.5*(c-1) : hmin
  hf = (c-1)*h0
  h1 = (c-1 < 1) ? hf : h0
  m0 = h0**(-2)
  m1 = h1**(-2)
  mf = hf**(-2)
  mx = hx**(-2)
  my = hy**(-2)
  print (!quarter ? 6+N : 4+N), 3
  if (!quarter) {
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
    for (i = 0; i <= N; ++i) {
      # TODO: N -> N/2
      theta = (i < N/2) ? 1.0*i/(0.5*N) : 1.0*(N-i)/(0.5*N); 
      hi = h1*theta + h0*(1-theta)
      mi = hi**(-2)
      printf ("  %.16g %.16g %.16g\n", mi, 0, mi)
    }
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
    printf ("  %.16g %.16g %.16g\n", mf, 0, mf)
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
  } else {
    for (i = 0; i <= N; ++i) {
      theta = 1.0*(N-i)/N; 
      hi = h1*theta + h0*(1-theta)
      mi = hi**(-2)
      printf ("  %.16g %.16g %.16g\n", mi, 0, mi)
    }
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
    printf ("  %.16g %.16g %.16g\n", mx, 0, my)
    printf ("  %.16g %.16g %.16g\n", mf, 0, mf)
  }
}
' \
> $name.mtr
$verbose && echo "! $name.mtr created" 1>&2

to_clean="$to_clean $name.mtr"

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
