text:
!> IMPACT
!! \author Rolf Henniger, Institute of Fluid Dynamics, ETH Zurich (henniger@ifd.mavt.ethz.ch)
!! \date Mai 2005 - Dec 2011

!> \brief module providing functions to initiliaze and apply RestrictionOp
module cmod_InterpolationOp

  use iso_c_binding
  use mpi

  implicit none

contains

  subroutine MG_getCIS( &
      Nc,               &
      Nf,               &
      bL,               &
      bU,               &
      xf,               &
      dd,               &
      cI ) bind(c,name='MG_getCIS')

    implicit none

    integer(c_int), intent(in)  :: Nc
    integer(c_int), intent(in)  :: Nf
    integer(c_int), intent(in)  :: bL, bU

    real(c_double), intent(in)  :: xf(bL:(Nf+bU))

    integer(c_int), intent(in)  :: dd

    real(c_double), intent(out) :: cI(1:2,1:Nc)

    integer(c_int)              :: i,ic
    real(c_double)              :: Dx12

    cI = 0.

    !====================================================================================
    !=== Interpolation, linienweise, 1d =================================================
    !====================================================================================
    ! coarse
    !      ic                             ic+1
    !  ----o------------------------------o------------------------------o---------------
    !      |------Dx10----|
    !      |------------Dx12--------------|
    !  ----o--------------o---------------o-----------o-----
    !     xf(i-1)        xf(i)           xf(i+1)
    ! fine

    !------------------------------------------------------------------------------------
    do ic = 1, Nc

      i = dd * ic
      Dx12 = xf(i+1)-xf(i-1)

      cI(1,ic) = ( xf(i+1)-xf(i  ) )/Dx12
      cI(2,ic) = ( xf(i  )-xf(i-1) )/Dx12

    end do
    !------------------------------------------------------------------------------------

  end subroutine MG_getCIS

  !> \brief gets coefficents for interpolation of velocity points
  !!
  !! \param[in] Nc grid size on coarse grid
  !! \param[in] bLc lower ghost width on coarse grid
  !! \param[in] bUc upper ghost width on coarse grid
  !! \param[in] SSc start index on coarse grid
  !! \param[in] NNc end index on coarse grid
  !! \param[in] BC_L lower boundary conditions
  !! \param[in] BC_U upper boundary conditions
  !! \param[in] Nf grid size on fine grid
  !! \param[in] bLf lower ghost width on fine grid
  !! \param[in] bUf upper ghost width on fine grid
  !! \param[in] offset offset
  !! \param[in] xc coarse coordinates of velocity
  !! \param[in] xf fine coordinates of velocity
  !! \param[in] dd coarsening factor
  !! \param[out] cIV interpolation coefficients
  subroutine MG_getCIV( &
      Nc,               &
      bLc,              &
      bUc,              &
      !SSc,              &
      !NNc,              &
      BC_L,             &
      BC_U,             &
      Nf,               &
      bLf,              &
      bUf,              &
      offset,           &
      xc,               &
      xf,               &
      dd,               &
      cIV ) bind(c,name='MG_getCIV')

    implicit none

    integer(c_int), intent(in)    :: Nc

    integer(c_int), intent(in)    :: bLc, bUc

    !integer(c_int), intent(in)    :: SSc
    !integer(c_int), intent(in)    :: NNc

    integer(c_int), intent(in)    :: BC_L,BC_U

    integer(c_int), intent(in)    :: Nf

    integer(c_int), intent(in)    :: bLf, bUf

    integer(c_int), intent(in)    :: offset

    real(c_double), intent(in)    :: xc( bLc:(Nc+bUc) )

    real(c_double), intent(in)    :: xf( bLf:(Nf+bUf) )

    integer(c_int), intent(in)    :: dd

    real(c_double), intent(inout) :: cIV( 1:2, 0:Nf )

    integer(c_int)                ::  i, ic
    !real(c_double)                ::  Dx1a
    real(c_double)                ::  Dx12

    !====================================================================================
    !=== Interpolation, linienweise, 1d =================================================
    !====================================================================================
    ! fine
    !    xf(i-1)        xf(i)         xf(i+1)       xf(i+1)      xf(i+2)
    !  ---->-------------->--------------->----------->----------->------------
    !             |-----------Dx12--------------|
    !  ----------->----------------------------->----------------------------->----------
    !          xc(ic)                      xc(ic+1)                        xc(ic+1)
    ! coarse

    cIV = 0.

    do i = 0, Nf
      ic = ( i )/dd + offset

      !Dx1a = xc(ic)-xf(i   )
      !Dx12 = xc(ic)-xc(ic+1)

      !cIV(1,i) = 1-Dx1a/Dx12
      !cIV(2,i) =   Dx1a/Dx12

      Dx12 = xc(ic+1) - xc(ic)

      cIV(1,i) = ( xc(ic+1) - xf(i ) )/Dx12
      cIV(2,i) = ( xf(i   ) - xc(ic) )/Dx12

    end do

    ! a little bit shaky, please verify this when IO is ready
    !if (BC_L > 0) then
    !cIV( 1,0) = 0.
    !cIV( 2,0) = 1.
    !end if

    !maybe false
    if( BC_L == -2 )then
      cIV(1:2,0) = 0.
    end if

    !if (BC_U > 0) then
    !cIV(1,Nf) = 1.
    !cIV(2,Nf) = 0.
    !end if
    if( BC_U == -2 )then
      cIV(1:2,Nf) = 0.
    end if

  end subroutine MG_getCIV

  !> \brief corrects values at corners through extrapolation
  !!
  !! \param[in] n
  !! \param[in] bL
  !! \param[in] bU
  !! \param[in] BCL
  !! \param[in] BCU
  !! \param[in] x1
  !! \param[in] x2
  !! \param[in] x3
  !! \param[inout] phi \f$\phi\f$
  !!
  !! \f[ \phi[2] = \frac{ \phi[4]-\phi[3] }{x[4]-x[3]} ( x[2]-x[3] )\f]
  !! \f[ \phi[n-1] = \frac{ \phi[n-2]-\phi[n-3] }{x[n-2]-x[n-3]} ( x[n-1]-x[n-2] )\f]
  subroutine MG_InterpolateCornersPost( &
      n,                                &
      bL,                               &
      bU,                               &
      BCL,                              &
      BCU,                              &
      x1,                               &
      x2,                               &
      x3,                               &
      phi ) bind(c,name='MG_InterpolateCornersPost')

    implicit none

    integer(c_int), intent(in)    :: n(3)

    integer(c_int), intent(in)    :: bL(3)
    integer(c_int), intent(in)    :: bU(3)

    integer(c_int), intent(in)    :: BCL(3)
    integer(c_int), intent(in)    :: BCU(3)

    real(c_double), intent(in)    :: x1( bL(1):(n(1)+bU(1)) )
    real(c_double), intent(in)    :: x2( bL(2):(n(2)+bU(2)) )
    real(c_double), intent(in)    :: x3( bL(3):(n(3)+bU(3)) )

    real(c_double), intent(inout) :: phi( bL(1):(n(1)+bU(1)), bL(2):(n(2)+bU(2)), bL(3):(n(3)+bU(3)) )

    integer(c_int) :: s(3), e(3), i

    do i = 1, 3

      if( BCL(i)>0 ) then
        s(i) = 2
      else
        s(i) = 1
      end if

      if( BCU(i)>0 ) then
        e(i) = n(i)-1
      else
        e(i) = n(i)
      end if

    end do

    if( BCL(1)>0 .and. BCL(2)>0 ) then
      phi(1, 1, s(3):e(3)) = 0.
      phi(2, 1, s(3):e(3)) = ( phi( 4, 1, s(3):e(3)) - phi( 3, 1, s(3):e(3)) )/( x1(4) - x1(3) )*( x1(2) - x1(3) ) + phi( 3, 1, s(3):e(3))
      phi(1, 2, s(3):e(3)) = ( phi( 1, 4, s(3):e(3)) - phi( 1, 3, s(3):e(3)) )/( x2(4) - x2(3) )*( x2(2) - x2(3) ) + phi( 1, 3, s(3):e(3))
      phi(2, 2, s(3):e(3)) = ( &
        phi( 1, 2, s(3):e(3))*( x1(3)-x1(2) )/( x1(3)-x1(1) ) + &
        phi( 3, 2, s(3):e(3))*( x1(2)-x1(1) )/( x1(3)-x1(1) ) + &
        phi( 2, 1, s(3):e(3))*( x2(3)-x2(2) )/( x2(3)-x2(1) ) + &
        phi( 2, 3, s(3):e(3))*( x2(2)-x2(1) )/( x2(3)-x2(1) )   &
        )/2.
    endif
    if( BCL(1)>0 .and. BCU(2)>0 ) then
      phi(1, n(2)  , s(3):e(3)) =  0.
      phi(2, n(2)  , s(3):e(3)) = ( phi(4, n(2)  , s(3):e(3)) - phi(3, n(2)  , s(3):e(3)) )/( x1(4)      - x1(3)      )*( x1(2)      - x1(3)      ) + phi(3, n(2)  , s(3):e(3))
      phi(1, n(2)-1, s(3):e(3)) = ( phi(1, n(2)-2, s(3):e(3)) - phi(1, n(2)-3, s(3):e(3)) )/( x2(n(2)-2) - x2(n(2)-3) )*( x2(n(2)-1) - x2(n(2)-2) ) + phi(1, n(2)-2, s(3):e(3))
      phi(2, n(2)-1, s(3):e(3)) = ( &
        phi(1, n(2)-1, s(3):e(3))*( x1(3)     -x1(2)      )/( x1(3)-x1(1)         ) + &
        phi(3, n(2)-1, s(3):e(3))*( x1(2)     -x1(1)      )/( x1(3)-x1(1)         ) + &
        phi(2, n(2)  , s(3):e(3))*( x2(n(2)-1)-x2(n(2)-2) )/( x2(n(2))-x2(n(2)-2) ) + &
        phi(2, n(2)-2, s(3):e(3))*( x2(n(2)  )-x2(n(2)-1) )/( x2(n(2))-x2(n(2)-2) )   &
        )/2.
    endif
    if( BCU(1)>0 .and. BCL(2)>0 ) then
      phi(n(1)  , 1, s(3):e(3)) = 0.
      phi(n(1)-1, 1, s(3):e(3)) = ( phi(n(1)-2, 1, s(3):e(3)) - phi(n(1)-3, 1, s(3):e(3)) )/( x1(n(1)-2) - x1(n(1)-3) )*( x1(n(1)-1) - x1(n(1)-2) ) + phi(n(1)-2, 1, s(3):e(3))
      phi(n(1)  , 2, s(3):e(3)) = ( phi(n(1)  , 4, s(3):e(3)) - phi(n(1)  , 3, s(3):e(3)) )/( x2(4) - x2(3) )*( x2(2) - x2(3) )                     + phi(n(1)  , 3, s(3):e(3))
      phi(n(1)-1, 2, s(3):e(3)) = ( &
        phi(n(1)  , 2, s(3):e(3))*( x1(n(1)-1)-x1(n(1)-2) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-2, 2, s(3):e(3))*( x1(n(1)  )-x1(n(1)-1) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-1, 1, s(3):e(3))*( x2(3)-x2(2) )/( x2(3)-x2(1) ) + &
        phi(n(1)-1, 3, s(3):e(3))*( x2(2)-x2(1) )/( x2(3)-x2(1) )   &
        )/2.
    endif
    if( BCU(1)>0 .and. BCU(2)>0 ) then
      phi(n(1)  , n(2)  , s(3):e(3)) = 0.
      phi(n(1)-1, n(2)  , s(3):e(3)) = ( phi(n(1)-2, n(2)  , s(3):e(3)) - phi(n(1)-3, n(2)  , s(3):e(3)) )/( x1(n(1)-2) - x1(n(1)-3) )*( x1(n(1)-1) - x1(n(1)-2) ) + phi(n(1)-2, n(2)  , s(3):e(3))
      phi(n(1)  , n(2)-1, s(3):e(3)) = ( phi(n(1)  , n(2)-2, s(3):e(3)) - phi(n(1)  , n(2)-3, s(3):e(3)) )/( x2(n(2)-2) - x2(n(2)-3) )*( x2(n(2)-1) - x2(n(2)-2) ) + phi(n(1)  , n(2)-2, s(3):e(3))
      phi(n(1)-1, n(2)-1, s(3):e(3)) = ( &
        phi(n(1)-1, n(2)  , s(3):e(3))*( x2(n(2)-1)-x2(n(2)-2) )/( x2(n(2))-x2(n(2)-2) ) + &
        phi(n(1)-1, n(2)-2, s(3):e(3))*( x2(n(2)  )-x2(n(2)-1) )/( x2(n(2))-x2(n(2)-2) ) + &
        phi(n(1)  , n(2)-1, s(3):e(3))*( x1(n(1)-1)-x1(n(1)-2) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-2, n(2)-1, s(3):e(3))*( x1(n(1)  )-x1(n(1)-1) )/( x1(n(1))-x1(n(1)-2) )   &
        )/2.
    endif

    if( BCL(1)>0 .and. BCL(3)>0 ) then
      phi(1, s(2):e(2), 1) = 0.
      phi(2, s(2):e(2), 1) = ( phi(4, s(2):e(2), 1)-phi(3, s(2):e(2), 1) )/( x1(4)-x1(3) )*( x1(2)-x1(3) ) + phi(3, s(2):e(2), 1)
      phi(1, s(2):e(2), 2) = ( phi(1, s(2):e(2), 4)-phi(1, s(2):e(2), 3) )/( x3(4)-x3(3) )*( x3(2)-x3(3) ) + phi(1, s(2):e(2), 3)
      phi(2, s(2):e(2), 2) = ( &
        phi(2, s(2):e(2), 1)*( x3(3)-x3(2) )/( x3(3)-x3(1) ) + &
        phi(2, s(2):e(2), 3)*( x3(2)-x3(1) )/( x3(3)-x3(1) ) + &
        phi(1, s(2):e(2), 2)*( x1(3)-x1(2) )/( x1(3)-x1(1) ) + &
        phi(3, s(2):e(2), 2)*( x1(2)-x1(1) )/( x1(3)-x1(1) )   &
        )/2.
    endif
    if( BCL(1)>0 .and. BCU(3)>0 ) then
      phi(1, s(2):e(2), n(3)  ) = 0
      phi(2, s(2):e(2), n(3)  ) = ( phi(4, s(2):e(2), n(3)  )-phi(3, s(2):e(2), n(3)  ) )/( x1(4)-x1(3) )*( x1(2)-x1(3) )                     + phi(3, s(2):e(2), n(3)  )
      phi(1, s(2):e(2), n(3)-1) = ( phi(1, s(2):e(2), n(3)-2)-phi(1, s(2):e(2), n(3)-3) )/( x3(n(3)-2)-x3(n(3)-3) )*( x3(n(3)-1)-x3(n(3)-2) ) + phi(1, s(2):e(2), n(3)-2)
      phi(2, s(2):e(2), n(3)-1) = ( &
        phi(2, s(2):e(2), n(3)  )*( x3(n(3)-1)-x3(n(3)-2) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(2, s(2):e(2), n(3)-2)*( x3(n(3)  )-x3(n(3)-1) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(1, s(2):e(2), n(3)-1)*( x1(3)-x1(2) )/( x1(3)-x1(1) ) + &
        phi(3, s(2):e(2), n(3)-1)*( x1(2)-x1(1) )/( x1(3)-x1(1) ) &
        )/2.
    endif
    if( BCU(1)>0 .and. BCL(3)>0 ) then
      phi(n(1)  , s(2):e(2), 1) = 0.
      phi(n(1)-1, s(2):e(2), 1) = ( phi(n(1)-2, s(2):e(2), 1) - phi(n(1)-3, s(2):e(2), 1) )/( x1(n(1)-2)-x1(n(1)-3) )*( x1(n(1)-1)-x1(n(1)-2) ) + phi(n(1)-2, s(2):e(2), 1)
      phi(n(1)  , s(2):e(2), 2) = ( phi(n(1)  , s(2):e(2), 4) - phi(n(1)  , s(2):e(2), 3) )/( x3(4)     -x3(3)      )*( x3(2)     -x3(3)      ) + phi(n(1)  , s(2):e(2), 3)
      phi(n(1)-1, s(2):e(2), 2) = ( &
        phi(n(1)  , s(2):e(2), 2)*( x1(n(1)-1)-x1(n(1)-2) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-2, s(2):e(2), 2)*( x1(n(1)  )-x1(n(1)-1) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-1, s(2):e(2), 1)*( x3(3)     -x3(2)      )/( x3(3)   -x3(1)      ) + &
        phi(n(1)-1, s(2):e(2), 3)*( x3(2)     -x3(1)      )/( x3(3)   -x3(1)      )   &
        )/2.
    endif
    if( BCU(1)>0 .and. BCU(3)>0 ) then
      phi(n(1)  , s(2):e(2)  , n(3)  ) = 0.
      phi(n(1)-1, s(2):e(2), n(3)  ) = ( phi(n(1)-2, s(2):e(2), n(3)  )-phi(n(1)-3, s(2):e(2), n(3)  ) )/( x1(n(1)-2)-x1(n(1)-3) )*( x1(n(1)-1)-x1(n(1)-2) ) + phi(n(1)-2, s(2):e(2), n(3)  )
      phi(n(1)  , s(2):e(2), n(3)-1) = ( phi(n(1)  , s(2):e(2), n(3)-2)-phi(n(1)  , s(2):e(2), n(3)-3) )/( x3(n(3)-2)-x3(n(3)-3) )*( x3(n(3)-1)-x3(n(3)-2) ) + phi(n(1)  , s(2):e(2), n(3)-2)
      phi(n(1)-1, s(2):e(2), n(3)-1) = ( &
        phi(n(1)-1, s(2):e(2), n(3)  )*( x3(n(3)-1)-x3(n(3)-2) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(n(1)-1, s(2):e(2), n(3)-2)*( x3(n(3)  )-x3(n(3)-1) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(n(1)  , s(2):e(2), n(3)-1)*( x1(n(1)-1)-x1(n(1)-2) )/( x1(n(1))-x1(n(1)-2) ) + &
        phi(n(1)-2, s(2):e(2), n(3)-1)*( x1(n(1)  )-x1(n(1)-1) )/( x1(n(1))-x1(n(1)-2) )   &
        )/2.
    endif

    if( BCL(2)>0 .and. BCL(3)>0 ) then
      phi(s(1):e(1), 1, 1 ) = 0.
      phi(s(1):e(1), 2, 1 ) = ( phi(s(1):e(1), 4, 1 )-phi(s(1):e(1), 3, 1 ) )/( x2(4)-x2(3) )*( x2(2)-x2(3) ) + phi(s(1):e(1), 3, 1 )
      phi(s(1):e(1), 1, 2 ) = ( phi(s(1):e(1), 1, 4 )-phi(s(1):e(1), 1, 3 ) )/( x3(4)-x3(3) )*( x3(2)-x3(3) ) + phi(s(1):e(1), 1, 3 )
      phi(s(1):e(1), 2, 2 ) = ( &
        phi(s(1):e(1), 2, 1 )*( x3(3)-x3(2) )/( x3(3)-x3(1) ) + &
        phi(s(1):e(1), 2, 3 )*( x3(2)-x3(1) )/( x3(3)-x3(1) ) + &
        phi(s(1):e(1), 1, 2 )*( x2(3)-x2(2) )/( x2(3)-x2(1) ) + &
        phi(s(1):e(1), 3, 2 )*( x2(2)-x2(1) )/( x2(3)-x2(1) )   &
        )/2.
    endif
    if( BCL(2)>0 .and. BCU(3)>0 ) then
      phi(s(1):e(1), 1, n(3)  ) = 0.
      phi(s(1):e(1), 2, n(3)  ) = ( phi(s(1):e(1), 4, n(3)  )-phi(s(1):e(1), 3, n(3)  ) )/( x2(4)     -x2(3)      )*( x2(2)     -x2(3)      ) + phi(s(1):e(1), 3, n(3)  )
      phi(s(1):e(1), 1, n(3)-1) = ( phi(s(1):e(1), 1, n(3)-2)-phi(s(1):e(1), 1, n(3)-3) )/( x3(n(3)-2)-x3(n(3)-3) )*( x3(n(3)-1)-x3(n(3)-2) ) + phi(s(1):e(1), 1, n(3)-2)
      phi(s(1):e(1), 2, n(3)-1) = ( &
        phi(s(1):e(1), 2, n(3)  )*( x3(n(3)-1)-x3(n(3)-2) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(s(1):e(1), 2, n(3)-2)*( x3(n(3)  )-x3(n(3)-1) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(s(1):e(1), 1, n(3)-1)*( x2(3)     -x2(2)      )/( x2(3)-x2(1)         ) + &
        phi(s(1):e(1), 3, n(3)-1)*( x2(2)     -x2(1)      )/( x2(3)-x2(1)         ) &
        )/2.
    endif
    if( BCU(2)>0 .and. BCL(3)>0 ) then
      phi(s(1):e(1), n(2)  , 1 ) = 0.
      phi(s(1):e(1), n(2)-1, 1 ) = ( phi(s(1):e(1), n(2)-2, 1)-phi(s(1):e(1), n(2)-3, 1) )/( x2(n(2)-2)-x2(n(2)-3) )*( x2(n(2)-1)-x2(n(2)-2) ) + phi(s(1):e(1), n(2)-2, 1)
      phi(s(1):e(1), n(2)  , 2 ) = ( phi(s(1):e(1), n(2)  , 4)-phi(s(1):e(1), n(2)  , 3) )/( x3(4)     -x3(3)      )*( x3(2)     -x3(3)      ) + phi(s(1):e(1), n(2),   3)
      phi(s(1):e(1), n(2)-1, 2 ) = ( &
        phi(s(1):e(1), n(2)-1, 1)*( x3(3)     -x3(2)      )/( x3(3)-x3(1)         ) + &
        phi(s(1):e(1), n(2)-1, 3)*( x3(2)     -x3(1)      )/( x3(3)-x3(1)         ) + &
        phi(s(1):e(1), n(2)  , 2)*( x2(n(2)-1)-x2(n(2)-2) )/( x2(n(2))-x2(n(2)-2) ) + &
        phi(s(1):e(1), n(2)-2, 2)*( x2(n(2)  )-x2(n(2)-1) )/( x2(n(2))-x2(n(2)-2) ) &
        )/2.
    endif
    if( BCU(2)>0 .and. BCU(3)>0 ) then
      phi(s(1):e(1), n(2)  , n(3)  ) = 0.
      phi(s(1):e(1), n(2)-1, n(3)  ) = ( phi(s(1):e(1), n(2)-2, n(3)  )-phi(s(1):e(1), n(2)-3, n(3)  ) )/( x2(n(2)-2)-x2(n(2)-3) )*( x2(n(2)-1)-x2(n(2)-2) ) + phi(s(1):e(1), n(2)-2, n(3)  )
      phi(s(1):e(1), n(2)  , n(3)-1) = ( phi(s(1):e(1), n(2)  , n(3)-2)-phi(s(1):e(1), n(2)  , n(3)-3) )/( x3(n(3)-2)-x3(n(3)-3) )*( x3(n(3)-1)-x3(n(3)-2) ) + phi(s(1):e(1), n(2)  , n(3)-2)
      phi(s(1):e(1), n(2)-1, n(3)-1) = (                                                   &
        phi(s(1):e(1), n(2)-1, n(3)  )*( x3(n(3)-1)-x3(n(3)-2) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(s(1):e(1), n(2)-1, n(3)-2)*( x3(n(3)  )-x3(n(3)-1) )/( x3(n(3))-x3(n(3)-2) ) + &
        phi(s(1):e(1), n(2),   n(3)-1)*( x2(n(2)-1)-x2(n(2)-2) )/( x2(n(2))-x2(n(2)-2) ) + &
        phi(s(1):e(1), n(2)-2, n(3)-1)*( x2(n(2)  )-x2(n(2)-1) )/( x2(n(2))-x2(n(2)-2) )   &
        )/2.
    endif

  end subroutine MG_InterpolateCornersPost

  !> \note:       - für allgemeine dd geeignet
  !!              - dd(i) /= 1 <===> N(i) /= 1
  !!              - es wird nur in eine Richung ausgetauscht
  !!              - Null-Setzen am Rand nicht notwendig
  !!              - Es wird sequentiell über alle Raumrichtungen interpoliert,
  !!                um keinen individuellen Interpolationsstencil für jeden
  !!                Punkt im Raum speichern zu müssen.
  !!              - Durch das sequentielle Interpolieren kann der
  !!                Interpolationsstencil klein und damit der Gesamtaufwand
  !!                minimiert werden (Alternative: 8- bzw. 26-Punkt Stencil(!)).
  !!              - Interpolationskoeffizienten werden auf dem jeweils feineren
  !!                Gitter gespeichert, um nicht auf die entsprechenden Indizes
  !!                des gröberen Gitters umrechnen zu müssen.
  !!              - Die Block-überlappenden Stirnflächen werden ebenfalls
  !!                mitverarbeitet, aber eigentlich nicht gebraucht (erleichtert
  !!                die Programmierung), so dass eigentlich eine Initialisierung
  !!                notwendig wäre. Dies wird jedoch zuvor schon in der
  !!                korrespondierenden Restriktions-Routine erledigt, so dass
  !!                dies hier nicht mehr notwendig ist.
  !!
  !! \attention: Verwendung von parametrischen Strides dd in Schleifen
  !! verhindert die Vektorisierung bzw. das Prefetching! Andererseits ist der
  !! Geschwindigkeitsgewinn nur sehr gering (schon getestet).
  !!
  !! - Geschwindigkeit ist trotz Prefetching / Vektorisierung stark durch die
  !!   Speicherzugriffszeit limitiert
  !! - Das wird z.B. deutlich bei Single- vs.  Dualcorebetrieb
  subroutine MG_interpolate(  &
      Nc,                     &
      bLc,                    &
      bUc,                    &
      Nf,                     &
      bLf,                    &
      bUf,                    &
      iimax,                  &
      dd,                     &
      cI1,                    &
      cI2,                    &
      cI3,                    &
      phic,                   &
      phif ) bind(c,name='MG_interpolate')

    implicit none

    integer(c_int), intent(in)     :: Nc(3)

    integer(c_int), intent(in)     :: bLc(3)
    integer(c_int), intent(in)     :: bUc(3)

    integer(c_int), intent(in)     :: Nf(3)

    integer(c_int), intent(in)     :: bLf(3)
    integer(c_int), intent(in)     :: bUf(3)

    integer(c_int), intent(in)     :: iimax(3)
    integer(c_int), intent(in)     :: dd(3)

    real(c_double),  intent(in)    :: cI1 ( 1:2, 1:Nc(1) )
    real(c_double),  intent(in)    :: cI2 ( 1:2, 1:Nc(2) )
    real(c_double),  intent(in)    :: cI3 ( 1:2, 1:Nc(3) )

    real(c_double),  intent(in) :: phic (bLc(1):(Nc(1)+bUc(1)),bLc(2):(Nc(2)+bUc(2)),bLc(3):(Nc(3)+bUc(3)))

    real(c_double),  intent(out)   :: phif (bLf(1):(Nf(1)+bUf(1)),bLf(2):(Nf(2)+bUf(2)),bLf(3):(Nf(3)+bUf(3)))

    integer(c_int)               ::  i,ic
    integer(c_int)               ::  j
    integer(c_int)               ::  k

    !************************************************************************************
    !pgi$ unroll = n:8
    phif( 1:Nf(1):dd(1), 1:Nf(2):dd(2), 1:Nf(3):dd(3) ) = phic( 1:iimax(1), 1:iimax(2), 1:iimax(3) )
    !************************************************************************************
    !====================================================================================
    if( dd(3) /= 1 ) then ! (dimens == 2) <==> (dd(3) == 1) automatisch erfüllt!

      do k = 1+1, Nf(3)-1, dd(3)
        ic = k/dd(3)
        do j = 1, Nf(2), dd(2)
          !pgi$ unroll = n:8
          do i = 1, Nf(1), dd(1)
            phif(i,j,k) = cI3(1,ic)*phif(i,j,k-1) + cI3(2,ic)*phif(i,j,k+1)
          end do
        end do
      end do

    end if
    !====================================================================================
    if( dd(2) /= 1 ) then ! TEST!!! in 2D wird hier doppelte Arbeit geleistet! (NNf(3) == 2??)

      do k = 1, Nf(3)
        do j = 1+1, Nf(2)-1, dd(2)
          ic = j/dd(2)
          !pgi$ unroll = n:8
          do i = 1, Nf(1), dd(1)
            phif(i,j,k) = cI2(1,ic)*phif(i,j-1,k) + cI2(2,ic)*phif(i,j+1,k)
          end do
        end do
      end do

    end if
    !====================================================================================
    if( dd(1) /= 1 ) then

      do k = 1, Nf(3)
        do j = 1, Nf(2)
          !pgi$ unroll = n:8
          do i = 1+1, Nf(1)-1, dd(1)
            ic = i/dd(1)
            phif(i,j,k) = cI1(1,ic)*phif(i-1,j,k) + cI1(2,ic)*phif(i+1,j,k)
          end do
        end do
      end do

    end if
    !====================================================================================

  end subroutine MG_interpolate

  !> \brief interpolating velocities
  subroutine MG_interpolateV( &
      dir,                    &
      Nc,                     &
      bLc,bUc,                &
      !SSc,NNc,                &
      Nf,                     &
      bLf,bUf,                &
      !SSf,NNf,                &
      iimax,                  &
      dd,                     &
      cIV,                    &
      cI1,                    &
      cI2,                    &
      cI3,                    &
      phic,                   &
      phif ) bind (c,name='MG_interpolateV')

    implicit none

    integer(c_int), intent(in)     :: dir

    integer(c_int), intent(in)     :: Nc(1:3)

    integer(c_int), intent(in)     :: bLc(1:3)
    integer(c_int), intent(in)     :: bUc(1:3)

    !integer(c_int), intent(in)     :: SSc(1:3)
    !integer(c_int), intent(in)     :: NNc(1:3)

    integer(c_int), intent(in)     :: Nf(1:3)

    integer(c_int), intent(in)     :: bLf(1:3)
    integer(c_int), intent(in)     :: bUf(1:3)

    !integer(c_int), intent(in)     :: SSf(1:3)
    !integer(c_int), intent(in)     :: NNf(1:3)

    !integer(c_int), intent(in)     :: BCL(1:3)
    !integer(c_int), intent(in)     :: BCU(1:3)

    integer(c_int), intent(in)     :: iimax(1:3)
    integer(c_int), intent(in)     :: dd(1:3)

    real(c_double),  intent(in)    :: cIV ( 1:2, 0:Nf(dir) )

    real(c_double),  intent(in)    :: cI1 ( 1:2, 1:Nc(1) )
    real(c_double),  intent(in)    :: cI2 ( 1:2, 1:Nc(2) )
    real(c_double),  intent(in)    :: cI3 ( 1:2, 1:Nc(3) )

    real(c_double),  intent(in)   :: phic (bLc(1):(Nc(1)+bUc(1)),bLc(2):(Nc(2)+bUc(2)),bLc(3):(Nc(3)+bUc(3)))

    real(c_double),  intent(out)  :: phif (bLf(1):(Nf(1)+bUf(1)),bLf(2):(Nf(2)+bUf(2)),bLf(3):(Nf(3)+bUf(3)))

    integer(c_int)                ::  i, ic
    integer(c_int)                ::  j, jc
    integer(c_int)                ::  k, kc

    !integer(c_int)                :: l

    !integer(c_int)                :: S(1:3)
    !integer(c_int)                :: N(1:3)
    integer(c_int)                :: ds(1:3)

    ds = dd

    do i = 1,3
      if( iimax(i)<Nc(i) ) ds = 1
    end do

    if( 1==dir ) then

      if( dd(1) /= 1 ) then

        do k = 1, Nf(3), dd(3)
          if( 1==dd(3) ) then
            kc = k
          else
            kc = ( k+1 )/dd(3)
          end if
          do j = 1, Nf(2), dd(2)
            if( 1==dd(2) ) then
              jc = j
            else
              jc = ( j+1 )/dd(2) ! holy shit
            end if
            do i = 0, Nf(1) ! zero for dirichlet
              ic = ( i )/dd(1)
              phif(i,j,k) = cIV(1,i)*phic(ic,jc,kc) + cIV(2,i)*phic(ic+1,jc,kc)
            end do
          end do
        end do

      else

        do k = 1, Nf(3), dd(3)
          if( 1==dd(3) ) then
            kc = k
          else
            kc = ( k+1 )/dd(3)
          end if
          do j = 1, Nf(2), dd(2)
            if( 1==dd(2) ) then
              jc = j
            else
              jc = ( j+1 )/dd(2) ! holy shit
            end if
            do i = 0, Nf(1)
              phif(i,j,k) = phic(i,jc,kc)
            end do
          end do
        end do

      end if

      if( dd(2) /= 1 ) then

        do k = 1, Nf(3), dd(3)
          do j = 2, Nf(2)-1, dd(2)
            jc = ( j )/dd(2)
            !pgi$ unroll = n:8
            do i = 0, Nf(1)
              phif(i,j,k) = cI2(1,jc)*phif(i,j-1,k) + cI2(2,jc)*phif(i,j+1,k)
            end do
          end do
        end do

      end if

      if( dd(3) /= 1 ) then

        do k = 2, Nf(3)-1, dd(3)
          kc = k/dd(3)
          do j = 1, Nf(2)
            !pgi$ unroll = n:8
            do i = 0, Nf(1)
              phif(i,j,k) = cI3(1,kc)*phif(i,j,k-1) + cI3(2,kc)*phif(i,j,k+1)
            end do
          end do
        end do

      end if

    end if

    if( 2==dir ) then

      if( dd(2) /= 1 ) then

        do k = 1, Nf(3), dd(3)
          if( 1==dd(3) ) then
            kc = k
          else
            kc = ( k+1 )/dd(3)
          end if
          do j = 0, Nf(2) ! zero for dirichlet
            jc = ( j )/dd(2)
            do i = 1, Nf(1), dd(1)
              if( 1==dd(1) ) then
                ic = i
              else
                ic = ( i+1 )/dd(1) ! holy shit
              end if
              phif(i,j,k) = cIV(1,j)*phic(ic,jc,kc)+cIV(2,j)*phic(ic,jc+1,kc)
            end do
          end do
        end do

      else

        do k = 1, Nf(3), dd(3)
          if( 1==dd(3) ) then
            kc = k
          else
            kc = ( k+1 )/dd(3)
          end if
          do j = 0, Nf(2)
            do i = 1, Nf(1), dd(1)
              if( 1==dd(1) ) then
                ic = i
              else
                ic = ( i+1 )/dd(1) ! holy shit
              end if
              phif(i,j,k) = phic(ic,j,kc)
            end do
          end do
        end do

      end if

      if( dd(1) /= 1 ) then

        do k = 1, Nf(3), dd(3)
          do j = 0, Nf(2)
            !pgi$ unroll = n:8
            do i = 2, Nf(1)-1, dd(1)
              ic = i/dd(1)
              phif(i,j,k) = cI1(1,ic)*phif(i-1,j,k) + cI1(2,ic)*phif(i+1,j,k)
            end do
          end do
        end do

      end if

      if( dd(3) /= 1 ) then

        do k = 2, Nf(3)-1, dd(3)
          kc = k/dd(3)
          do j = 0, Nf(2)
            !pgi$ unroll = n:8
            do i = 1, Nf(1)
              phif(i,j,k) = cI3(1,kc)*phif(i,j,k-1) + cI3(2,kc)*phif(i,j,k+1)
            end do
          end do
        end do

      end if

    end if

    if( 3==dir ) then

      if( dd(3) /= 1 ) then

        do k = 0, Nf(3)
          kc = ( k )/dd(3)
          do j = 1, Nf(2), dd(2) ! zero for dirichlet
            if( 1==dd(2) ) then
              jc = j
            else
              jc = ( j+1 )/dd(2)
            end if
            do i = 1, Nf(1), dd(1)
              if( 1==dd(1) ) then
                ic = i
              else
                ic = ( i+1 )/dd(1) ! holy shit
              end if
              phif(i,j,k) = cIV(1,k)*phic(ic,jc,kc)+cIV(2,k)*phic(ic,jc,kc+1)
            end do
          end do
        end do

      else

        do k = 0, Nf(3)
          do j = 1, Nf(2), dd(2)
            if( 1==dd(2) ) then
              jc = j
            else
              jc = ( j+1 )/dd(2)
            end if
            do i = 1, Nf(1), dd(1)
              if( 1==dd(1) ) then
                ic = i
              else
                ic = ( i+1 )/dd(1) ! holy shit
              end if
              phif(i,j,k) = phic(ic,jc,k)
            end do
          end do
        end do

      end if

      if( dd(1) /= 1 ) then

        do k = 0, Nf(3)
          do j = 1, Nf(2), dd(2)
            !pgi$ unroll = n:8
            do i = 2, Nf(1)-1, dd(1)
              ic = i/dd(1)
              phif(i,j,k) = cI1(1,ic)*phif(i-1,j,k) + cI1(2,ic)*phif(i+1,j,k)
            end do
          end do
        end do

      end if

      if( dd(2) /= 1 ) then

        do k = 0, Nf(3)
          do j = 2, Nf(2)-1, dd(2)
            jc = j/dd(2)
            !pgi$ unroll = n:8
            do i = 1, Nf(1)
              phif(i,j,k) = cI2(1,jc)*phif(i,j-1,k) + cI2(2,jc)*phif(i,j+1,k)
            end do
          end do
        end do

      end if

    end if

  end subroutine MG_interpolateV

  !> \brief
  !! \note todo understand recevbuf size not correct
  subroutine MG_InterpolateScatter( &
      Nc,                           &
      bLc,                          &
      bUc,                          &
      !np,                           &
      iimax,                        &
      n_gather,                     &
      participate_yes,              &
      rankc2,                       &
      comm2,                        &
      dispI,                        &
      offsI,                        &
      phic ) bind(c,name='MG_InterpolateScatter')

    implicit none

    integer(c_int), intent(in)     :: Nc(3)

    integer(c_int), intent(in)     :: bLc(3)
    integer(c_int), intent(in)     :: bUc(3)

    !integer(c_int), intent(in)     :: np(3)

    integer(c_int), intent(in)     :: iimax(3)

    integer(c_int), intent(in)     :: n_gather(3)

    logical(c_bool), intent(in)    :: participate_yes

    integer(c_int), intent(in)     :: rankc2

    integer(c_int), intent(in)     :: comm2

    integer(c_int), intent(in)     :: dispI(     1:n_gather(1)*n_gather(2)*n_gather(3) )
    integer(c_int), intent(in)     :: offsI(1:3, 1:n_gather(1)*n_gather(2)*n_gather(3) )

    real(c_double),intent(inout)   :: phic (bLc(1):(Nc(1)+bUc(1)),bLc(2):(Nc(2)+bUc(2)),bLc(3):(Nc(3)+bUc(3)))

    integer(c_int)                ::  i, ii
    integer(c_int)                ::  j, jj
    integer(c_int)                ::  k, kk

    integer(c_int)                ::  merror
    real(c_double), allocatable   ::  recvbuf(:,:,:)
    integer(c_int)                ::  offsg(1:3), dispg
    !real(c_double)                ::  sendbuf( 1:( ( Nc(1)+np(1)-1 )*( Nc(2)+np(2)-1 )*( Nc(3)+np(3)-1 ) ) ) ! from IMPACT
    !real(c_double)                ::  sendbuf( 0:( ( Nc(1)+np(1)*3 )*( Nc(2)+np(2)*3 )*( Nc(3)+np(3)*3 ) ) )
    real(c_double)                ::  sendbuf( 0:( (iimax(1)+1)*(iimax(2)+1)*(iimax(3)+1)*n_gather(1)*n_gather(2)*n_gather(3) ) )
    !real(c_double)                ::  sendbuf( 0:( ( Nc(1)+12 )*( Nc(2)+12 )*( Nc(3)+12 ) ) )

    if( participate_yes ) then

      do k = 1, n_gather(3)
        do j = 1, n_gather(2)
          do i = 1, n_gather(1)

            dispg      = dispI(    i+(j-1)*n_gather(1)+(k-1)*n_gather(1)*n_gather(2))
            offsg(1:3) = offsI(1:3,i+(j-1)*n_gather(1)+(k-1)*n_gather(1)*n_gather(2))

            do kk = 0, iimax(3)
              do jj = 0, iimax(2)
                do ii = 0, iimax(1)
                  sendbuf( dispg + ii + jj*(iimax(1)+1) + kk*(iimax(1)+1)*(iimax(2)+1) ) = phic(ii+offsg(1),jj+offsg(2),kk+offsg(3))
                end do
              end do
            end do

          end do
        end do
      end do

    end if

    allocate(recvbuf(0:iimax(1),0:iimax(2),0:iimax(3))) ! Anmerkung: Besser nicht fest allocieren um Speicherplatz zu sparen, ODER gleich "phic" verwenden!

    call MPI_SCATTER(sendbuf,(iimax(1)+1)*(iimax(2)+1)*(iimax(3)+1),MPI_REAL8,recvbuf,(iimax(1)+1)*(iimax(2)+1)*(iimax(3)+1),MPI_REAL8,rankc2,comm2,merror)

    phic(0:iimax(1),0:iimax(2),0:iimax(3)) = recvbuf(0:iimax(1),0:iimax(2),0:iimax(3))

    deallocate(recvbuf)
    !------------------------------------------------------------------------------

  end subroutine MG_InterpolateScatter

end module cmod_InterpolationOp

max_stars_repo_name: huppd/PINTimpact

max_issues_repo_name: huppd/PINTimpact

max_forks_repo_name: huppd/PINTimpact


text:
      MODULE itjj_I
      INTERFACE
!
      INTEGER FUNCTION ITJJ(IK,ID,KG,BK,BD,IBT,BT,KG1,ITP,ITG,IQ)
      USE vast_kind_param, ONLY:  DOUBLE
      INTEGER,      INTENT(IN)                :: KG, KG1, IQ
      INTEGER,      INTENT(OUT)               :: ITP, ITG
      INTEGER,      INTENT(IN),  DIMENSION(7) :: IK, ID
      INTEGER,      INTENT(OUT), DIMENSION(7) :: IBT
      REAL(DOUBLE), INTENT(IN),  DIMENSION(3) :: BK, BD
      REAL(DOUBLE), INTENT(OUT), DIMENSION(3) :: BT
      END FUNCTION
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


! { dg-do compile }
!
! Error checking for the SELECT TYPE statement
!
! Contributed by Janus Weil <janus@gcc.gnu.org>

  type :: t1
    integer :: i = 42
    class(t1),pointer :: cp
  end type

  type, extends(t1) :: t2
    integer :: j = 99
  end type

  type :: t3
    real :: r
  end type

  type :: ts
    sequence
    integer :: k = 5
  end type

  class(t1), pointer :: a => NULL()
  type(t1), target :: b
  type(t2), target :: c
  a => b
  print *, a%i

  type is (t1)  ! { dg-error "Unexpected TYPE IS statement" }

  select type (3.5)  ! { dg-error "is not a named variable" }
  select type (a%cp) ! { dg-error "is not a named variable" }
  select type (b)    ! { dg-error "Selector shall be polymorphic" }
  end select

  select type (a)
    print *,"hello world!"  ! { dg-error "Expected TYPE IS, CLASS IS or END SELECT" }
  type is (t1)
    print *,"a is TYPE(t1)"
  type is (t2)
    print *,"a is TYPE(t2)"
  class is (ts)  ! { dg-error "must be extensible" }
    print *,"a is TYPE(ts)"
  type is (t3)   ! { dg-error "must be an extension of" }
    print *,"a is TYPE(t3)"
  type is (t4)   ! { dg-error "error in TYPE IS specification" }
    print *,"a is TYPE(t3)"
  class is (t1)
    print *,"a is CLASS(t1)"
  class is (t2) label  ! { dg-error "Syntax error" }
    print *,"a is CLASS(t2)"
  class default  ! { dg-error "cannot be followed by a second DEFAULT CASE" }
    print *,"default"
  class default  ! { dg-error "cannot be followed by a second DEFAULT CASE" }
    print *,"default2"
  end select

label: select type (a)
  type is (t1) label
    print *,"a is TYPE(t1)"
  type is (t2)  ! { dg-error "overlaps with TYPE IS" }
    print *,"a is TYPE(t2)"
  type is (t2)  ! { dg-error "overlaps with TYPE IS" }
    print *,"a is still TYPE(t2)"
  class is (t1) labe   ! { dg-error "Expected block name" }
    print *,"a is CLASS(t1)"
  end select label

end


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program 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(version 2);
!
!
!     This program 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 this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine midexternalfaces(iexternfa,nexternfa,ifacext,
     &   ifreefa,itetfa,ifac,kontet,kontetor,ialsetexternel,
     &   nexternel,iedgextfa,ifacexted,ipoed,iedg,iexternedg)
!
!     stores the nodes belonging to the external faces of the
!     unrefined mesh
!
      implicit none
!
      integer nexternfa,iexternfa(*),ifacext(6,*),i,j,ifreefa,
     &  iel,itetfa(2,*),ifac(4,*),kontetor(6,*),nofa(6,4),
     &  kontet(4,*),node,k,ialsetexternel(*),nexternel,id,
     &  iedgextfa(2,*),ifacexted(3,*),n1,n2,n3,indexe,iedg(3,*),
     &  iexternedg(*),ipoed(*),iedge
!
      intent(in) ifreefa,itetfa,ifac,kontet,kontetor,ipoed,iedg,
     &  iexternedg
!
      intent(inout) nexternfa,iexternfa,ifacext,ialsetexternel,
     &  nexternel,iedgextfa,ifacexted
!
!     local node numbers of the nodes belonging to the 4 faces
!     of a tetrahedron. From the middle nodes 4 is subtracted,
!     since they are stored in position 1..6 in field kontetor
!
      data nofa /2,3,4,2,6,5,
     &           3,1,4,3,4,6,
     &           1,2,4,1,5,4,
     &           1,3,2,3,2,1/
!
      nexternfa=0
      nexternel=0
!
!     loop over all faces. For the unrefined mesh the faces are stored
!     in ifac in a consecutive order.
!     A face i is external if iexternfa(i)!=0; at exit the value
!     of iexternfa(i) is the number of the external face in field
!     ifacext
!
      do i=1,ifreefa-1
         if(iexternfa(i).ne.0) then
            nexternfa=nexternfa+1
            iexternfa(i)=nexternfa
!
!           recovering the element containing the face
!
            iel=itetfa(1,i)
!
!           the local number of the face is the local number of the
!           vertex node not belonging to the face
!
            loop: do j=1,4
               node=kontet(j,iel)
               do k=1,3
                  if(ifac(k,i).eq.node) cycle loop
               enddo
               exit
            enddo loop
!
            do k=1,3
               ifacext(k,nexternfa)=kontet(nofa(k,j),iel)
               ifacext(k+3,nexternfa)=kontetor(nofa(k+3,j),iel)
            enddo
!
!           catalogueing the elements adjacent to an external face =
!           external elements
!
!           these elements are needed for the projection on the
!           external surface for newly generated nodes in the
!           refined mesh
!
c            call nident(ialsetexternel,iel,nexternel,id)
c            if(id.gt.0) then
c               if(ialsetexternel(id).eq.iel) cycle
c            endif
c            nexternel=nexternel+1
c            do j=nexternel,id+2,-1
c               ialsetexternel(j)=ialsetexternel(j-1)
c            enddo
c            ialsetexternel(id+1)=iel
!
!           catalogueing the (external) edges per external face in
!           ifacexted and the external faces per external edge in
!           iedgextfa
!
!           first (external) edge of external face
!
            n1=ifacext(1,nexternfa)
            n2=ifacext(2,nexternfa)
            if(n2.lt.n1) then
               n3=n1
               n1=n2
               n2=n3
            endif
c            write(*,*) 'midexternafaces1 ',n1,n2,nexternfa
            indexe=ipoed(n1)
            do
               if(iedg(2,indexe).eq.n2) then
                  iedge=iexternedg(indexe)
                  ifacexted(1,nexternfa)=iedge
                  if(iedgextfa(1,iedge).eq.0) then
                     iedgextfa(1,iedge)=nexternfa
                  else
                     iedgextfa(2,iedge)=nexternfa
                  endif
c                  write(*,*) iedge,iedgextfa(1,iedge),iedgextfa(2,iedge)
                  exit
               else
                  indexe=iedg(3,indexe)
               endif
            enddo
!
!           second (external) edge of external face
!
            n1=ifacext(2,nexternfa)
            n2=ifacext(3,nexternfa)
            if(n2.lt.n1) then
               n3=n1
               n1=n2
               n2=n3
            endif
c            write(*,*) 'midexternafaces2 ',n1,n2,nexternfa
            indexe=ipoed(n1)
            do
               if(iedg(2,indexe).eq.n2) then
                  iedge=iexternedg(indexe)
                  ifacexted(2,nexternfa)=iedge
                  if(iedgextfa(1,iedge).eq.0) then
                     iedgextfa(1,iedge)=nexternfa
                  else
                     iedgextfa(2,iedge)=nexternfa
                  endif
c                  write(*,*) iedge,iedgextfa(1,iedge),iedgextfa(2,iedge)
                  exit
               else
                  indexe=iedg(3,indexe)
               endif
            enddo
!
!           third (external) edge of external face
!
            n1=ifacext(3,nexternfa)
            n2=ifacext(1,nexternfa)
            if(n2.lt.n1) then
               n3=n1
               n1=n2
               n2=n3
            endif
c            write(*,*) 'midexternafaces3 ',n1,n2,nexternfa
            indexe=ipoed(n1)
            do
               if(iedg(2,indexe).eq.n2) then
                  iedge=iexternedg(indexe)
                  ifacexted(3,nexternfa)=iedge
                  if(iedgextfa(1,iedge).eq.0) then
                     iedgextfa(1,iedge)=nexternfa
                  else
                     iedgextfa(2,iedge)=nexternfa
                  endif
c                  write(*,*) iedge,iedgextfa(1,iedge),iedgextfa(2,iedge)
                  exit
               else
                  indexe=iedg(3,indexe)
               endif
            enddo
!
         endif
      enddo
!
c      do i=1,54
c         write(*,*) 'iedgext/fa',
c     &            (iedgextfa(j,i),j=1,2)
c      enddo
c      do i=1,36
c         write(*,*) 'ifacext/ed',(ifacext(j,i),j=1,6),
c     &              (ifacexted(j,i),j=1,3)
c      enddo
c      do i=1,69
c         write(*,*) 'iedg',(iedg(j,i),j=1,2),iexternedg(i)
c      enddo
c      do i=1,80
c         write(*,*) 'ifac',(ifac(j,i),j=1,3),iexternfa(i)
c      enddo
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop