text:
!> Orchid -- 2D/3D Euler/MagnetoHydroDynamics solver.
!> Copyright (C) Butakov Oleg 2019.

Module orchid_solver_grid_gauss_const
Implicit None
Contains
!########################################################################################################
!########################################################################################################
!########################################################################################################
Function mhd_grid_gauss_point(Npoints, k) Result(x)
    !> Points of the Gauss quadrature.
    !> {{{
    Integer, Intent(In) :: Npoints
    Integer, Intent(In) :: k
    Real(8) :: x
    !> }}}
    Real(8), Dimension(1:1), Parameter :: x1 = [ +0.000000000D0 ]
    Real(8), Dimension(1:2), Parameter :: x2 = [ -0.577350269D0, +0.577350269D0 ]
    Real(8), Dimension(1:3), Parameter :: x3 = [ -0.774596669D0, +0.000000000D0, +0.774596669D0 ]
    Real(8), Dimension(1:4), Parameter :: x4 = [ -0.861136312D0, -0.339981044D0, +0.339981044D0, &
                                                 +0.861136312D0 ]
    Real(8), Dimension(1:5), Parameter :: x5 = [ -0.906179846D0, -0.538469310D0, +0.000000000D0, &
                                                 +0.538469310D0, +0.906179846D0 ]
    Select Case ( Npoints )
        Case ( 1 )
            x = x1(k)
        Case ( 2 )
            x = x2(k)
        Case ( 3 )
            x = x3(k)
        Case ( 4 )
            x = x4(k)
        Case ( 5 )
            x = x5(k)
        Case Default
            Write (0,*) 'Too many points for Gauss quadratures.'
            Error Stop -1
    End Select
End Function mhd_grid_gauss_point
!########################################################################################################
!########################################################################################################
!########################################################################################################
Function mhd_grid_gauss_weight(Npoints, k) Result(w)
    !> Weights of the Gauss quadrature.
    !> {{{
    Integer, Intent(In) :: Npoints
    Integer, Intent(In) :: k
    Real(8) :: w
    !> }}}
    Real(8), Dimension(1:1), Parameter :: w1 = [ +1.000000000D0 ]
    Real(8), Dimension(1:2), Parameter :: w2 = [ +0.500000000D0, +0.50000000D0 ]
    Real(8), Dimension(1:3), Parameter :: w3 = [ +0.277777778D0, +0.444444444D0, +0.277777778D0 ]
    Real(8), Dimension(1:4), Parameter :: w4 = [ +0.173927423D0, +0.326072577D0, +0.326072577D0, &
                                                 +0.173927423D0 ]
    Real(8), Dimension(1:5), Parameter :: w5 = [ +0.118463443D0, +0.239314335D0, +0.284444444D0, &
                                                 +0.239314335D0, +0.118463443D0 ]
    Select Case ( Npoints )
        Case ( 1 )
            w = w1(k)
        Case ( 2 )
            w = w2(k)
        Case ( 3 )
            w = w3(k)
        Case ( 4 )
            w = w4(k)
        Case ( 5 )
            w = w5(k)
        Case Default
            Write (0,*) 'Too many points for Gauss quadratures.'
            Error Stop -1
    End Select
End Function mhd_grid_gauss_weight
!########################################################################################################
!########################################################################################################
!########################################################################################################
End Module orchid_solver_grid_gauss_const

Module orchid_solver_grid_gauss
Use orchid_solver_params
Use orchid_solver_grid
Use orchid_solver_grid_gauss_const
Implicit None
Type :: MhdGridGaussNode
    Real(8) :: x, y, z, w
    Contains
End Type MhdGridGaussNode
Type :: MhdGridGaussCell
    Integer :: nnode, nnode_end
    Contains
End Type MhdGridGaussCell
Type :: MhdGridGaussFace
    Integer :: nnode, nnode_end
    Contains
End Type MhdGridGaussFace
Type, Extends(MhdGrid) :: MhdGridGauss
    Integer :: ncell_nodes_min, ncell_nodes_max
    Integer :: nface_nodes_min, nface_nodes_max
    Type(MhdGridGaussCell), Dimension(:), Allocatable :: gcells
    Type(MhdGridGaussFace), Dimension(:), Allocatable :: gfaces
    Type(MhdGridGaussNode), Dimension(:), Allocatable :: cell_nodes, face_nodes
    Contains
    Procedure, Public, Non_Overridable :: init_gauss => mhd_grid_gauss_init
    Procedure, Public, Non_Overridable :: init_gauss1D => mhd_grid_gauss_init1D
    Procedure, Public, Non_Overridable :: init_gauss_dummy => mhd_grid_gauss_init_dummy
End Type MhdGridGauss
Private :: mhd_grid_gauss_init1D, &
           mhd_grid_gauss_init_dummy
Contains
!########################################################################################################
!########################################################################################################
!########################################################################################################
Subroutine mhd_grid_gauss_init(This, Npoints)
    !> Initialize a Gauss grid.
    !> {{{
    Class(MhdGridGauss), Intent(InOut) :: This
    Integer, Intent(In), Optional :: Npoints
    !> }}}
    If ( Present(Npoints) ) Then
        Select Case ( This%ndims )
            Case ( 1 )
                Call This%init_gauss1D(Npoints)
            !Case ( 2 )
            !    Call This%init_gauss2D(Npoints)
            !Case ( 3 )
            !    Call This%init_gauss3D(Npoints)
            Case Default
                Write (0,*) 'Invalid amount of the grid dimensions. ', &
                            'Was grid initialized?'
                Error Stop -1
        End Select
    Else
        Call This%init_gauss_dummy()
    End If
End Subroutine mhd_grid_gauss_init
!########################################################################################################
!########################################################################################################
!########################################################################################################
Subroutine mhd_grid_gauss_init1D(This, Npoints)
    !> Initialize a 1D Gauss grid.
    !> {{{
    Class(MhdGridGauss), Intent(InOut) :: This
    Integer, Intent(In) :: Npoints
    !> }}}
    Integer :: n, k, kk
    Real(8) :: h, x
    !>-------------------------------------------------------------------------------
    !> Allocate the Gauss Grid.
    This%ncell_nodes_min = This%ncells_min
    This%ncell_nodes_max = This%ncells_max*Npoints
    This%nface_nodes_min = This%nfaces_min
    This%nface_nodes_max = This%nfaces_max
    Allocate(This%gcells(This%ncells_min:This%ncells_max))
    Allocate(This%gfaces(This%nfaces_min:This%nfaces_max))
    Allocate(This%cell_nodes(This%ncell_nodes_min:This%ncell_nodes_max))
    Allocate(This%face_nodes(This%nface_nodes_min:This%nface_nodes_max))
    !>-------------------------------------------------------------------------------

    !>-------------------------------------------------------------------------------
    !> Initialize Cells and CellNodes.
    !$OMP Parallel Do Private(k, kk, h, x)
    Do n = This%ncells_min, This%ncells_max
        h = This%cells(n)%Vcell
        x = This%nodes(This%cell2node(This%cells(n)%nnode))%x
        This%gcells(n)%nnode = Npoints*( n - 1 ) + 1
        This%gcells(n)%nnode_end = Npoints*n
        Do k = This%gcells(n)%nnode, This%gcells(n)%nnode_end
            kk = k - This%gcells(n)%nnode + 1
            This%cell_nodes(k)%w = mhd_grid_gauss_weight(Npoints, kk)
            This%cell_nodes(k)%x = mhd_grid_gauss_point(Npoints, kk)
            This%cell_nodes(k)%x = x + 0.5D0*h*( This%cell_nodes(k)%x + 1.0D0 )
            This%cell_nodes(k)%y = 0.0D0
            This%cell_nodes(k)%z = 0.0D0
        End Do
    End Do
    !$OMP End Parallel Do
    !> Initialize Faces and FaceNodes.
    !$OMP Parallel Do
    Do n = This%nfaces_min, This%nfaces_max
        This%gfaces(n)%nnode = n
        This%gfaces(n)%nnode_end = n
        This%face_nodes(n)%x = This%faces(n)%x
        This%face_nodes(n)%y = This%faces(n)%y
        This%face_nodes(n)%z = This%faces(n)%z
        This%face_nodes(n)%w = 1.0D0
    End Do
    !$OMP End Parallel Do
    !>-------------------------------------------------------------------------------
End Subroutine mhd_grid_gauss_init1D
!########################################################################################################
!########################################################################################################
!########################################################################################################
Subroutine mhd_grid_gauss_init_dummy(This)
    !> Initialize a dummy Gauss grid (for first-order integration).
    !> {{{
    Class(MhdGridGauss), Intent(InOut) :: This
    !> }}}
    Integer :: n
    !>-------------------------------------------------------------------------------
    !> Allocate the Gauss Grid.
    This%ncell_nodes_min = This%ncells_min
    This%ncell_nodes_max = This%ncells_max
    This%nface_nodes_min = This%nfaces_min
    This%nface_nodes_max = This%nfaces_max
    Allocate(This%gcells(This%ncells_min:This%ncells_max))
    Allocate(This%gfaces(This%nfaces_min:This%nfaces_max))
    Allocate(This%cell_nodes(This%ncell_nodes_min:This%ncell_nodes_max))
    Allocate(This%face_nodes(This%nface_nodes_min:This%nface_nodes_max))
    !>-------------------------------------------------------------------------------

    !>-------------------------------------------------------------------------------
    !> Initialize Cells and CellNodes.
    !$OMP Parallel Do
    Do n = This%ncells_min, This%ncells_max
        This%gcells(n)%nnode = n
        This%gcells(n)%nnode_end = n
        This%cell_nodes(n)%x = This%cells(n)%x
        This%cell_nodes(n)%y = This%cells(n)%y
        This%cell_nodes(n)%z = This%cells(n)%z
        This%cell_nodes(n)%w = 1.0D0
    End Do
    !$OMP End Parallel Do
    !> Initialize Faces and FaceNodes.
    !$OMP Parallel Do
    Do n = This%nfaces_min, This%nfaces_max
        This%gfaces(n)%nnode = n
        This%gfaces(n)%nnode_end = n
        This%face_nodes(n)%x = This%faces(n)%x
        This%face_nodes(n)%y = This%faces(n)%y
        This%face_nodes(n)%z = This%faces(n)%z
        This%face_nodes(n)%w = 1.0D0
    End Do
    !$OMP End Parallel Do
    !>-------------------------------------------------------------------------------
End Subroutine mhd_grid_gauss_init_dummy
!########################################################################################################
!########################################################################################################
!########################################################################################################
End Module orchid_solver_grid_gauss

max_stars_repo_name: Jhuighuy/OrchidSolver

max_issues_repo_name: Jhuighuy/OrchidSolver

max_forks_repo_name: Jhuighuy/OrchidSolver


text:
Hyun Do Ahn was a South Korean exchange student here at UC Davis who died during an 11day backpacking trip in Yosemite National Park. He studied English and learned how to surf in the States.

He was reported missing when he missed his flight back home to South Korea.

He died of an apparent fall 20 miles from the beginning of his trek. A ranger said that of the last five hikers to go missing in Yosemite, four (including Ahn) traveled alone. http://www.dailyrepublic.com/articles/2005/08/11/local_news/news07.txt Its not something we encourage, the ranger said.

I was actually a friend of Hyun Do. I know it was not his original intention of going on this expedition alone, but decided to continue it even after not finding anyone to go with him. According to Yosemite officials he may have sustained a fall, but he was cause of death was due to multiple blunt force injuries most likely sustained after slipping while attempting to cross a swollen creek. According to journal entries, he most likely died 5 days after his missed flight.  Users/CarlosBarahona

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
These are unicode dice for the dice page.  You can add a rolled die to any entry by using {{{RandomQuote(Dice/d6)}}}

  Image(d601.png, noborder)
  Image(d602.png, noborder)
  Image(d603.png, noborder)
  Image(d604.png, noborder)
  Image(d605.png, noborder)
  Image(d606.png, noborder)


text:
C ======================================================================
C  Diffusion tensor
C ======================================================================
      Integer Function Ddiff(x, y, z, label, DDATA, IDATA, iSYS, Coef)
      include 'fem3Dtet.fd'
C ======================================================================
C  Diffusion tensor is passed through DDATA array. DDATA is filled
C  in main.f using subroutine Drotate at the end of this file.
C ======================================================================
      Real*8  x, y, z, DDATA(3,*), Coef(9, 9)
      Integer label, IDATA(*), iSYS(*)
C ======================================================================
      iSYS(1) = 3
      iSYS(2) = 3

