! Copyright (c) 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019 Harald Klimach <harald.klimach@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! Copyright (c) 2012-2013 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2016 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012, 2015-2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013 Daniel Harlacher <d.harlacher@grs-sim.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2014 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2015-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2015-2016 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017 Daniel Petró <daniel.petro@student.uni-siegen.de>
!
! Parts of this file were written by Harald Klimach, Simon Zimny and Manuel
! Hasert for German Research School for Simulation Sciences GmbH.
!
! Parts of this file were written by Harald Klimach, Kannan Masilamani,
! Daniel Harlacher, Kartik Jain, Verena Krupp, Jiaxing Qi, Peter Vitt,
! Daniel Fleischer, Tobias Girresser and Daniel Petró for University Siegen.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! This file contains the source code for growing and dynamic arrays.
! This is used for arrays of primitives (int, long_int, real, ...) as well as
! for arrays of derived datatypes (tem_variable_type,...).
!
! To use these macros include the following to your source file.
!
! Smart growing array (GA) for ?tstring?
! Growing Arrays:
!
! declaration
!
!
! implementation
!

! -----------------------------------------------------------------
! 2d Array, which can grow in second dimension only (GA2d)
! tname ... indicates type of dynamic array (long, int, real, ...)


!
!------------------------------------------------------------------------------
!
! dynamic Arrays:
!
! declaration
!
!
! implementation
!
!> \brief Module to describe triangles. Contains triangle type
!! with vertices of the triangle
module tem_triangle_module
  use env_module,                only: rk, minLength, zeroLength, eps
  use tem_aux_module,            only: tem_abort
  use tem_logging_module,        only: logunit
  use tem_math_module,           only: cross_product3D
  use tem_cube_module,           only: tem_cube_type
  use tem_transformation_module, only: tem_transformation_type

  ! include aotus modules
  use aotus_module,     only: aot_get_val, aoterr_Fatal, aoterr_WrongType,     &
    &                         flu_State
  use aot_table_module, only: aot_table_open, aot_table_close, aot_table_length
  use aot_out_module,   only: aot_out_type, aot_out_val,                       &
    &                         aot_out_open_table, aot_out_close_table

  implicit none

  private

  public :: grw_triangleArray_type
  public :: tem_triangle_type
  public :: init, append, truncate, destroy, empty, placeAt
  public :: tem_triangleCubeOverlap
  public :: tem_load_triangle
  public :: tem_triangle_out

  !> This type contains three vertices of the triangle
  type tem_triangle_type
    !> 1st index contains x,y,z coordinates
    !! and 2nd index the vertex number.
    real(kind=rk) :: nodes(3,3)
  end type tem_triangle_type

  !> growing array type for type(tem_triangle_type)
  type grw_trianglearray_type
    integer :: nvals = 0
    integer :: containersize = 0
    type(tem_triangle_type), allocatable :: val(:)
  end type

  !> initialize the dynamic array
  interface init
    module procedure init_ga_triangle
  end interface

  !> truncate the array, meaning
  !! cut off the trailing empty entries
  interface truncate
    module procedure truncate_ga_triangle
  end interface

  !> empty the entries  without changing arrays
  interface empty
    module procedure empty_ga_triangle
  end interface

  !> destroy the dynamic array
  interface destroy
    module procedure destroy_ga_triangle
  end interface

  !> insert an element at a given position
  interface placeat
    module procedure placeat_ga_triangle
    module procedure placeat_ga_triangle_vec
  end interface

  !> append a value to the dynamic array
  !! and return its position.
  interface append
    module procedure append_ga_triangle
    module procedure append_ga_triangle_vec
  end interface

  !> increase the size of the container
  !! for the array.
  interface expand
    module procedure expand_ga_triangle
  end interface


  !> interface to write out triangles in lua format to a file
  interface tem_triangle_out
    module procedure tem_triangle_out_scal
    module procedure tem_triangle_out_vec
  end interface tem_triangle_out

  !> interface to load triangles
  interface tem_load_triangle
    module procedure tem_load_triangle
    module procedure tem_load_triangle_single
  end interface tem_load_triangle