c HMFEM needs inverted diffusion tensor; we use Lapack for this
      Do i = 1, 3
         Do j = i, 3
            Coef(i, j) = DDATA(i, j)
         End do
      End do

      Call invertSPDmatrix(3, Coef, 9)

      Ddiff = TENSOR_SYMMETRIC

      Return
      End

C ======================================================================
C Boundary condition
C ======================================================================
      Integer Function Dbc(x, y, z, label, DDATA, IDATA, iSYS, eBC)
      Include 'assemble.fd'
C ======================================================================
      Real*8  x, y, z, DDATA(*), eBC(*)
      Integer label, IDATA(*), iSYS(*)
C ======================================================================
      iSYS(1) = 1
      iSYS(2) = 1

      If (label.ge.1.and.label.le.6) then
          eBC(1) = 0D0
          Dbc = BC_DIRICHLET
      Else if (label.ge.7.and.label.le.12) then
          eBC(1) = 2D0
          Dbc = BC_DIRICHLET
      Else
          stop 'wrong labelF'
      End if

      Return
      End

C ======================================================================
C Right hand side
C ======================================================================
      Integer Function Drhs(x, y, z, label, DDATA, IDATA, iSYS, F)
      Include 'fem3Dtet.fd'
C ======================================================================
      Real*8  x, y, z, DDATA(*), F(*)
      Integer label, IDATA(*), iSYS(*)
C ======================================================================
      iSYS(1) = 1
      iSYS(2) = 1

      F(1) = 0D0
      Drhs = TENSOR_SCALAR

      Return
      End

C ======================================================================
c Templated routine for elemental matrix
C ======================================================================
      Subroutine FEM3Dext(
     &           XY1, XY2, XY3, XY4,
     &           lbE, lbF, lbR, lbP, DDATA, IDATA, iSYS,
     &           LDA, A, F, nRow, nCol,
     &           templateR, templateC)
      Implicit none
      Include 'fem3Dtet.fd'
C ======================================================================
      Real*8  XY1(*), XY2(*), XY3(*), XY4(*)

      Integer lbE, lbF(4), lbR(6), lbP(4)
      Real*8  DDATA(*)
      Integer IDATA(*), iSYS(*)

      Integer LDA, nRow, nCol
      Real*8  A(LDA, *), F(*)
      Integer templateR(*), templateC(*)

C Local variables
      Integer  Ddiff, Drhs, Dbc
      External Ddiff, Drhs, Dbc
      Real*8   tri_area
      External tri_area

      Real*8   Q(4, 4), B(4), L(4), G(1), det, alpha, s
      Real*8   x, y, z, eBC(1)
      Integer  i, j, k, ir,ic, INFO
C ======================================================================
      nRow = 4
      nCol = 4

c ... set up templates (unknowns are located on faces)
      Do i = 1, 4
         templateR(i) = Fdof
         templateC(i) = Fdof
      End do

c ... compute the mass matrix (Q)
c     pass a copy of iSYS which is changed inside fem3Dtet
      Call fem3Dtet(XY1, XY2, XY3, XY4,
     &              IDEN, FEM_RT0, IDEN, FEM_RT0,
     &              lbE, Ddiff, DDATA, IDATA, iSYS, 2,
     &              4, Q, ir, ic)

c ... compute right hand side
c     pass a copy of iSYS which is changed inside fem3Dtet
      Call fem3Dtet(XY1, XY2, XY3, XY4,
     &              IDEN, FEM_P0, IDEN, FEM_P0,
     &              lbE, Drhs, DDATA, IDATA, iSYS, 2,
     &              1, G, ir, ic)

c ... compute the constraint matrix
      B(1) = tri_area(XY1, XY2, XY3)
      B(2) = tri_area(XY2, XY3, XY4)
      B(3) = tri_area(XY3, XY4, XY1)
      B(4) = tri_area(XY4, XY1, XY2)

c ... elliminate fluxes: invert the mass matrix first
      call invertSPDmatrix(4,Q,4)

      Do i = 1, 4
         Do j = 1, 4
            Q(i, j) = Q(i, j) * B(j)
         End do
      End do

      Do i = 1, 4
         L(i) = Q(1, i) * B(1) + Q(2, i) * B(2)
     &        + Q(3, i) * B(3) + Q(4, i) * B(4)
      End do

c ... elliminate central pressure
      alpha = L(1) + L(2) + L(3) + L(4)

      Do i = 1, 4
         Do j = i, 4
            A(i, j) = B(i) * Q(i, j) - L(i) * L(j) / alpha
            A(j, i) = A(i, j)
         End do
         F(i) = L(i) * G(1) / alpha
      End do

c ... impose boundary conditions (assume nRow = nCol)
      Do k = 1, 4
         If(lbF(k).GT.0) Then
            If(k.EQ.1) Then
               x = (XY1(1) + XY2(1) + XY3(1)) / 3
               y = (XY1(2) + XY2(2) + XY3(2)) / 3
               z = (XY1(3) + XY2(3) + XY3(3)) / 3
            ElseIf(k.EQ.2) Then
               x = (XY2(1) + XY3(1) + XY4(1)) / 3
               y = (XY2(2) + XY3(2) + XY4(2)) / 3
               z = (XY2(3) + XY3(3) + XY4(3)) / 3
            ElseIf(k.EQ.3) Then
               x = (XY3(1) + XY4(1) + XY1(1)) / 3
               y = (XY3(2) + XY4(2) + XY1(2)) / 3
               z = (XY3(3) + XY4(3) + XY1(3)) / 3
            ElseIf(k.EQ.4) Then
               x = (XY4(1) + XY1(1) + XY2(1)) / 3
               y = (XY4(2) + XY1(2) + XY2(2)) / 3
               z = (XY4(3) + XY1(3) + XY2(3)) / 3
            End if

c           pass of copy of iSYS which is changed inside Dbc
            i = Dbc(x, y, z, lbF(k), DDATA, IDATA, iSYS, eBC)

            Do i = 1, nRow
               F(i) = F(i) - A(i, k) * eBC(1)
               A(k, i) = 0
               A(i, k) = 0
            End do

            A(k, k) = 1
            F(k) = eBC(1)
         End if
      End do

      Return
      End

C ======================================================================
C  The routine takes an input matrix D and rotates it about i-th
C  coordinate axis by angle theta.
C ======================================================================
      Subroutine Drotate(D, i, theta)
C ======================================================================
      Real*8  D(3, 3), theta, c, s, a, b
C ======================================================================
      j = i + 1
      If( j.GT.3 ) j = 1

      c = dcos(theta)
      s = dsin(theta)

      Do k = 1, 3
         a = D(k, i)
         b = D(k, j)

         D(k, i) = a * c - b * s
         D(k, j) = a * s + b * c
      End do

      Do k = 1, 3
         a = D(i, k)
         b = D(j, k)

         D(i, k) = a * c - b * s
         D(j, k) = a * s + b * c
      End do

      Return
      End

c ==========================================================
      Subroutine HMFEMrecoverUP(nv, vrt, nt, tet, material, DDATA,IDATA,
     &                          SOL, SOLatTets, VECatTets, iW, MaxWi)
c ==========================================================
c  Recovers pressures SOLatTets at cell centers and velocities
c  VECatTets at cells centers from Lagrange multipliers SOL
c  defined on mesh faces
c ==========================================================
      Implicit none
      include 'fem3Dtet.fd'

      Integer  nv, nt, MaxWi
      Real*8   vrt(3,*), DDATA(*)
      Integer  tet(4,*), material(*), IDATA(*),  iW(*)
      Real*8   SOL(*), SOLatTets(*), VECatTets(3,*)

      EXTERNAL FEM3Dext, Ddiff, Drhs
      Integer  Ddiff, Drhs

      Real*8    tri_area, calVol
      External  tri_area, calVol

C=====LOCAL VARIABLES======
      Integer   n,i,j,k, i1,i2,i3,i4, ir,ic, ifc, info
      Integer   nf, iIFE, inEP, iIEP, iiW, iSYS(5)

      Real*8    localA(5, 5), localF(5), xyzc(3), rW(100), s
      Real*8    face_area(4), vol

      Integer   Face2OpPnt(4)
      Data      Face2OpPnt / 4, 1, 2, 3 /

c ==========================================================
c     make a list E -> F
      iIFE = 1
      inEP = iIFE + 4 * nt
      iIEP = inEP + nv
      iiW  = iIEP + 4 * nt
      If(iiW .gt.MaxWi) Stop 'Please increase MaxWi'

      Call listE2F(nv, nf, nt, tet, iW(iIFE), iW(inEP), iW(iIEP))

c     loop over mesh cells
      Do n = 1, nt
         i1 = tet(1, n)
         i2 = tet(2, n)
         i3 = tet(3, n)
         i4 = tet(4, n)

c        compute mass matrix using material properties
c        we use the user-specified routine Ddiff()
         iSYS(1) = n
         iSYS(2) = i1
         iSYS(3) = i2
         iSYS(4) = i3
         iSYS(5) = i4

         Call fem3Dtet(
     &        vrt(1,i1), vrt(1,i2), vrt(1,i3), vrt(1,i4),
     &        IDEN, FEM_RT0, IDEN, FEM_RT0,
     &        material(n), Ddiff, DDATA, IDATA, iSYS, 2,
     &        5, localA, ir, ic)

c        add mass conservation equation (constraint matrix)
c        we use the divergent theorem here as we did in forlibfem.f
         localA(1, 5) = tri_area(vrt(1,i1), vrt(1,i2), vrt(1,i3))
         localA(2, 5) = tri_area(vrt(1,i2), vrt(1,i3), vrt(1,i4))
         localA(3, 5) = tri_area(vrt(1,i3), vrt(1,i4), vrt(1,i1))
         localA(4, 5) = tri_area(vrt(1,i4), vrt(1,i1), vrt(1,i2))

         localA(5, 5) = 0D0

         Do i = 1, 4
            face_area(i) = localA(i, 5)
         End do

c        compute right hand side (Lagrange muptipliers are in SOL)
  !       if(n.eq.227) then

!         write(*,*)iW(iIFE + 4*(n-1) + 0),iW(iIFE + 4*(n-1) + 1),
!     & iW(iIFE + 4*(n-1) + 2),iW(iIFE + 4*(n-1) + 3)
 !        write(*,*) tet(1,n),tet(2,n),tet(3,n),tet(4,n)
 !       endif

         Do i = 1, 4
            ifc = iW(iIFE + 4*(n-1) + i-1)

            localF(i) = localA(i, 5) * SOL(ifc)
         End do

c        we use the user-specified routine Drhs()
         Call fem3Dtet(
     &        vrt(1,i1), vrt(1,i2), vrt(1,i3), vrt(1,i4),
     &        IDEN, FEM_P0, IDEN, FEM_P0,
     &        material(n), Drhs, DDATA, IDATA, iSYS, 2,
     &        1, localF(5), ir, ic)

c        solve the 5x5 system localA x = localF. The first four
c        entries in the solution x are velocities on tetrahedron faces.
c        The 5-th entry is the solution p at the mass center of the
c        tetrahedron. The user may use the LAPACK routine dsysv.

         Call dsysv('U', 5, 1, localA, 5, iW, localF, 5, rW, 100, info)
         If (info.ne.0) stop 'error in dsysv'

c        recover pressure
         SOLatTets(n) = localF(5)

c        recover velocity vector
         vol = calVol(vrt(1,i1), vrt(1,i2), vrt(1,i3), vrt(1,i4))
         vol = dabs(vol)

         Do i = 1, 3
            xyzc(i) = (vrt(i,i1) + vrt(i,i2) + vrt(i,i3) + vrt(i,i4))/4
            VECatTets(i,n) = 0d0
         End do

         Do j = 1, 4 ! loop over faces
            k = Face2OpPnt(j)

            s = face_area(k) / (3*vol) * localF(j)

            i1 = tet(k, n)
            Do i = 1, 3
               VECatTets(i,n) = VECatTets(i,n) + s * (xyzc(i)-vrt(i,i1))
            End do
         End do
      End do

      Return
      End

max_stars_repo_name: INM-RAS/INMOST

max_issues_repo_name: INMOST-DEV/INMOST

max_forks_repo_name: INM-RAS/INMOST


text:
      subroutine elseif (lab)
      integer lab
      call outgo (lab+1)
      call indent (-1)
      call outcon (lab)
      call indent (1)
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


program test_io_params

  use ioparams, only: group1_t, group2_t
  use ioparams, only: read_nml, write_nml, read_nml_file
  use ioparams, only: set_param, get_param
  use ioparams, only: parse_command_args, print_help, join_array
  use ioparams, only: command_argument_as_array
  use ioparams, only: has_param, set_param_string ! low-level

  implicit none

  type(group1_t) :: group1
  type(group2_t) :: group2
  character(len=256) :: filename, filename2
  integer :: iounit = 88
  double precision :: test_dp
  logical :: test_l
  integer :: test_int, test_int_arr(7)
  character(len=50) :: test_s
  character(len=50) :: args(14)
  character(len=50), allocatable :: unmatched(:)

  integer :: i, parsed, io, iostat, iostats(2), stat
  character(len=256) :: arg, argv

  filename = "namelist.nml"

  write(*,*) " "
  write(*,*) "++++++++++++++++++++++++++++++++++++"
  write(*,*) "Testing nml2f90.py"
  write(*,*) "++++++++++++++++++++++++++++++++++++"
  write(*,*) " "
  write(*,*) "Test reading namelist"
  write(*,*) "---------------------"
  ! Read parameters from file
  write(*,*) "READ from ",trim(filename)
  open(iounit, file=filename, status="OLD")
  write(*,*) "...group1...(read_nml)"
  call read_nml(iounit, group1)
  close(iounit)
  write(*,*) "...group2...(read_nml_file)"
  call read_nml_file(filename, group2)
  write(*,*) "Done. Namelist read successfully."
  write(*,*) " "
  write(*,*) "Test set_param / get_param "
  write(*,*) "---------------------------"
  call get_param(group2, "string1", test_s)
  call set_param(group2, 'string1', "this is a new string set via set_param")
  write(*,*) "group2%string1:", trim(test_s), " ==> ", trim(group2%string1), " (this is ...)"

  call get_param(group2, "double1", test_dp)
  call set_param(group2, 'double1', 11111111111.11111111111d0)
  write(*,*) "group2%double1:", test_dp, " ==> ", group2%double1, " (111...)"

  call get_param(group2, 'integer1', test_int)
  call set_param(group2, 'integer1', 777777777)
  write(*,*) "group2%integer1:", test_int, " ==> ", group2%integer1, " (777..)"

  call get_param(group2, "intarr1", test_int_arr)
  call set_param(group2, "intarr1", [7,6,5,4,3,2,1])
  write(*,*) "group2%intarr1:", test_int_arr, " ==> ", group2%intarr1, " (7,6,5...)"

  call get_param(group2, "logical1", test_l)
  call set_param(group2, "logical1", .false.)
  write(*,*) "group2%longical1:", test_l, " ==> ", group2%logical1, " (false)"

  write(*,*) " "
  write(*,*) "Test has_param"
  write(*,*) "--------------"

  write(*,*) "has_param(group1, 'integer1') ? ", has_param(group1, 'integer1')
  write(*,*) "has_param(group1, 'thisdoesnotexists') ? ", has_param(group1, 'thisdoesnotexists')
  write(*,*)
  write(*,*) "Test set_param_string"
  write(*,*) "---------------------"
  call set_param_string(group2, 'string1', "another string")
  call set_param_string(group2, 'double1', "3333.44444")
  call set_param_string(group2, 'integer1', "3333")
  call set_param_string(group2, "intarr1", "(/ 7,7,7,7,7,7,7 /)")
  call set_param_string(group2, "logical1", ".true.")

  write(*,*)
  write(*,*) "Test command-line arguments passed as string array"
  write(*,*) "--------------------------------------------------"
  args(1) = "--string1"
  args(2) = "newtext.txt"
  args(3) = "--logical1"
  args(4) = "F"
  args(5) = "--integer1"
  args(6) = "2"
  args(7) = "--double1"
  args(8) = "777"
  args(9) = "--stringarr1"
  args(10) = "bb,cc,dd"
  args(11) = "--intarr1"
  args(12) = "1, 2, 3, 4, 5, 6, 7"
  args(13) = "--dblarr1"
  args(14) = "5.,4,3,2,1"
  call parse_command_args(group1, stop_on_help=.false., args=args, unmatched=unmatched)
  call parse_command_args(group2, args=unmatched) !, unmatched=unmatched)
  call assert_cmd("group1%string1", "newtext.txt"== group1%string1)
  call assert_cmd("group1%logical1", .false. .eqv. group1%logical1)
  call assert_cmd("group1%integer1", 2== group1%integer1)
  call assert_cmd("group2%double1", 777.d0== group2%double1)
  call assert_cmd("group1%stringarr1(2)", "cc"== group1%stringarr1(2))
  call assert_cmd("group2%intarr1(3)", 3== group2%intarr1(3))
  call assert_cmd("group2%dblarr1(1)", 5.d0== group2%dblarr1(1))

  !
  ! Retrieve command-line arguments
  !
  write(*,*) " "
  write(*,*) "Test command line parameters  "
  write(*,*) "----------------------------- "
  write(*,*) "Type ./test.x -h for help on how to do that."
  i = 1
  io=0
  parsed = 0
  call command_argument_as_array(unmatched)
  call parse_command_args(group1, args=unmatched, unmatched=unmatched, stop_on_help=.false.)
  call parse_command_args(group2, args=unmatched, unmatched=unmatched)
  if (size(unmatched) > 0) then
    write(*,*) "Invalid parameters: ", trim(join_array(unmatched))
    stop
  endif

  ! Print namelist to screen
  write(*,*) " "
  write(*,*) "Test writing to namelist_check.nml "

  open(88, file="namelist_check.nml", action="WRITE")
  call write_nml(88, group1)
  call write_nml(88, group2)
  close(88)
  write(*,*) "Done. All tests terminated successfully"

contains

  subroutine assert_cmd(nm, test)
    character(len=*), intent(in) :: nm
    logical, intent(in) :: test
    if (.not.test) then
      write(*,*) "ERROR: failed cmd",nm
      stop
    endif
  end subroutine
end program


text:
module prop

  complex*16, save, public                 ::  de, e3, e0
  complex*16, save, public, allocatable    ::  mut0(:,:), muts(:,:), mutl(:,:)

  public :: calc_mut0
  public :: calc_muts
  public :: calc_mutl
  public :: calc_de
  public :: calc_3de
  public :: calc_e0

  contains

  subroutine calc_mut0(dnst)

    use param,         only : nsite, norb, sigma_x, sigma_y, sigma_z
    use constants,     only : nspin, cmplx_0, cmplx_i
    use utilities,     only : trace

    implicit none

    integer                :: isite
    complex*16, intent(in) :: dnst(norb,norb,nsite,nspin,nspin)

    allocate(mut0(3,nsite))

    mut0 = cmplx_0
    do isite = 1, nsite
       mut0(1,isite) = mut0(1,isite) +         trace(dnst(:,:,isite,1,2)+dnst(:,:,isite,2,1),norb)
       mut0(2,isite) = mut0(2,isite) + cmplx_i*trace(dnst(:,:,isite,2,1)-dnst(:,:,isite,1,2),norb)
       mut0(3,isite) = mut0(3,isite) +         trace(dnst(:,:,isite,1,1)-dnst(:,:,isite,2,2),norb)
    end do

  end subroutine calc_mut0

  subroutine calc_muts(ddens)

    use param,          only : nsite, norb, sigma_x, sigma_y, sigma_z
    use constants,      only : nspin, cmplx_0, cmplx_i
    use utilities,      only : trace

    implicit none

    integer                :: isite
    complex*16, intent(in) :: ddens(norb,norb,nsite,nspin,nspin)
    complex*16             :: dmuts(3,nsite)

    allocate(muts(3,nsite))

    muts = cmplx_0
    dmuts = cmplx_0
    do isite = 1, nsite
       dmuts(1,isite) = dmuts(1,isite) +         trace(ddens(:,:,isite,1,2)+ddens(:,:,isite,2,1),norb)
       dmuts(2,isite) = dmuts(2,isite) + cmplx_i*trace(ddens(:,:,isite,2,1)-ddens(:,:,isite,1,2),norb)
       dmuts(3,isite) = dmuts(3,isite) +         trace(ddens(:,:,isite,1,1)-ddens(:,:,isite,2,2),norb)
    end do

    muts(:,:) = mut0(:,:) + dmuts(:,:)

  end subroutine calc_muts

  subroutine calc_mutl(ddens)

    use param,          only : nsite, norb, lx, ly, lz
    use constants,      only : nspin, cmplx_0
    use utilities,      only : trace

    implicit none

    integer                :: isite, iorb, iorb1, ispin, ispin1
    complex*16, intent(in) :: ddens(norb,norb,nsite,nspin,nspin)

    allocate(mutl(3,nsite))

    mutl = cmplx_0
    do isite = 1, nsite
       do iorb = 1, norb
          do iorb1 = 1, norb
                mutl(1,isite) = mutl(1,isite) + lx(iorb,iorb1)*(ddens(iorb,iorb1,isite,1,1) + ddens(iorb,iorb1,isite,2,2))
                mutl(2,isite) = mutl(2,isite) + ly(iorb,iorb1)*(ddens(iorb,iorb1,isite,1,1) + ddens(iorb,iorb1,isite,2,2))
                mutl(3,isite) = mutl(3,isite) + lz(iorb,iorb1)*(ddens(iorb,iorb1,isite,1,1) + ddens(iorb,iorb1,isite,2,2))
          end do
       end do
    end do

  end subroutine calc_mutl

  subroutine calc_e0(dens)

    use param,          only : nsite, norb, hf_val, hf_pot, nkpt, efermi
    use constants,      only : nspin, cmplx_0
    use utilities,      only : trace

    implicit none

    integer                :: itot, ispin, ispin1, ikpt, iorb, iorb1, isite
    complex*16, intent(in) :: dens(norb,norb,nsite,nspin,nspin)

    e0 = 0.0
    do itot = 1, norb*nsite
       do ispin = 1, nspin
          do ikpt = 1, nkpt
             if (hf_val(itot,ispin,ikpt) > efermi) then
                cycle
             end if

             e0 = e0 + hf_val(itot,ispin,ikpt)/real(nkpt)
          end do
       end do
    end do

    do isite = 1, nsite
       do iorb = 1, norb
          do iorb1 = 1, norb
             e0 = e0 - 0.5*hf_pot(iorb,iorb1,isite,1)*dens(iorb,iorb1,isite,1,1)
             e0 = e0 - 0.5*hf_pot(iorb,iorb1,isite,2)*dens(iorb,iorb1,isite,2,2)
          end do
       end do
    end do

  end subroutine calc_e0

  subroutine calc_de(dvpot,ddens)

    use param,          only : nsite, norb
    use constants,      only : nspin, cmplx_0
    use utilities,      only : trace

    implicit none

    integer                :: isite, ispin, ispin1, iorb, iorb1, icomb, icomb1

    complex*16, intent(in) :: dvpot(norb,norb,nsite,nspin,nspin),&
                              ddens(norb,norb,nsite,nspin,nspin)
    de = 0.0
    do ispin = 1, nspin
       do ispin1 = 1, nspin
          do isite = 1, nsite
             do iorb = 1, norb
                do iorb1 = 1, norb
                   de = de + 0.5*dvpot(iorb,iorb1,isite,ispin,ispin1)*ddens(iorb,iorb1,isite,ispin,ispin1)
                end do
             end do
          end do
       end do
    end do

  end subroutine calc_de

  subroutine calc_3de(dvpot,pot)

    use param,                only : nrpt, ntot, norb, nsite, r_bond, nkpt, kpt_cart, &
                                     fdtemp, efermi, nelec
    use constants,            only : nspin, hamname, cmplx_i, kboltz, cmplx_0
    use utilities,            only : cutility_diag
    use fractions,            only : fermi_dirac

    implicit none

    integer                :: i, ikpt, irpt, itot, iat, iat1, dum, ispin, ispin1, isite, isite1, iorb, iorb1
    logical                :: io_find
    real*8                 :: eigval3(ntot,nkpt), emax, emin, numax, numin, num, e
    complex*16             :: ham(norb,norb,nspin,nspin), ham_k(norb,norb,nsite,nsite,nkpt,nspin,nspin), &
                              vnew(norb,norb,nsite,nspin,nspin), ham_full(ntot,ntot,nkpt), temp(ntot,ntot), &
                              eigvec3(ntot,ntot,nkpt), newdens(norb,norb,nsite,nspin,nspin)
    complex*16             :: rdotk, fac
    complex*16, intent(in) :: dvpot(norb,norb,nsite,nspin,nspin), pot(norb,norb,nsite,nspin,nspin)

    inquire(file=hamname, exist=io_find)
    if (io_find) then
       ham_k = cmplx_0
       do ikpt = 1, nkpt
          open(2, file=hamname, status='old')
          ham = cmplx_0
          do irpt = 1, nrpt
             read(2,*) iat, iat1, dum, dum
             read(2,*)
             do iorb = 1, norb
                read(2,*) ham(iorb,:,1,1)
             end do

             !Expand to spindown subspace
             ham(:,:,2,2) = ham(:,:,1,1)

             rdotk = dot_product(kpt_cart(:,ikpt),r_bond(:,irpt))
             fac = exp(cmplx_i*rdotk)
             ham_k(:,:,iat,iat1,ikpt,1,1) = ham_k(:,:,iat,iat1,ikpt,1,1) + fac*ham(:,:,1,1)
             ham_k(:,:,iat,iat1,ikpt,2,2) = ham_k(:,:,iat,iat1,ikpt,2,2) + fac*ham(:,:,2,2)
             ham_k(:,:,iat,iat1,ikpt,1,2) = ham_k(:,:,iat,iat1,ikpt,1,2) + fac*ham(:,:,1,2)
             ham_k(:,:,iat,iat1,ikpt,2,1) = ham_k(:,:,iat,iat1,ikpt,2,1) + fac*ham(:,:,2,1)
          end do
          close(2)
       end do
    else
       print*, 'No sclr_ham.in file found'
       stop
    end if

    ! New hf potential as V_hf + dvp

    vnew = pot + dvpot

    !Sum kinetic part to the hamiltonian
    ham_full = cmplx_0
    do ikpt = 1, nkpt
       do ispin = 1, nspin
          do ispin1 = 1, nspin
             do isite = 1, nsite
                do isite1 = 1, nsite
                   ham_full((isite-1)*norb+1+(ispin-1)*norb*nsite:isite*norb+(ispin-1)*norb*nsite, &
                            (isite1-1)*norb+1+(ispin1-1)*norb*nsite:isite1*norb+(ispin1-1)*norb*nsite,ikpt) = &
                            ham_k(:,:,isite,isite1,ikpt,ispin,ispin1)
                end do
             end do
          end do
       end do
    end do

    !Sum hf_potential
    do ikpt = 1, nkpt
       do ispin = 1, nspin
          do ispin1 = 1, nspin
             do isite = 1, nsite
                ham_full((isite-1)*norb+1+(ispin-1)*norb*nsite:isite*norb+(ispin-1)*norb*nsite, &
                         (isite-1)*norb+1+(ispin1-1)*norb*nsite:isite*norb+(ispin1-1)*norb*nsite,ikpt) =   &
                ham_full((isite-1)*norb+1+(ispin-1)*norb*nsite:isite*norb+(ispin-1)*norb*nsite, &
                         (isite-1)*norb+1+(ispin1-1)*norb*nsite:isite*norb+(ispin1-1)*norb*nsite,ikpt) +   &
                vnew(:,:,isite,ispin,ispin1)
             end do
          end do
       end do
    end do

    ! Diagonalize
    do ikpt = 1, nkpt
       temp(:,:) = ham_full(:,:,ikpt)
       call cutility_diag(temp(:,:),eigval3(:,ikpt),ntot)
       eigvec3(:,:,ikpt) = temp(:,:)
    end do

    ! Calculate new efermi
    !emax = maxval(eigval3(ntot,:))
    !emin = minval(eigval3(1,:))
    !numax = 0.0
    !numin = 0.0
    !do itot = 1, ntot
    !   do ikpt = 1, nkpt
    !      numax = numax + fermi_dirac((eigval3(itot,ikpt)-emax))/real(nkpt)
    !      numin = numin + fermi_dirac((eigval3(itot,ikpt)-emin))/real(nkpt)
    !   end do
    !end do