contains

  ! ****************************************************************************
  !> Load triangle information from config file.
  subroutine tem_load_triangle(me, transform, conf, thandle )
    !--------------------------------------------------------------------------!
    !inferface variables
    !> array of triangles
    type(tem_triangle_type), allocatable, intent(out) :: me(:)
    !> transformation for spatial object
    type(tem_transformation_type), intent(in) :: transform
    !> lua state
    type(flu_state) :: conf
    integer, intent(in) :: thandle !< handle for canonical objects
    !--------------------------------------------------------------------------!
    ! local varaibles
    integer :: tri_handle, tri_subHandle
    integer :: iObj, nObjects
    !--------------------------------------------------------------------------!

    write(logunit(1),*) 'Loading triangle: '

    call aot_table_open(L = conf, parent = thandle, thandle = tri_handle, &
      &                 key = 'object')
    call aot_table_open(L=conf, parent = tri_handle, thandle = tri_subHandle, &
      & pos = 1 )

    if ( tri_subHandle .eq. 0) then
      !object is a single table
      allocate(me(1))
      call aot_table_close(L=conf, thandle=tri_subHandle)
      call tem_load_triangle_single(me(1), transform, conf, tri_handle)
    else
      !object is a multiple table
      call aot_table_close(L=conf, thandle=tri_subHandle)
      nObjects = aot_table_length(L=conf, thandle=tri_handle)
      allocate(me(nObjects))
      do iObj=1,nObjects
        write(logunit(2),*) '  triangle ', iObj
        call aot_table_open(L=conf, parent=tri_handle, thandle=tri_suBHandle,&
          & pos=iObj)
        call tem_load_triangle_single(me(iObj), transform, conf, tri_Subhandle)
        call aot_table_close(L=conf, thandle=tri_subHandle)
        write(logunit(2),*) ''
      end do
    end if

    call aot_table_close(L=conf, thandle=tri_Handle)

  end subroutine tem_load_triangle
  ! ****************************************************************************


  ! ****************************************************************************
  !> Load single triangle from config file
  subroutine tem_load_triangle_single(me, transform, conf, thandle )
    !--------------------------------------------------------------------------!
    !inferface variables
    !> single triangle
    type(tem_triangle_type), intent(out) :: me
    !> transformation for spatial object
    type(tem_transformation_type), intent(in) :: transform
    !> lua state
    type(flu_state) :: conf
    integer, intent(in) :: thandle !< handle for canonical objects
    !--------------------------------------------------------------------------!
    ! local varaibles
    integer :: vError(3), errFatal(3)
    integer :: iNode
    integer :: sub_handle
    !--------------------------------------------------------------------------!
    errFatal = aoterr_fatal

    call aot_table_open(L=conf, parent=thandle, thandle=sub_handle, &
      &                 key='nodes')
    do iNode=1,3
      call aot_get_val(L=conf, thandle=sub_handle, &
        &              val=me%nodes(:,iNode), &
        &              pos=iNode, ErrCode=vError)

      if (any(btest(vError, errFatal))) then
        write(logunit(0),*) ' Error in configuration: triangle node nr ', iNode
        write(logunit(0),*) ' is not defined'
        call tem_abort()
      end if
    end do
    call aot_table_close(L=conf, thandle=sub_handle)

    write(logunit(2),*) '   node 1: ', me%nodes(:,1)
    write(logunit(2),*) '   node 2: ', me%nodes(:,2)
    write(logunit(2),*) '   node 3: ', me%nodes(:,3)

    !apply transformation
    if(transform%active) then
      if(transform%deform%active) then
        do iNode=1,3
          me%nodes(:,iNode) = matmul( transform%deform%matrix, &
            &                                  me%nodes(:,iNode) )
        enddo
      endif
      if(transform%translate%active) then
        do iNode=1,3
          me%nodes(:,iNode) = transform%translate%vec &
            &                         + me%nodes(:,iNode)
        enddo
      endif
    endif

  end subroutine tem_load_triangle_single
  ! ****************************************************************************

  ! ****************************************************************************
  !> Compute, if the triangle intersects the cube.
  function tem_triangleCubeOverlap(triangle, cube) result(overlaps)
    !--------------------------------------------------------------------------!
    type(tem_triangle_type), intent(in) :: triangle
    type(tem_cube_type), intent(in) :: cube
    logical :: overlaps
    !--------------------------------------------------------------------------!
    real(kind=rk) :: halfwidth
    !--------------------------------------------------------------------------!

    ! halfwidth is increased by eps to avoid precision problem in
    ! triangöe box intersection.
    halfwidth = cube%halfwidth + eps
    overlaps = triBoxOverlap_loc( cube%center,                       &
      &                           [halfwidth, halfwidth, halfwidth], &
      &                           triangle%nodes )

  end function tem_triangleCubeOverlap
  ! ****************************************************************************

  ! ***************************************************************************
  !> This routine checks for triangle box overlap
  !!
  !! this routine is conversion of c-code tribox3.c triBoxOverlap function.
  !! use separating axis theorem to test overlap between triangle and box
  !!  need to test for overlap in these directions:
  !!  1) the {x,y,z}-directions (actually, since we use the AABB of the triangle
  !!     we do not even need to test these)
  !!  2) normal of the triangle
  !!  3) separating axis test. crossproduct(edge from tri, {x,y,z}-directin)
  !!     this gives 3x3=9 more 7tests
  !! This code is available at:
  !! http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/code/tribox3.txt
  function triBoxOverlap_loc( boxCenter, boxHalfwidth, triNodes ) &
    &                         result(overlaps)
    !--------------------------------------------------------------------------!
    !> box center
    real(kind=rk), intent(in) :: boxCenter(3)
    !> halfwidth of the box
    real(kind=rk), intent(in) :: boxHalfwidth(3)
    !> nodes of the triangle
    !! 1st index denote x,y,z coordinates and
    !! 2nd index denote nodes
    real(kind=rk), intent(in) :: triNodes(3,3)
    logical :: overlaps
    !--------------------------------------------------------------------------!
    real(kind=rk) :: nodes(3,3)
    real(kind=rk) :: dirVec1(3), dirVec2(3), dirVec3(3)
    real(kind=rk) :: mmin, mmax
    real(kind=rk) :: normal(3), edge1(3), edge2(3), edge3(3)
    !--------------------------------------------------------------------------!
    overlaps = .false.

    dirVec1 = (/ 1.0_rk, 0.0_rk, 0.0_rk /)
    dirVec2 = (/ 0.0_rk, 1.0_rk, 0.0_rk /)
    dirVec3 = (/ 0.0_rk, 0.0_rk, 1.0_rk /)

    ! move everything so that the box center is in (0,0,0)
    nodes(:,1) = triNodes(:,1) - boxCenter
    nodes(:,2) = triNodes(:,2) - boxCenter
    nodes(:,3) = triNodes(:,3) - boxCenter

    ! compute triangle edges
    edge1 = nodes(:,2) - nodes(:,1)
    edge2 = nodes(:,3) - nodes(:,2)
    edge3 = nodes(:,1) - nodes(:,3)


    ! separating axis test
    if(axistest(dirVec1, edge1, nodes, boxhalfwidth)) return
    if(axistest(dirVec1, edge2, nodes, boxhalfwidth)) return
    if(axistest(dirVec1, edge3, nodes, boxhalfwidth)) return

    if(axistest(dirVec2, edge1, nodes, boxhalfwidth)) return
    if(axistest(dirVec2, edge2, nodes, boxhalfwidth)) return
    if(axistest(dirVec2, edge3, nodes, boxhalfwidth)) return

    if(axistest(dirVec3, edge1, nodes, boxhalfwidth)) return
    if(axistest(dirVec3, edge2, nodes, boxhalfwidth)) return
    if(axistest(dirVec3, edge3, nodes, boxhalfwidth)) return


    ! Bullet 1:
    !  first test overlap in the {x,y,z}-directions
    !  find min, max of the triangle each direction, and test for overlap in
    !  that direction -- this is equivalent to testing a minimal AABB around
    !  the triangle against the AABB
    ! test in X-direction
    mmin = min( nodes(1,1), nodes(1,2), nodes(1,3) )
    mmax = max( nodes(1,1), nodes(1,2), nodes(1,3) )
    if ( mmin > boxhalfwidth(1) .or. mmax < -boxhalfwidth(1) ) return

    ! test in Y-direction
    mmin = min( nodes(2,1), nodes(2,2), nodes(2,3) )
    mmax = max( nodes(2,1), nodes(2,2), nodes(2,3) )
    if ( mmin > boxhalfwidth(2) .or. mmax < -boxhalfwidth(2) ) return

    ! test in Z-direction
    mmin = min( nodes(3,1), nodes(3,2), nodes(3,3) )
    mmax = max( nodes(3,1), nodes(3,2), nodes(3,3) )
    if ( mmin > boxhalfwidth(3) .or. mmax < -boxhalfwidth(3) ) return


    ! Bullet 2:
    !  test if the box intersects the plane of the triangle
    !  compute plane equation of triangle: normal*x+d=0
    normal = cross_product3D(edge1,edge2)

    if( .not. planeBoxOverlap(normal, nodes(:,1), boxhalfwidth)) return

    overlaps = .true.

  end function triBoxOverlap_loc
  ! ***************************************************************************

  !> This function check whether there is separating axis between triangle
  !! and box.
  !! This function can be optimized by explicitly providing cross_product result
  !! see example:
  !! http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/code/tribox3.txt
  function Axistest( dirVec, edge, nodes, boxhalfwidth ) result (success)
    !--------------------------------------------------------------------------!
    real(kind=rk), intent(in) :: dirVec(3)
    real(kind=rk), intent(in) :: edge(3)
    real(kind=rk), intent(in) :: nodes(3,3)
    real(kind=rk), intent(in) :: boxhalfwidth(3)
    logical :: success
    !--------------------------------------------------------------------------!
    real(kind=rk) :: vecA(3), p1, p2, p3, pmin, pmax, rad
    !--------------------------------------------------------------------------!

    success = .true.

    ! a_(i,j) = e_i x f_j
    vecA = cross_product3D( dirVec, edge )

    p1 = dot_product( vecA, nodes(:,1) )
    p2 = dot_product( vecA, nodes(:,2) )
    p3 = dot_product( vecA, nodes(:,3) )

    pmin = min(p1, p2, p3)
    pmax = max(p1, p2, p3)

    rad = abs(vecA(1)) * boxhalfwidth(1) &
      & + abs(vecA(2)) * boxhalfwidth(2) &
      & + abs(vecA(3)) * boxhalfwidth(3)

    if (pmin > rad .or. pmax < -rad) return

    success = .false.

  end function Axistest
  ! ***************************************************************************

  ! ***************************************************************************
  !> This routine checks for plane box overlap
  !> this routine is conversion of c-code tribox3.c planeBoxOverlap function
  function planeBoxOverlap( normal, origin, boxhalfwidth ) result(overlaps)
    !--------------------------------------------------------------------------!
    !> normal direction of the plane
    real(kind=rk), intent(in) :: normal(3)
    !> origin of the plane
    real(kind=rk), intent(in) :: origin(3)
    !> halfwidth of the box
    real(kind=rk), intent(in) :: boxhalfwidth(3)
    logical :: overlaps
    !--------------------------------------------------------------------------!
    integer :: iDir
    real(kind=rk) :: vmin(3), vmax(3), tmp
    !--------------------------------------------------------------------------!
    overlaps = .false.
    ! find min and max of x,y,z of distance between boxhalfwidth and origin
    ! depends on the direction of the plane normal
    do iDir=1,3
      tmp = origin(iDir)
      if(normal(iDir) > 0.0_rk) then
        vmin(iDir) = - boxhalfwidth(iDir) - tmp
        vmax(iDir) =   boxhalfwidth(iDir) - tmp
      else
        vmin(iDir) =   boxhalfwidth(iDir) - tmp
        vmax(iDir) = - boxhalfwidth(iDir) - tmp
      endif
    end do

    if(dot_product(normal, vmin) > 0.0_rk) return

    if(dot_product(normal, vmax) >= 0.0_rk) overlaps = .true.

  end function planeBoxOverlap
  ! ***************************************************************************

  ! ************************************************************************** !
  !> Write out an array of triangles in lua format
  !! Only if nTriangles <= 10
  !!
  subroutine tem_triangle_out_vec( me, conf )
    ! --------------------------------------------------------------------------
    !> triangle types to write out
    type( tem_triangle_type ), intent(in) :: me(:)
    !> Aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf
    ! --------------------------------------------------------------------------
    ! counter
    integer :: i
    ! --------------------------------------------------------------------------

    ! create a table with name triangle
    call aot_out_open_table( put_conf = conf, tname = 'triangle' )

    if (size(me)<=10) then
      do i = 1, size(me)
        call tem_triangle_out_scal( me(i), conf )
      end do
    else
      call aot_out_val( put_conf = conf, vname = 'nTriangles', val = size(me) )
      call aot_out_val( put_conf = conf,                                  &
        &               val = 'Write triangle is limited to 10 triangles' )
    end if

    call aot_out_close_table( put_conf = conf )

  end subroutine tem_triangle_out_vec
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Write out a triangle shape in lua format
  !!
  subroutine tem_triangle_out_scal( me, conf )
    ! --------------------------------------------------------------------------
    !> triangle types to write out
    type( tem_triangle_type ), intent(in) :: me
    !> Aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf
    ! --------------------------------------------------------------------------
    integer :: iNode
    ! ---------------------------------------------------------------------------

    ! create a table with name triangle if not exist
    if( conf%level == 0 ) then
      call aot_out_open_table( put_conf = conf, tname = 'triangle' )
    else
      call aot_out_open_table( put_conf = conf )
    end if

    !write nodes
    call aot_out_open_table( put_conf = conf, tname = 'nodes' )
    do iNode=1,3
      call aot_out_val( put_conf = conf, val = me%nodes(:,iNode) )
    enddo
    call aot_out_close_table( put_conf = conf )

    call aot_out_close_table( put_conf = conf )

  end subroutine tem_triangle_out_scal
  ! ************************************************************************** !

  subroutine init_ga_triangle(me, length)
    type(grw_trianglearray_type), intent(out) :: me !< dynamic array to init
    integer, intent(in), optional :: length !< initial length of the container

    if (present(length)) then
      me%containersize = length
    else
      me%containersize = zerolength
    end if
    ! deallocate ...
    if( allocated( me%val ))     &
      deallocate(me%val)
    ! ... and reallocate
    allocate(me%val(me%containersize))
    me%nvals = 0

  end subroutine init_ga_triangle

  subroutine destroy_ga_triangle(me)
    type(grw_trianglearray_type), intent(inout) :: me !< dynamic array to destroy

    me%containersize = 0
    me%nvals = 0
    if( allocated( me%val ) ) deallocate(me%val)

  end subroutine destroy_ga_triangle


  subroutine truncate_ga_triangle(me)
    !------------------------------------------------------------------------
    type(grw_trianglearray_type) :: me !< array to truncate
    !------------------------------------------------------------------------
    type(tem_triangle_type), allocatable :: tarray(:)
    !------------------------------------------------------------------------
    integer :: ii
    !------------------------------------------------------------------------

    ! nothing to do if container size is not larger than the number of values
    ! in the array.
    if (me%containersize > me%nvals) then
      allocate(tarray(me%nvals))
      do ii = 1, me%nvals
        tarray(ii) = me%val(ii)
      end do
      call move_alloc(tarray, me%val)
      me%containersize = me%nvals
    end if

  end subroutine truncate_ga_triangle


  subroutine empty_ga_triangle(me)
    !------------------------------------------------------------------------
    type(grw_trianglearray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------

    me%nvals = 0

  end subroutine empty_ga_triangle

  !> adds the value to a given position inside the growing array.
  !!
  !! if the requested position is outside the current array bounds, the array
  !! will be resized accordingly. if it is inside the current array bounds, the
  !! element at the requested position will be replaced.
  subroutine placeat_ga_triangle(me, val, pos, length)
    type(grw_trianglearray_type) :: me !< array to place the value into
    type(tem_triangle_type), intent(in) :: val !< value to place at the given position
    integer, intent(in) :: pos !< predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length


    ! value to append is larger than all existing ones,
    ! just put it to the end of the list, this captures
    ! also the case of empty lists.
    ! in this case foundpos = me%nvals + 1 holds.
    if (pos > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = pos, length = length)
    end if

    me%nvals = max( pos, me%nvals )
    me%val(pos) = val

  end subroutine placeat_ga_triangle


  !> adds the values starting from a given position inside the growing array.
  !!
  !! if the requested position is outside the current array bounds, the array
  !! will be resized accordingly. if it is inside the current array bounds, the
  !! elements starting from the requested position will be replaced up to
  !! the element at position `pos + size(val) - 1`.
  subroutine placeat_ga_triangle_vec(me, val, pos, length)
    type(grw_trianglearray_type) :: me !< array to append the value to
    type(tem_triangle_type), intent(in) :: val(:) !< values to append
    integer, intent(in) :: pos !< predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: ub, ii

    if (me%nvals == huge(me%nvals)) then
      write(*,*) "reached end of integer range for growing array!"
      write(*,*) "aborting!!"
      stop
    end if

    ub = pos + size(val) - 1

    if (ub > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = ub, length = length)
    end if

    me%nvals = max( ub, me%nvals )
    do ii = pos, ub
      me%val(ii) = val(1+ii-pos)
    end do

  end subroutine placeat_ga_triangle_vec


  subroutine append_ga_triangle(me, val, length)
    type(grw_trianglearray_type) :: me !< array to append the value to
    type(tem_triangle_type), intent(in) :: val !< value to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    ! value to append is larger than all existing ones,
    ! just put it to the end of the list, this captures
    ! also the case of empty lists.
    ! in this case foundpos = me%nvals + 1 holds.
    if (me%nvals+1 > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, length = length)
    end if

    me%nvals = me%nvals+1
    me%val(me%nvals) = val

  end subroutine append_ga_triangle

  subroutine append_ga_triangle_vec(me, val, length)
    type(grw_trianglearray_type) :: me !< array to append the value to
    type(tem_triangle_type), intent(in) :: val(:) !< values to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: lb, ub, ii

    if (me%nvals == huge(me%nvals)) then
      write(*,*) "reached end of integer range for growing array!"
      write(*,*) "aborting!!"
      stop
    end if

    lb = me%nvals + 1
    ub = lb + size(val) - 1

    if (ub > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = ub, length = length)
    end if

    me%nvals = max( ub, me%nvals )
    do ii = lb, ub
      me%val(ii) = val(1+ii-lb)
    end do

  end subroutine append_ga_triangle_vec


  subroutine expand_ga_triangle(me, pos, length)
    type(grw_trianglearray_type) :: me !< array to resize
    integer, intent(in), optional :: pos !< optional predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length

    type(tem_triangle_type), allocatable :: swpval(:)
    integer :: explen, ii

    explen = 0
    ! increase the container by the requested length of double it
    if( present(length) ) then
      explen = max( length, minlength )
    else
      ! set the global minimum length, if doubling would be smaller than that
      explen = max(me%containersize, minlength)
    end if

    ! if a position is given, increase the container to at least the size to
    ! fit the position.
    if( present(pos) ) explen = max(explen, pos-me%containersize)

    ! if the current size plus explen exceeds the max container size,
    ! reduce the size to the max container size.
    if( (huge(me%containersize) - explen) <= me%containersize) then
      ! set max container size
      me%containersize = huge(me%containersize)
    else
      ! set the new container size
      me%containersize = me%containersize + explen
    end if

    if ( me%nvals > 0 ) then
      allocate(swpval(me%containersize))
      do ii = 1, me%nvals
        swpval(ii) = me%val(ii)
      end do
      call move_alloc( swpval, me%val )
    else ! me%nvals == 0
      if ( allocated(me%val) ) deallocate( me%val )
      allocate( me%val(me%containersize) )
    end if

  end subroutine expand_ga_triangle


end module tem_triangle_module

!> \page triangle Triangle
!! Triangles are defined by three nodes. It is possible to define triangles
!! as a geometry kind. Also, the planes and hollow boxes are converted to
!! triangles internally.
!!
!! Valid definition:
!! \li Single triangle
!! \verbatim
!! geometry={
!!   kind='triangle',
!!     object={
!!       nodes = {             -- triangle 1
!!         {0.0, 0.0, 0.5},    -- node 1
!!         {-0.5, -0.5, -0.5}, -- node 2
!!         {0.5, -0.5, -0.5}      -- node 3
!!       }
!!     }
!! }
!! \endverbatim
!!
!! \li Multiple triangles
!! \verbatim
!! geometry={
!!   kind='triangle',
!!     object={
!!       {
!!         nodes = {             -- triangle 1
!!           {0.0, 0.0, 0.5},    -- node 1
!!           {-0.5, -0.5, -0.5}, -- node 2
!!           {0.5, -0.5, -0.5}      -- node 3
!!         }
!!       },
!!       {
!!         nodes = {             -- triangle 2
!!           {0.0, 0.0, 0.5},    -- node 1
!!           {-0.5, -0.5, -0.5}, -- node 2
!!           {-0.5, 0.5, -0.5}      -- node 3
!!         }
!!       }
!!     }
!! }
!! \endverbatim
!! \n\n
!! A pyramid has created with four triangles and a plane. The corresponding
!! seeder code is given below:
!! \include testsuite/triangle/seeder.lua
!! \n\n
!! Mesh with pyramid generated by the seeder file:
!! \image html triangle.png
!! \image html triangle_withedges.png