!
!    if ((numax > nelec) .and. (numin < nelec)) then
!       do i = 1, 30
!          num = 0
!          e = (emax-emin)/2.0 + emin
!          do itot = 1, ntot
!             do ikpt = 1, nkpt
!                num = num + fermi_dirac((eigval3(itot,ikpt)-e))/real(nkpt)
!             end do
!          end do
!          if (abs(num-nelec) < 1e-6) then
!             efermi = e
!             exit
!          else if (num < nelec) then
!             emin = e
!          else if (num > nelec) then
!             emax = e
!          end if
!       end do
!    end if

    ! Calculate new density matrix
    newdens = cmplx_0
    do itot = 1, ntot
       do ikpt = 1, nkpt
          !if (eigenval(itot,ikpt) > efermi) then
          !   cycle
          !end if
          do isite = 1, nsite
             do ispin = 1, nspin
                do ispin1 = 1, nspin
                   do iorb = 1, norb
                      do iorb1 = 1, norb
                         newdens(iorb,iorb1,isite,ispin,ispin1) = &
                         newdens(iorb,iorb1,isite,ispin,ispin1) + fermi_dirac(eigval3(itot,ikpt)-efermi)*  &
                                                                   conjg(eigvec3((isite-1)*norb+iorb+(ispin-1)*norb*nsite,itot,ikpt))* &
                                                                   eigvec3((isite-1)*norb+iorb1+(ispin1-1)*norb*nsite,itot,ikpt)/nkpt
                      end do
                   end do
                end do
             end do
          end do
       end do
    end do

    ! Calculate energy
    e3 = 0.0
    do itot = 1, ntot
       do ikpt = 1, nkpt
          if (eigval3(itot,ikpt) > efermi) then
             cycle
          end if

          e3 = e3 + eigval3(itot,ikpt)/nkpt

       end do
    end do

    do ispin = 1, nspin
       do ispin1 = 1, nspin
          do isite = 1, nsite
             do iorb = 1, norb
                do iorb1 = 1, norb
                   e3 = e3 - 0.5*vnew(iorb,iorb1,isite,ispin,ispin1)*newdens(iorb,iorb1,isite,ispin,ispin1)
                end do
             end do
          end do
       end do
    end do

  end subroutine calc_3de

end module prop

max_stars_repo_name: oarcelus/HFPP

max_issues_repo_name: oarcelus/HFPP

max_forks_repo_name: oarcelus/HFPP


text:
module init_fluid_module

      use amrex_fort_module, only: rt => amrex_real
      use iso_c_binding,     only: c_int

contains
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: init_fluid                                              !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
   subroutine init_fluid(slo, shi, lo, hi, &
                         domlo, domhi, p, vel, density, tracer, eta, &
                         dx, dy, dz, xlo, ylo, zlo, probtype) &
      bind(C, name="init_fluid")

      use constant, only: ro_0, mu
      use constant, only: ic_u, ic_v, ic_w, ntrac

      implicit none

! Dummy arguments .....................................................//
      integer(c_int), intent(in   ) ::  lo(3),  hi(3)
      integer(c_int), intent(in   ) :: slo(3), shi(3)
      integer(c_int), intent(in   ) :: domlo(3),domhi(3)

      real(rt), intent(inout) ::       p(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3)  )
      real(rt), intent(inout) ::     vel(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3),3)
      real(rt), intent(inout) :: density(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3))
      real(rt), intent(inout) ::  tracer(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3),ntrac)
      real(rt), intent(inout) ::     eta(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3)  )

      real(rt), intent(in   ) :: dx, dy, dz
      real(rt), intent(in   ) :: xlo,ylo,zlo
      integer,  intent(in   ) :: probtype

      ! Set the initial fluid density and viscosity
      vel(:,:,:,1) = ic_u
      vel(:,:,:,2) = ic_v
      vel(:,:,:,3) = ic_w

      density = ro_0
      tracer  = 0.0
      eta     = mu

      if (probtype ==  5) call tracer_advection(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo)

      if (probtype ==  1) call taylor_green(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
      if (probtype ==  2) call double_shear_layer(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
      if (probtype ==  3) call taylor_green3d(lo, hi, vel, slo, shi, dx, dy, dz, domlo)
      if (probtype == 31) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
      if (probtype == 32) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
      if (probtype == 33) call plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)
      if (probtype ==  4) call couette(lo, hi, vel, slo, shi, dx, dy, dz, domlo, domhi)
      if (probtype == 11) call tuscan(lo, hi, vel, density, tracer, slo, shi, dx, dy, dz, domlo, domhi)
      if (probtype == 34) call channel_flow(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, xlo, ylo, zlo)

   end subroutine init_fluid

   subroutine tracer_advection(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo)

      use constant, only: zero, half, one

      implicit none

      ! Array bounds
      integer(c_int), intent(in   ) :: lo(3), hi(3)
      integer(c_int), intent(in   ) :: domlo(3)
      integer(c_int), intent(in   ) :: slo(3), shi(3)

      ! Arrays
      real(rt),       intent(inout) ::    vel(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)
      real(rt),       intent(inout) :: tracer(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)

      real(rt),       intent(in   ) :: dx, dy, dz

      ! Local variables
      integer(c_int)                  :: i, j, k
      real(rt)                        :: x, y
      real(rt)                        :: twopi = 8.0_rt * atan(one)

      ! x-direction
      do j = lo(2), hi(2)
         y =  ( real(j,rt) + half ) * dy
         do i = lo(1), hi(1)
            x =  ( real(i,rt) + half ) * dx
            do k = lo(3), hi(3)
               tracer(i,j,k,:) = 1.d0
               if (i.eq.20 .and. j.eq.20 .and. k.eq.20) tracer(i,j,k,1) = 1.2d0
            end do
         end do
      end do

   end subroutine tracer_advection

   subroutine taylor_green(lo, hi, vel, slo, shi, dx, dy, dz, domlo)

      use constant, only: zero, half, one

      implicit none

      ! Array bounds
      integer(c_int), intent(in   ) :: lo(3), hi(3)
      integer(c_int), intent(in   ) :: domlo(3)
      integer(c_int), intent(in   ) :: slo(3), shi(3)

      ! Arrays
      real(rt),       intent(inout) :: vel(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)

      real(rt),       intent(in   ) :: dx, dy, dz

      ! Local variables
      integer(c_int)                  :: i, j, k
      real(rt)                        :: x, y
      real(rt)                        :: twopi = 8.0_rt * atan(one)

      ! x-direction
      do j = lo(2), hi(2)
         y =  ( real(j,rt) + half ) * dy
         do i = lo(1), hi(1)
            x =  ( real(i,rt) + half ) * dx
            do k = lo(3), hi(3)
               vel(i,j,k,1) =  sin(twopi * x) * cos(twopi * y)
               vel(i,j,k,2) = -cos(twopi * x) * sin(twopi * y)
               vel(i,j,k,3) = zero
            end do
         end do
      end do

   end subroutine taylor_green

    subroutine taylor_green3d(lo, hi, vel, slo, shi, dx, dy, dz, domlo)

       use constant, only: zero, half, one

       implicit none

       ! Array bounds
       integer(c_int), intent(in   ) :: lo(3), hi(3)
       integer(c_int), intent(in   ) :: domlo(3)
       integer(c_int), intent(in   ) :: slo(3), shi(3)

       ! Arrays
       real(rt),       intent(inout) :: vel(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)

       real(rt),       intent(in   ) :: dx, dy, dz

       ! Local variables
       integer(c_int)                  :: i, j, k
       real(rt)                        :: x, y, z
       real(rt)                        :: twopi = 8.0_rt * atan(one)

       do j = lo(2), hi(2)
          y =  ( real(j,rt) + half ) * dy
          do i = lo(1), hi(1)
             x =  ( real(i,rt) + half ) * dx
             do k = lo(3), hi(3)
                z =  ( real(k,rt) + half ) * dz
                vel(i,j,k,1) =  sin(twopi * x) * cos(twopi * y) * cos(twopi * z)
                vel(i,j,k,2) = -cos(twopi * x) * sin(twopi * y) * cos(twopi * z)
                vel(i,j,k,3) = zero
             end do
          end do
       end do

    end subroutine taylor_green3d

   subroutine double_shear_layer(lo, hi, vel, slo, shi, dx, dy, dz, domlo)

      use constant, only: zero, half, one

      implicit none

      ! Array bounds
      integer(c_int), intent(in   ) :: lo(3), hi(3)
      integer(c_int), intent(in   ) :: domlo(3)
      integer(c_int), intent(in   ) :: slo(3), shi(3)

      ! Arrays
      real(rt),       intent(inout) :: vel(slo(1):shi(1),slo(2):shi(2),slo(3):shi(3),3)

      real(rt),       intent(in   ) :: dx, dy, dz

      ! Local variables
      integer(c_int)                :: i, j, k, plane
      real(rt)                      :: x, y, z
      real(rt)                      :: twopi = 8.0_rt * atan(one)

      plane = 1

      select case ( plane )

      case (1)  ! x-y plane

         ! x-direction
         do j = lo(2), hi(2)
            y =  ( real(j,rt) + half ) * dy
            do i = lo(1), hi(1)
               x =  ( real(i,rt) + half ) * dx
               do k = lo(3), hi(3)
                  vel(i,j,k,1) = tanh ( 30.0_rt * (0.25_rt - abs ( y - 0.5_rt ) ) )
                  vel(i,j,k,2) = 0.05_rt * sin ( twopi * x )
                  vel(i,j,k,3) = zero
               end do
            end do
         end do

      case (2)  ! x-z plane

         ! x-direction
         do k = lo(3), hi(3)
            z =  ( real(k,rt) + half ) * dz
            do i = lo(1), hi(1)
               x =  ( real(i,rt) + half ) * dx
               do j = lo(2), hi(2)
                  vel(i,j,k,1) = tanh ( 30.0_rt * (0.25_rt - abs ( z - 0.5_rt ) ) )
                  vel(i,j,k,2) = zero
                  vel(i,j,k,3) = 0.05_rt * sin ( twopi * x )
               end do
            end do
         end do

      case (3)  ! y-z plane

         ! x-direction
         do k = lo(3), hi(3)
            z =  ( real(k,rt) + half ) * dz
            do j = lo(2), hi(2)
               y =  ( real(j,rt) + half ) * dy
               do i = lo(1), hi(1)
                  vel(i,j,k,1) = zero
                  vel(i,j,k,2) = tanh ( 30.0_rt * (0.25_rt - abs ( z - 0.5_rt ) ) )
                  vel(i,j,k,3) = 0.05_rt * sin ( twopi * y )
               end do
            end do
         end do

      end select

   end subroutine double_shear_layer

   subroutine plane_poiseuille(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, probtype)

      use constant, only: ic_u, ic_v, ic_w, ntrac

      implicit none

      integer(c_int),   intent(in   ) ::    lo(3),    hi(3)
      integer(c_int),   intent(in   ) :: domlo(3), domhi(3)
      integer(c_int),   intent(in   ) ::   slo(3),   shi(3)
      integer(c_int),   intent(in   ) :: probtype

      real(rt),         intent(in   ) :: dx, dy, dz
      real(rt),         intent(inout) ::    vel(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3), 3)
      real(rt),         intent(inout) :: tracer(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3), ntrac)

      ! Local variables
      integer(c_int)                  :: i, j, k
      integer(c_int)                  :: num_cells
      real(rt)                        :: x, y, z

      if (probtype .eq. 31) then

         num_cells = domhi(2) - domlo(2) + 1

!        if (ntrac .ne. 3) then
!           print *,"ntrac here is ", ntrac
!           print *,"probtype = 31 is supposed to have 3 tracers!"
!           stop
!        end if

         do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            do i = lo(1), hi(1)
               y =  (real(j,rt) + 0.5d0) / dble(num_cells)
               x =  (real(i,rt) + 0.5d0) / dble(num_cells)
               vel(i,j,k,1) = 6.0d0 * ic_u * y * (1.0d0 - y)
               vel(i,j,k,2) = 0.0d0
               vel(i,j,k,3) = 0.0d0
               tracer(i,j,k,1:ntrac) = 0.0d0
               if (ntrac .gt. 0 .and. i .le. domhi(1)/8) &
                  tracer(i,j,k,1) = 1.d0
               if (ntrac .gt. 1 .and. i .le. domhi(1)/2) &
                  tracer(i,j,k,2) = 2.d0
               if (ntrac .gt. 2 .and. i .le. 3*domhi(1)/4) &
                  tracer(i,j,k,3) = 3.d0
            end do
         end do
         end do

      else if (probtype .eq. 32) then

         num_cells = domhi(3) - domlo(3) + 1

         do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            do i = lo(1), hi(1)
               z =  (real(k,rt) + 0.5d0) / dble(num_cells)
               vel(i,j,k,1) = 0.0d0
               vel(i,j,k,2) = 6.0d0 * ic_v * z * (1.0d0 - z)
               vel(i,j,k,3) = 0.0d0
               tracer(i,j,k,1:ntrac) = 0.0d0
               if (ntrac .gt. 0 .and. j .le. domhi(2)/8) &
                  tracer(i,j,k,1) = 1.d0
               if (ntrac .gt. 1 .and. j .le. domhi(2)/2) &
                  tracer(i,j,k,2) = 2.d0
               if (ntrac .gt. 2 .and. j .le. 3*domhi(2)/4) &
                  tracer(i,j,k,3) = 3.d0
            end do
         end do
         end do

      else if (probtype .eq. 33) then

         num_cells = domhi(1) - domlo(1) + 1

         do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            do i = lo(1), hi(1)
               x =  (real(i,rt) + 0.5d0) / dble(num_cells)
               vel(i,j,k,1) = 0.0d0
               vel(i,j,k,2) = 0.0d0
               vel(i,j,k,3) = 6.0d0 * ic_w * x * (1.0d0 - x)
               tracer(i,j,k,1:ntrac) = 0.0d0
               if (ntrac .gt. 0 .and. k .le. domhi(3)/8) &
                  tracer(i,j,k,1) = 1.d0
               if (ntrac .gt. 1 .and. k .le. domhi(3)/2) &
                  tracer(i,j,k,2) = 2.d0
               if (ntrac .gt. 2 .and. k .le. 3*domhi(3)/4) &
                  tracer(i,j,k,3) = 3.d0
            end do
         end do
         end do

      else
          print *,'Dont know this probtype in plane_poiseuille: ', probtype
          stop
      end if

   end subroutine plane_poiseuille

   subroutine couette(lo, hi, vel, slo, shi, dx, dy, dz, domlo, domhi)

      use constant,          only: zero, half, one
      use constant,          only: ic_u

      implicit none

      integer(c_int),   intent(in   ) ::    lo(3),    hi(3)
      integer(c_int),   intent(in   ) :: domlo(3), domhi(3)
      integer(c_int),   intent(in   ) ::   slo(3),   shi(3)

      real(rt),         intent(inout) :: vel(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3), 3)

      real(rt),         intent(in   ) :: dx, dy, dz

      ! Local variables
      integer(c_int)                  :: i, j, k
      integer(c_int)                  :: num_cells_y
      real(rt)                        :: y

      num_cells_y = domhi(2) - domlo(2) + 1

      do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            y =  (real(j,rt) + half) / num_cells_y
            do i = lo(1), hi(1)
               vel(i,j,k,1) = ic_u * (y - half)
               vel(i,j,k,2) = zero
               vel(i,j,k,3) = zero
            end do
         end do
      end do

   end subroutine couette

   subroutine tuscan(lo, hi, vel, density, tracer, slo, shi, dx, dy, dz, domlo, domhi)

      use constant,          only: zero, half, one

      implicit none

      integer(c_int),   intent(in   ) ::    lo(3),    hi(3)
      integer(c_int),   intent(in   ) :: domlo(3), domhi(3)
      integer(c_int),   intent(in   ) ::   slo(3),   shi(3)

      real(rt),         intent(inout) ::     vel(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3), 3)
      real(rt),         intent(inout) :: density(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3))
      real(rt),         intent(inout) ::  tracer(slo(1):shi(1), slo(2):shi(2), slo(3):shi(3))

      real(rt),         intent(in   ) :: dx, dy, dz

      ! Local variables
      integer(c_int)                  :: i, j, k
      integer(c_int)                  :: num_cells
      real(rt)                        :: pert

      num_cells = domhi(3) - domlo(3) + 1

      pert = 0.01d0

      ! Start the flow at rest
      vel = zero

      do k = lo(3), hi(3)
         do j = lo(2), hi(2)
            do i = lo(1), hi(1)
               if (k .le. num_cells/2) then
                  density(i,j,k) = one
                  tracer (i,j,k) = zero
               else
                  density(i,j,k) = one
                  tracer (i,j,k) = zero + pert
               endif
            end do
         end do
      end do

   end subroutine tuscan

   subroutine channel_flow(lo, hi, vel, tracer, slo, shi, dx, dy, dz, domlo, domhi, xlo, ylo, zlo)

      use constant, only: ic_u, ic_v, ic_w, ntrac,one

      implicit none

      integer(c_int),   intent(in   ) ::    lo(3),    hi(3)
      integer(c_int),   intent(in   ) :: domlo(3), domhi(3)
      integer(c_int),   intent(in   ) ::   slo(3),   shi(3)

      real(rt),         intent(in   ) :: dx, dy, dz, xlo, ylo, zlo
      real(rt),         intent(inout) ::    vel(slo(1):shi(1), slo(2):shi(2),slo(3):shi(3), 3)
      real(rt),         intent(inout) :: tracer(slo(1):shi(1), slo(2):shi(2),slo(3):shi(3), ntrac)

      ! Local variables
      integer(c_int)                  :: i, j, k
      integer(c_int)                  :: num_cells, num_cells_x, num_cells_y,num_cells_z
      real(rt)                        :: x, y, z, yval, zval, pi

      pi = 4.0_rt * atan(one)

      do k = lo(3), hi(3)
         z = zlo + dz*(real(k,rt) + 0.5d0)
         do j = lo(2), hi(2)
            y = ylo + dy*(real(j,rt) + 0.5d0)
            do i = lo(1), hi(1)
               x = xlo + dx*(real(i,rt) + 0.5d0)

               vel(i,j,k,1) = 8.0d0*(5.0d0*(1.0d0-z**4)/4.0d0 + 0.3d0*cos(5.0d0*pi*y)*exp(0.5-32.4*(1.0d0-abs(z))**2)*(1.0d0-abs(z)))
               vel(i,j,k,2) = 8.0d0*21.6d0*sin(4.0d0*pi*x)*exp(-32.4d0*(1.0d0-abs(z))**2)*(1.0d0-abs(z))
               vel(i,j,k,3) = 0.0d0
               tracer(i,j,k,1:ntrac) = 0.0d0
               if (ntrac .gt. 0 .and. i .le. domhi(1)/4) &
                  tracer(i,j,k,1) = 1.d0
               if (ntrac .gt. 1 .and. i .le. domhi(1)/2) &
                  tracer(i,j,k,2) = 2.d0
               if (ntrac .gt. 2 .and. i .le. 3*domhi(1)/4) &
                  tracer(i,j,k,3) = 3.d0
            end do
         end do
      end do

   end subroutine channel_flow

end module init_fluid_module

max_stars_repo_name: psakievich/amr-wind

max_issues_repo_name: psakievich/amr-wind

max_forks_repo_name: psakievich/amr-wind


text:
      subroutine runY_000000l(k,l,Xtwiddle,Gtwiddle,Shat7zzzz,N0)
      implicit none
C---  Expression for extension of Eq. 5.60a
C---  Calculates D00llll
C---  Small terms of order Xtwiddle(0,k)*Diiiii,Xtwiddle(0,0)*Diiiiii
C---  Denominator Gtwiddle(k,l)
      include 'pvDnames.f'
      include 'pvDv.f'
      include 'Darraydef.f'
      include 'Darrays.f'
      integer ep,N0,k,l,np
      parameter(np=3)
      double precision Xtwiddle(0:np,0:np),Gtwiddle(np,np)
      double complex Shat7zzzz(np,z2max,-2:0)

      do ep=-2,0
      Dv(dzzzzzzi(l)+N0,ep)=
     . (Gtwiddle(k,1)*Shat7zzzz(1,z2(l,l),ep)
     . +Gtwiddle(k,2)*Shat7zzzz(2,z2(l,l),ep)
     . +Gtwiddle(k,3)*Shat7zzzz(3,z2(l,l),ep)
     . +Xtwiddle(k,0)*Dv(dzzzzii(z2(l,l))+N0,ep)
     . -Xtwiddle(0,0)*Dv(dzzzziii(z3(k,l,l))+N0,ep))
     . /(4d0*Gtwiddle(k,l))
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      program t
      print *,"don't"
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
c     UC Copyright Notice
c
c     This software is Copyright (c) 2014-2015 The Regents of the
c     University of California. All Rights Reserved.
c
c     Permission to copy and modify this software and its documentation
c     for educational, research and non-profit purposes, without fee,
c     and without a written agreement is hereby granted, provided that
c     the above copyright notice, this paragraph and the following three
c     paragraphs appear in all copies.
c
c     Permission to make commercial use of this software may be obtained
c     by contacting:
c
c     Technology Transfer Office
c     9500 Gilman Drive, Mail Code 0910
c     University of California
c     La Jolla, CA 92093-0910
c     (858) 534-5815
c     invent@ucsd.edu
c
c     This software program and documentation are copyrighted by The
c     Regents of the University of California. The software program and
c     documentation are supplied "as is", without any accompanying
c     services from The Regents. The Regents does not warrant that the
c     operation of the program will be uninterrupted or error-free. The
c     end-user understands that the program was developed for research
c     purposes and is advised not to rely exclusively on the program for
c     any reason.
c
c     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
c     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
c     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
c     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
c     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
c     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY
c     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
c     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
c     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE
c     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
c     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

#include "cvFlowsolverOptions.h"

#if (VER_CLOSEDLOOP == 1)

!> This subroutine initiate the General boundary condition

      SUBROUTINE initGenBC

      USE GeneralBC
      INCLUDE "global.h"
      INCLUDE "mpif.h"
      INCLUDE "common_blocks/workfc.h"
      INCLUDE "common_blocks/conpar.h"
      INCLUDE "common_blocks/nomodule.h"

      REAL*8 temp(nshg)
      INTEGER loopA
      LOGICAL ierr

      ALLOCATE(nsrflistCoupled(0:MAXSURF), ECoupled(0:MAXSURF),
     2   QCoupled(0:MAXSURF), QnCoupled(0:MAXSURF), ACoupled(0:MAXSURF),
     2   PCoupled(0:MAXSURF), PnCoupled(0:MAXSURF), PGenDer(0:MAXSURF))

      IF (ipvsq .GE. 2) THEN
         PGenDerFlag = .TRUE.
      ELSE
         PGenDerFlag = .FALSE.
      END IF
      PGenDer = 0D0

      nsrflistCoupled = nsrflistDirichlet
      nsrflistCoupled(numDirichletSrfs+1:numCoupledSrfs) =
     2   nsrflistNeumann(1:numNeumannSrfs)

c     Eval Area of Coupled Surfaces
      temp = 1D0
      CALL integrScalar(ACoupled, temp, nsrflistCoupled, numCoupledSrfs)

      IF (myrank .EQ. master) THEN

c       CHECK AREA
        DO loopA = 1,numCoupledSrfs
          IF(ACoupled(loopA).LE.0.0D0)THEN
            PRINT *,''
            PRINT *,'WARNING: Zero Area for Coupled Surface ',loopA
            PRINT *,''
          ENDIF
        ENDDO

         IF (iGenFromFile .EQ. 1) THEN
            INQUIRE(FILE='GenBC',EXIST=ierr)
            IF (.NOT.ierr) THEN
               INQUIRE(FILE='../GenBC',EXIST=ierr)
               IF (ierr) THEN
                  CALL system('cp ../GenBC ./')
               ELSE
                  PRINT *, 'ERROR: No executable GenBC has been found'
                  PRINT *, 'Please, make sure that you had this',
     2               ' executable file inside runnig directory'
                  STOP
               END IF
            END IF
         END IF

         INQUIRE(FILE='InitialData',EXIST=ierr)
         IF (ierr) THEN
            PRINT *, ''
            PRINT *, 'NOTE: Initializing General BC form previous simulation'
            PRINT *, ''
         ELSE
            INQUIRE(FILE='../InitialData',EXIST=ierr)
            IF (ierr) THEN
               CALL system('cp ../InitialData ./')
               PRINT *, ''
               PRINT *,'NOTE: General BC has been initialized from provided',
     2            ' file'
               PRINT *, ''
            ELSE
               PRINT *, ''
               PRINT *, 'NOTE: Self GenBC initializiation'
               PRINT *, ''
            END IF
         END IF
      END IF

      RETURN
      END SUBROUTINE initGenBC

!> Main subroutine for gathering required information and
!! communicating with 0D to acquire necessary information

      SUBROUTINE calcGenBC (y, yold)

      USE GeneralBC
      INCLUDE "global.h"
      INCLUDE "mpif.h"
      INCLUDE "common_blocks/conpar.h"
      INCLUDE "common_blocks/workfc.h"
      INCLUDE "common_blocks/inpdat.h"
      INCLUDE "common_blocks/matdat.h"
      INCLUDE "common_blocks/nomodule.h"

      INTEGER MPIstat(MPI_STATUS_SIZE), i, ia

      REAL*8 y(nshg,ndof), yold(nshg,ndof), QInitial, en(nshg,ndof),
     2       enOfEle,rho

      REAL*8 testnorm
      INTEGER ierr

c     Get Density
      rho = datmat(1,1,1)

      DO ia=1, nshg
         enOfEle = yold(ia,4)
         DO i=1, nsd
           enOfEle = enOfEle + 5D-1*rho*yold(ia,i)*yold(ia,i)
         END DO
         DO i=1, nsd
            en(ia,i) = enOfEle*yold(ia,i)
         END DO
      END DO

      CALL GetFlowQ (QCoupled,  y,    nsrflistCoupled, numCoupledSrfs)
      CALL GetFlowQ (QnCoupled, yold, nsrflistCoupled, numCoupledSrfs)
      CALL GetFlowQ (ECoupled,  en,   nsrflistCoupled, numCoupledSrfs)

      CALL integrScalar (PCoupled, y(:,4),     nsrflistCoupled,
     2   numCoupledSrfs)

      CALL integrScalar (PnCoupled, yold(:,4), nsrflistCoupled,
     2   numCoupledSrfs)

      DO i=1, numCoupledSrfs
         PnCoupled(i) = PnCoupled(i)/ACoupled(i)
         PCoupled(i) = PCoupled(i)/ACoupled(i)
      END DO

      IF (myrank .EQ. master) THEN
         IF (GenFlag .EQ. 'L') THEN
            i = numCoupledSrfs
            CALL printGen (i, nsrflistCoupled(1:i), QCoupled(1:i),
     2         PCoupled(1:i), ECoupled(1:i))
         END IF

         OPEN (1,FILE='GenBC.int', STATUS='UNKNOWN', FORM='UNFORMATTED')
         WRITE (1) GenFlag
         WRITE (1) Delt(1)
         WRITE (1) numDirichletSrfs
         WRITE (1) numNeumannSrfs
         DO i=1, numCoupledSrfs
            IF (i .LE. numDirichletSrfs) THEN
               WRITE (1) PnCoupled(i), PCoupled(i)
            ELSE
               WRITE (1) QnCoupled(i), QCoupled(i)
            END IF
         END DO
         CLOSE (1)

c        EXTERNAL CALL TO GENBC
         IF (iGenFromFile .EQ. 1) THEN
            CALL system('./GenBC')
         ELSE
            CALL system('GenBC')
         END IF

         OPEN (1,FILE='GenBC.int',STATUS='OLD', FORM='UNFORMATTED')
         DO i=1, numCoupledSrfs
            IF (i .LE. numDirichletSrfs) THEN
               READ (1) QCoupled(i)
            ELSE
               READ (1) PCoupled(i)
            END IF
         END DO
         CLOSE(1)
      END IF

      i = MAXSURF + 1
      CALL MPI_BCAST(QCoupled, i, MPI_DOUBLE_PRECISION, master,
     2   MPI_COMM_WORLD, ierr)

      CALL MPI_BCAST(PCoupled, i, MPI_DOUBLE_PRECISION, master,
     2   MPI_COMM_WORLD, ierr)

      RETURN
      END SUBROUTINE calcGenBC

!> Main subroutine for calculating the general surface pressure
!! derivative with respect to flow rate

      SUBROUTINE calcGenBCDerivative

      USE GeneralBC
      INCLUDE "global.h"
      INCLUDE "mpif.h"
      INCLUDE "common_blocks/workfc.h"
      INCLUDE "common_blocks/inpdat.h"
      INCLUDE "common_blocks/nomodule.h"

      INTEGER MPIstat(MPI_STATUS_SIZE), i0
      REAL*8 diff, PBase(numCoupledSrfs), garbage
      REAL*8, PARAMETER :: absTol = 1D-8, relTol = 1D-5
      INTEGER ierr,i,j

      IF (numNeumannSrfs .EQ. 0) RETURN

      PGenDerFlag = .FALSE.
      i0 = numDirichletSrfs

      diff = SUM(ABS(QCoupled(i0+1:numCoupledSrfs)))
     2   /REAL(numNeumannSrfs,8)

      IF (diff*relTol .LT. absTol) THEN
         diff = absTol
      ELSE
         diff = diff*relTol
      END IF

      IF (myrank .EQ. master) THEN
         DO j=i0, numCoupledSrfs
            OPEN (1, FILE='GenBC.int', STATUS='UNKNOWN',
     2         FORM='UNFORMATTED')
            WRITE (1) 'D'
            WRITE (1) Delt(1)
            WRITE (1) numDirichletSrfs
            WRITE (1) numNeumannSrfs
            DO i=1, numCoupledSrfs
               IF (i .LE. i0) THEN
                  WRITE (1) PnCoupled(i), PCoupled(i)
               ELSE
                  IF (i .NE. j) THEN
                     WRITE (1) QnCoupled(i), QCoupled(i)
                  ELSE
                     WRITE (1) QnCoupled(i), QCoupled(i) + diff
                  END IF
               END IF
            END DO
            CLOSE (1)

            IF (iGenFromFile .EQ. 1) THEN
               CALL system('./GenBC')
            ELSE
               CALL system('GenBC')
            END IF

            OPEN (1,FILE='GenBC.int',STATUS='OLD', FORM='UNFORMATTED')
            DO i=1, numCoupledSrfs
               IF (i .LE. i0) THEN
                  READ (1) garbage
               ELSE
                  IF (i .NE. j) THEN
                     IF (j .EQ. i0) THEN
                        READ (1) PBase(i)
                     ELSE
                        READ (1) garbage
                     END IF
                  ELSE
                     READ (1) PGenDer(i)
                     PGenDer(i) = (PGenDer(i) - PBase(i))/diff
                  END IF
               END IF
            END DO
            CLOSE(1)
         END DO
      END IF

      i = MAXSURF + 1
      CALL MPI_BCAST(PGenDer, i, MPI_DOUBLE_PRECISION, master,
     2   MPI_COMM_WORLD, ierr)

      RETURN
      END SUBROUTINE calcGenBCDerivative

!> Writting results to the disk

      SUBROUTINE printGen (nSrfs, surfID, Q, P, E)

      IMPLICIT NONE

      LOGICAL ierr

      INTEGER nSrfs, nParam, i
      PARAMETER (nParam = 3)
      INTEGER surfID(nSrfs)

      REAL*8 Q(nSrfs), P(nSrfs), E(nSrfs), R(nParam,nSrfs)

      CHARACTER*32 fileNames(nParam)
      CHARACTER(len=32) :: myFMT1,myFMT2

      IF (nSrfs .EQ. 0) RETURN

      fileNames = (/'QGeneral','PGeneral','EGeneral'/)

      R(1,:) = Q
      R(2,:) = P
      R(3,:) = E

c     SET FORMATS
      WRITE(myFMT1, '("(",I0,"(E13.5))")') nSrfs
      WRITE(myFMT2, '("(",I0,"(I13))")') nSrfs

      DO i=1, nParam
         INQUIRE(FILE=TRIM(fileNames(i)), EXIST=ierr)
         IF (ierr) THEN
            OPEN (1, FILE=TRIM(fileNames(i)), STATUS='OLD',
     2         ACCESS='APPEND')

            WRITE (1,fmt=myFMT1) R(i,:)
            CLOSE (1)
         ELSE
            OPEN (1, FILE=TRIM(fileNames(i)), STATUS='NEW')
            WRITE (1,fmt=myFMT2) surfID
            WRITE (1,fmt=myFMT1) R(i,:)
            CLOSE (1)
         END IF
      END DO

      RETURN
      END SUBROUTINE printGen

!> This subroutine applies the normal velocity constraint
!! for the selected surfaces

      SUBROUTINE normCon (iens, y, xKebe, xGoC, rl)

! This is how the xKebe is formed
!  1 2 3
!  4 5 6
!  7 8 9

      USE pvsQbi

      INCLUDE "global.h"
      INCLUDE "common_blocks/propar.h"
      INCLUDE "common_blocks/shpdat.h"
      INCLUDE "common_blocks/conpar.h"
      INCLUDE "common_blocks/elmpar.h"
      INCLUDE "common_blocks/nomodule.h"

      INTEGER iens(npro,nshape), ie, ia, ib, iAc, i, j, l, iSrf

      REAL*8 norm(nsd,nsd), ndy, y(nshg, ndofl),
     2   xKebe(npro,9,nshl,nshl),
     3   K(nsd,nsd), nK(nsd,nsd), nKnt(nsd,nsd),
     4   xGoC(npro,4,nshl,nshl), G(nsd), nG(nsd),
     5   rl(npro,nshl,nflow), R(nsd), nR(nsd)

      INTENT(IN) iens, y
      INTENT(INOUT) xKebe, xGoC, rl

c      return
      DO iSrf = 1, numNormalSrfs
         DO ie = 1, npro
            DO ia = 1, nshl
               iAc = ABS(iens(ie,ia))
               IF (ndsurf(iAc) .EQ. nsrflistNormal(iSrf)) THEN
                  CALL getNorm (iAc, norm)

                  ndy = 0D0
                  DO i = 1, nsd
                     ndy = ndy + norm(1,i)*y(iA,i)
                  END DO
                  IF (.true.)then!ndy .LT. 0D0) THEN
c********************************************************************
c First modifing RHS matrix
                     DO i = 1, nsd
                        R(i) = rl(ie,ia,i)
                     END DO
c Here setting nR = n*R
                     DO i=1, nsd
                        nR(i) = 0D0
                        DO j=1,nsd
                           nR(i) = nR(i) + norm(i,j)*R(j)
                        END DO
                     END DO

c Now setting second and third row equal to zero, since the normal to the face is
c in the first direction
                     DO i=2, nsd
                        nR(i) = 0D0
                     END DO

c Rotating back the coordinate R = n^t*nR
                     DO i=1, nsd
                        R(i) = 0D0
                        DO j=1,nsd
                           R(i) = R(i) + norm(j,i)*nR(j)
                        END DO
                     END DO

c Assigning the new residual vector
                     DO i = 1, nsd
                        rl(ie,ia,i) = R(i)
                     END DO

c********************************************************************
c Now modifing LHS matrices

                     DO ib = 1, nshl
                        DO i=1, nsd
                           G(i) = xGoC(ie,i,ia,ib)
                           DO j=1,nsd
                              K(i,j) = xKebe(ie,j+nsd*(i-1),ia,ib)
                           END DO
                        END DO
c
c Here setting nK = n*K and nG = n*G
c
                        DO i=1, nsd
                           nG(i) = 0D0
                           DO j=1,nsd
                              nG(i) = nG(i) + norm(i,j)*G(j)
                              nK(i,j) = 0D0
                              DO l=1, nsd
                                 nK(i,j) = nK(i,j) + norm(i,l)*K(l,j)
                              END DO
                           END DO
                        END DO
c
c Here setting nKnt = nK*n^t
c
                        DO i=1, nsd
                           DO j=1,nsd
                              nKnt(i,j) = 0D0
                              DO l=1, nsd
                                 nKnt(i,j) = nKnt(i,j) +
     2                              nK(i,l)*norm(j,l)
                              END DO
                           END DO
                        END DO
c
c Now setting second and third row equal to zero, since the normal to the face is
c in the first direction and we dont want to constrain the velocity in normal direction,
c first row should be conserved. Note that diagonal element is not chaged to prevent
c singularity of stiffness matrix
c
                        DO i=2, nsd
                           IF (ia .NE. ib) THEN
                              nG(i) = 0D0
                              DO j=1,nsd
                                 nKnt(i,j) = 0D0
                              END DO
                           ELSE
                              nG(i) = 0D0
                              DO j=1,nsd
                                 IF (i.NE.j) THEN
                                    nKnt(i,j) = 0D0
                                 END IF
                              END DO
                           END IF
                        END DO

c Rotating back the coordinate nK = nKnt*n and G = nt*nG

                        DO i=1, nsd
                           G(i) = 0D0
                           DO j=1,nsd
                              G(i) = G(i) + norm(j,i)*nG(j)
                              nK(i,j) = 0D0
                              DO l=1, nsd
                                 nK(i,j) = nK(i,j) + nKnt(i,l)*norm(l,j)
                              END DO
                           END DO
                        END DO
c
c Rotating back the coordinate K = n^t*nK
c
                        DO i=1, nsd
                           DO j=1,nsd
                              K(i,j) = 0D0
                              DO l=1, nsd
                                 K(i,j) = K(i,j) + norm(l,i)*nK(l,j)
                              END DO
                           END DO
                        END DO
c
c Assigning the new stiffness matrix
c
                        DO i=1, nsd
                           xGoC(ie, i, ia, ib) = G(i)
                           DO j=1,nsd
                              xKebe(ie, j + nsd*(i-1), ia, ib) = K(i,j)
                           END DO
                        END DO
                     END DO ! loop for another local node
                  END IF ! if clause for inward velocity
               END IF ! if clause for nodes inside the face
            END DO ! loop for the node
         END DO ! loop over the elements
      END DO ! loop over the faces

      RETURN

      END SUBROUTINE normCon

!> This subroutine is intended for calculating normal vector

      SUBROUTINE getNorm(iAc, norm)

      USE pvsQbi

      INCLUDE "global.h"

      INTEGER iAc, i
      REAL*8 normal(nsd), tang1(nsd), tang2(nsd), temp, norm(nsd, nsd)

      INTENT(IN) iAc
      INTENT(OUT) norm

      temp = 0D0
      DO i=1,nsd
         normal(i) = NABI(iAc,i)
         temp = temp + normal(i)**2D0
      END DO
      temp = SQRT(temp)
      normal = normal/temp

      IF (ABS(normal(1)) .NE. 1D0) THEN
         tang1(1) = 0D0
         tang1(2) = normal(3)
         tang1(3) = -normal(2)

         temp = 0D0
         DO i=1,nsd
            temp = temp + tang1(i)**2D0
         END DO
         temp = SQRT(temp)
         tang1 = tang1/temp
      ELSE
         tang1 = (/0D0,1D0,0D0/)
      END IF

      tang2(1) = normal(2)*tang1(3) - normal(3)*tang1(2)
      tang2(2) = normal(3)*tang1(1) - normal(1)*tang1(3)
      tang2(3) = normal(1)*tang1(2) - normal(2)*tang1(1)

      norm(1,:) = normal
      norm(2,:) = tang1
      norm(3,:) = tang2

      RETURN
      END SUBROUTINE getNorm

#endif


! { dg-do run { target fd_truncate } }
! { dg-options "-std=gnu" }
! Test case prepared by Jerry DeLisle  <jvdelisle@gcc.gnu.org>
real :: a(4), b(4)
real :: c
integer :: istat, j
character(25) :: msg

a = 23.45
b = 0.0
open(10, file='mydata_f2003_io_1', asynchronous="yes", blank="null")

write(10,'(10f8.3)', asynchronous="yes", decimal="comma", id=j) a
rewind(10)
read(10,'(10f8.3)', asynchronous="yes", decimal="comma", blank="zero") b
if (any(b.ne.23.45)) STOP 1

c = 3.14
write(msg, *, decimal="comma") c
if (msg(1:7).ne."   3,14") STOP 2

b = 0.0
rewind(10)
write(10,'(10f8.3)', asynchronous="yes", decimal="point") a
rewind(10)
read(10,'(10f8.3)', asynchronous="yes", decimal="point") b
if (any(b.ne.23.45)) STOP 3

wait(unit=10, err=25, iostat=istat, iomsg=msg, end=35, id=j)

! do some stuff with a
25 continue

35 continue

close(10, status="delete")
end


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

module hecmw_array_util
  use hecmw_util
  implicit none

  private
  public :: hecmw_qsort_int_array
  public :: hecmw_uniq_int_array
  public :: hecmw_bsearch_int_array

contains

  recursive subroutine hecmw_qsort_int_array(array, istart, iend)
    implicit none
    integer(kind=kint), intent(inout) :: array(:)
    integer(kind=kint), intent(in) :: istart, iend
    integer(kind=kint) :: left, right, center
    integer(kind=kint) :: pivot, tmp
    if (istart >= iend) return
    center = (istart + iend) / 2
    pivot = array(center)
    left = istart
    right = iend
    do
      do while (array(left) < pivot)
        left = left + 1
      end do
      do while (pivot < array(right))
        right = right - 1
      end do
      if (left >= right) exit
      tmp = array(left)
      array(left) = array(right)
      array(right) = tmp
      left = left + 1
      right = right - 1
    end do
    if (istart < left-1) call hecmw_qsort_int_array(array, istart, left-1)
    if (right+1 < iend) call hecmw_qsort_int_array(array, right+1, iend)
  end subroutine hecmw_qsort_int_array

  subroutine hecmw_uniq_int_array(array, istart, iend, ndup)
    implicit none
    integer(kind=kint), intent(inout) :: array(:)
    integer(kind=kint), intent(in) :: istart, iend
    integer(kind=kint), intent(out) :: ndup
    integer(kind=kint) :: i
    ndup = 0
    do i = istart+1, iend
      if (array(i) == array(i - 1 - ndup)) then
        ndup = ndup + 1
      else if (ndup > 0) then
        array(i - ndup) = array(i)
      endif
    end do
  end subroutine hecmw_uniq_int_array

  subroutine hecmw_bsearch_int_array(array, istart, iend, val, idx)
    implicit none
    integer(kind=kint), intent(in) :: array(:)
    integer(kind=kint), intent(in) :: istart, iend
    integer(kind=kint), intent(in) :: val
    integer(kind=kint), intent(out) :: idx
    integer(kind=kint) :: center, left, right, pivot
    left = istart
    right = iend
    do
      if (left > right) then
        idx = -1
        exit
      end if
      center = (left + right) / 2
      pivot = array(center)
      if (val < pivot) then
        right = center - 1
        cycle
      else if (pivot < val) then
        left = center + 1
        cycle
      else ! if (pivot == val) then
        idx = center
        exit
      end if
    end do
  end subroutine hecmw_bsearch_int_array

end module hecmw_array_util

max_stars_repo_name: fujitsu/FrontISTR

max_issues_repo_name: fujitsu/FrontISTR

max_forks_repo_name: fujitsu/FrontISTR


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 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 gen3dfrom1d(i,kon,ipkon,lakon,ne,iponor,xnor,knor,
     &  thicke,ntrans,inotr,trab,nk,nk_,co,offset,mi)
!
!     expands 1d element i into a 3d element
!
      implicit none
!
      character*8 lakon(*)
!
      integer mi(*),i,kon(*),ipkon(*),ne,iponor(2,*),knor(*),ntrans,
     &  inotr(2,*),nk,nk_,indexe,j,nodel(8),indexx,indexk,k,nodeb(8,3),
     &  nope,ishift
!
      real*8 xnor(*),thicke(mi(3),*),trab(7,*),co(3,*),offset(2,*),
     &  thickb(2,3),xnorb(6,3),sc
!
      indexe=ipkon(i)
!
!     check whether linear or quadratic
!
      if((lakon(i)(3:3).eq.'1').or.(lakon(i)(4:4).eq.'2')) then
         nope=2
         if(lakon(i)(4:4).eq.'R') then
            ishift=8
         else
            ishift=11
         endif
      elseif((lakon(i)(3:3).eq.'2').or.(lakon(i)(4:5).eq.'3')) then
         nope=3
         ishift=20
      endif
!
!     localizing the nodes, thicknesses and normals for the
!     beam element
!
c      do j=1,3
      do j=1,nope
         nodel(j)=kon(indexe+j)
         kon(indexe+ishift+j)=nodel(j)
         indexx=iponor(1,indexe+j)
         indexk=iponor(2,indexe+j)
         thickb(1,j)=thicke(1,indexe+j)
         thickb(2,j)=thicke(2,indexe+j)
         do k=1,6
            xnorb(k,j)=xnor(indexx+k)
         enddo
         do k=1,8
            nodeb(k,j)=knor(indexk+k)
         enddo
         if(ntrans.gt.0) then
            do k=1,8
               inotr(1,nodeb(k,j))=inotr(1,nodel(j))
            enddo
         endif
      enddo
!
!     generating the 3-D element topology for beam elements
!
!     rectangular cross section (or parent section)
!
      if(lakon(i)(8:8).ne.'C') then
         kon(indexe+1)=nodeb(1,1)
         do j=1,3
            co(j,nodeb(1,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &           +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
         enddo
         kon(indexe+2)=nodeb(1,nope)
         do j=1,3
            co(j,nodeb(1,nope))=co(j,nodel(nope))
     &           -thickb(1,nope)*xnorb(j,nope)*(.5d0+offset(1,i))
     &           +thickb(2,nope)*xnorb(j+3,nope)*(.5d0-offset(2,i))
         enddo
         kon(indexe+3)=nodeb(2,nope)
         do j=1,3
            co(j,nodeb(2,nope))=co(j,nodel(nope))
     &           -thickb(1,nope)*xnorb(j,nope)*(.5d0+offset(1,i))
     &           -thickb(2,nope)*xnorb(j+3,nope)*(.5d0+offset(2,i))
         enddo
         kon(indexe+4)=nodeb(2,1)
         do j=1,3
            co(j,nodeb(2,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
         enddo
         kon(indexe+5)=nodeb(4,1)
         do j=1,3
            co(j,nodeb(4,1))=co(j,nodel(1))
     &           +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &           +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
         enddo
         kon(indexe+6)=nodeb(4,nope)
         do j=1,3
            co(j,nodeb(4,nope))=co(j,nodel(nope))
     &           +thickb(1,nope)*xnorb(j,nope)*(.5d0-offset(1,i))
     &           +thickb(2,nope)*xnorb(j+3,nope)*(.5d0-offset(2,i))
         enddo
         kon(indexe+7)=nodeb(3,nope)
         do j=1,3
            co(j,nodeb(3,nope))=co(j,nodel(nope))
     &           +thickb(1,nope)*xnorb(j,nope)*(.5d0-offset(1,i))
     &           -thickb(2,nope)*xnorb(j+3,nope)*(.5d0+offset(2,i))
         enddo
         kon(indexe+8)=nodeb(3,1)
         do j=1,3
            co(j,nodeb(3,1))=co(j,nodel(1))
     &           +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
         enddo
!
!        middle nodes for quadratic elements
!
         if(nope.eq.3) then
            kon(indexe+9)=nodeb(1,2)
            do j=1,3
               co(j,nodeb(1,2))=co(j,nodel(2))
     &              -thickb(1,2)*xnorb(j,2)*(.5d0+offset(1,i))
     &              +thickb(2,2)*xnorb(j+3,2)*(.5d0-offset(2,i))
            enddo
            kon(indexe+10)=nodeb(5,3)
            do j=1,3
               co(j,nodeb(5,3))=co(j,nodel(3))
     &              -thickb(1,3)*xnorb(j,3)*(.5d0+offset(1,i))
     &              -thickb(2,3)*xnorb(j+3,3)*offset(2,i)
            enddo
            kon(indexe+11)=nodeb(2,2)
            do j=1,3
               co(j,nodeb(2,2))=co(j,nodel(2))
     &              -thickb(1,2)*xnorb(j,2)*(.5d0+offset(1,i))
     &              -thickb(2,2)*xnorb(j+3,2)*(.5d0+offset(2,i))
            enddo
            kon(indexe+12)=nodeb(5,1)
            do j=1,3
               co(j,nodeb(5,1))=co(j,nodel(1))
     &              -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &              -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
            enddo
            kon(indexe+13)=nodeb(4,2)
            do j=1,3
               co(j,nodeb(4,2))=co(j,nodel(2))
     &              +thickb(1,2)*xnorb(j,2)*(.5d0-offset(1,i))
     &              +thickb(2,2)*xnorb(j+3,2)*(.5d0-offset(2,i))
            enddo
            kon(indexe+14)=nodeb(7,3)
            do j=1,3
               co(j,nodeb(7,3))=co(j,nodel(3))
     &              +thickb(1,3)*xnorb(j,3)*(.5d0-offset(1,i))
     &              -thickb(2,3)*xnorb(j+3,3)*offset(2,i)
            enddo
            kon(indexe+15)=nodeb(3,2)
            do j=1,3
               co(j,nodeb(3,2))=co(j,nodel(2))
     &              +thickb(1,2)*xnorb(j,2)*(.5d0-offset(1,i))
     &              -thickb(2,2)*xnorb(j+3,2)*(.5d0+offset(2,i))
            enddo
            kon(indexe+16)=nodeb(7,1)
            do j=1,3
               co(j,nodeb(7,1))=co(j,nodel(1))
     &              +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &              -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
            enddo
            kon(indexe+17)=nodeb(8,1)
            do j=1,3
               co(j,nodeb(8,1))=co(j,nodel(1))
     &              -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &              +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
            enddo
            kon(indexe+18)=nodeb(8,3)
            do j=1,3
               co(j,nodeb(8,3))=co(j,nodel(3))
     &              -thickb(1,3)*xnorb(j,3)*offset(1,i)
     &              +thickb(2,3)*xnorb(j+3,3)*(.5d0-offset(2,i))
            enddo
            kon(indexe+19)=nodeb(6,3)
            do j=1,3
               co(j,nodeb(6,3))=co(j,nodel(3))
     &              -thickb(1,3)*xnorb(j,3)*offset(1,i)
     &              -thickb(2,3)*xnorb(j+3,3)*(.5d0+offset(2,i))
            enddo
            kon(indexe+20)=nodeb(6,1)
            do j=1,3
               co(j,nodeb(6,1))=co(j,nodel(1))
     &              -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &              -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
            enddo
!
!           generating coordinates for the expanded nodes which
!           are not used by the C3D20(R) element (needed for the
!           determination of the knot dimension)
!
            do j=1,3
               co(j,nodeb(5,2))=co(j,nodel(2))
     &              -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &              -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
            enddo
            do j=1,3
               co(j,nodeb(7,2))=co(j,nodel(2))
     &              +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &              -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
            enddo
            do j=1,3
               co(j,nodeb(8,2))=co(j,nodel(2))
     &              -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &              +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
            enddo
            do j=1,3
               co(j,nodeb(6,2))=co(j,nodel(2))
     &              -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &              -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
            enddo
         endif
      else
!
!                 circular cross section
!
         if(nope.eq.2) then
            write(*,*) '*ERROR in gen3dfrom1d: element ',i,
     &        'is a linear beam element with circular cross section'
            write(*,*) '       Please use quadratic elements for beams
     &with circular cross section.'
            call exit(201)
         endif
!
         sc=.5d0/dsqrt(2.d0)
         kon(indexe+1)=nodeb(1,1)
         do j=1,3
            co(j,nodeb(1,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*(sc+offset(1,i))
     &           +thickb(2,1)*xnorb(j+3,1)*(sc-offset(2,i))
         enddo
         kon(indexe+2)=nodeb(1,3)
         do j=1,3
            co(j,nodeb(1,3))=co(j,nodel(3))
     &           -thickb(1,3)*xnorb(j,3)*(sc+offset(1,i))
     &           +thickb(2,3)*xnorb(j+3,3)*(sc-offset(2,i))
         enddo
         kon(indexe+3)=nodeb(2,3)
         do j=1,3
            co(j,nodeb(2,3))=co(j,nodel(3))
     &           -thickb(1,3)*xnorb(j,3)*(sc+offset(1,i))
     &           -thickb(2,3)*xnorb(j+3,3)*(sc+offset(2,i))
         enddo
         kon(indexe+4)=nodeb(2,1)
         do j=1,3
            co(j,nodeb(2,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*(sc+offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*(sc+offset(2,i))
         enddo
         kon(indexe+5)=nodeb(4,1)
         do j=1,3
            co(j,nodeb(4,1))=co(j,nodel(1))
     &           +thickb(1,1)*xnorb(j,1)*(sc-offset(1,i))
     &           +thickb(2,1)*xnorb(j+3,1)*(sc-offset(2,i))
         enddo
         kon(indexe+6)=nodeb(4,3)
         do j=1,3
            co(j,nodeb(4,3))=co(j,nodel(3))
     &           +thickb(1,3)*xnorb(j,3)*(sc-offset(1,i))
     &           +thickb(2,3)*xnorb(j+3,3)*(sc-offset(2,i))
         enddo
         kon(indexe+7)=nodeb(3,3)
         do j=1,3
            co(j,nodeb(3,3))=co(j,nodel(3))
     &           +thickb(1,3)*xnorb(j,3)*(sc-offset(1,i))
     &           -thickb(2,3)*xnorb(j+3,3)*(sc+offset(2,i))
         enddo
         kon(indexe+8)=nodeb(3,1)
         do j=1,3
            co(j,nodeb(3,1))=co(j,nodel(1))
     &           +thickb(1,1)*xnorb(j,1)*(sc-offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*(sc+offset(2,i))
         enddo
         kon(indexe+9)=nodeb(1,2)
         do j=1,3
            co(j,nodeb(1,2))=co(j,nodel(2))
     &           -thickb(1,2)*xnorb(j,2)*(sc+offset(1,i))
     &           +thickb(2,2)*xnorb(j+3,2)*(sc-offset(2,i))
         enddo
         kon(indexe+10)=nodeb(5,3)
         do j=1,3
            co(j,nodeb(5,3))=co(j,nodel(3))
     &           -thickb(1,3)*xnorb(j,3)*(.5d0+offset(1,i))
     &           -thickb(2,3)*xnorb(j+3,3)*offset(2,i)
         enddo
         kon(indexe+11)=nodeb(2,2)
         do j=1,3
            co(j,nodeb(2,2))=co(j,nodel(2))
     &           -thickb(1,2)*xnorb(j,2)*(sc+offset(1,i))
     &           -thickb(2,2)*xnorb(j+3,2)*(sc+offset(2,i))
         enddo
         kon(indexe+12)=nodeb(5,1)
         do j=1,3
            co(j,nodeb(5,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
         enddo
         kon(indexe+13)=nodeb(4,2)
         do j=1,3
            co(j,nodeb(4,2))=co(j,nodel(2))
     &           +thickb(1,2)*xnorb(j,2)*(sc-offset(1,i))
     &           +thickb(2,2)*xnorb(j+3,2)*(sc-offset(2,i))
         enddo
         kon(indexe+14)=nodeb(7,3)
         do j=1,3
            co(j,nodeb(7,3))=co(j,nodel(3))
     &           +thickb(1,3)*xnorb(j,3)*(.5d0-offset(1,i))
     &           -thickb(2,3)*xnorb(j+3,3)*offset(2,i)
         enddo
         kon(indexe+15)=nodeb(3,2)
         do j=1,3
            co(j,nodeb(3,2))=co(j,nodel(2))
     &           +thickb(1,2)*xnorb(j,2)*(sc-offset(1,i))
     &           -thickb(2,2)*xnorb(j+3,2)*(sc+offset(2,i))
         enddo
         kon(indexe+16)=nodeb(7,1)
         do j=1,3
            co(j,nodeb(7,1))=co(j,nodel(1))
     &           +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
         enddo
         kon(indexe+17)=nodeb(8,1)
         do j=1,3
            co(j,nodeb(8,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &           +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
         enddo
         kon(indexe+18)=nodeb(8,3)
         do j=1,3
            co(j,nodeb(8,3))=co(j,nodel(3))
     &           -thickb(1,3)*xnorb(j,3)*offset(1,i)
     &           +thickb(2,3)*xnorb(j+3,3)*(.5d0-offset(2,i))
         enddo
         kon(indexe+19)=nodeb(6,3)
         do j=1,3
            co(j,nodeb(6,3))=co(j,nodel(3))
     &           -thickb(1,3)*xnorb(j,3)*offset(1,i)
     &           -thickb(2,3)*xnorb(j+3,3)*(.5d0+offset(2,i))
         enddo
         kon(indexe+20)=nodeb(6,1)
         do j=1,3
            co(j,nodeb(6,1))=co(j,nodel(1))
     &           -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &           -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
         enddo
!
!           generating coordinates for the expanded nodes which
!           are not used by the C3D20(R) element (needed for the
!           determination of the knot dimension)
!
         do j=1,3
            co(j,nodeb(5,2))=co(j,nodel(2))
     &           -thickb(1,1)*xnorb(j,1)*(.5d0+offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
         enddo
         do j=1,3
            co(j,nodeb(7,2))=co(j,nodel(2))
     &           +thickb(1,1)*xnorb(j,1)*(.5d0-offset(1,i))
     &           -thickb(2,1)*xnorb(j+3,1)*offset(2,i)
         enddo
         do j=1,3
            co(j,nodeb(8,2))=co(j,nodel(2))
     &           -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &           +thickb(2,1)*xnorb(j+3,1)*(.5d0-offset(2,i))
         enddo
         do j=1,3
            co(j,nodeb(6,2))=co(j,nodel(2))
     &           -thickb(1,1)*xnorb(j,1)*offset(1,i)
     &           -thickb(2,1)*xnorb(j+3,1)*(.5d0+offset(2,i))
         enddo
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop