!moduleMesh1DRad: 1D radial module
!              x == r
!              y == theta (unused)
!              z == unused
MODULE moduleMesh1DRad
  USE moduleMesh
  USE moduleMeshBoundary
  IMPLICIT NONE

  REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
  REAL(8), PARAMETER:: wSeg(1:3)  = (/        5.D0/9.D0 , 8.D0/9.D0,        5.D0/9.D0  /)

  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
    !Element coordinates
    REAL(8):: r = 0.D0
    CONTAINS
      PROCEDURE, PASS:: init => initNode1DRad
      PROCEDURE, PASS:: getCoordinates => getCoord1DRad

  END TYPE meshNode1DRad

  TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge1DRad
    !Element coordinates
    REAL(8):: r = 0.D0
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL()
    CONTAINS
      PROCEDURE, PASS:: init         => initEdge1DRad
      PROCEDURE, PASS:: getNodes     => getNodes1DRad
      PROCEDURE, PASS:: intersection => intersection1DRad
      PROCEDURE, PASS:: randPos      => randPos1DRad

  END TYPE meshEdge1DRad

  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad
    CONTAINS
      PROCEDURE, PASS:: detJac => detJ1DRad
      PROCEDURE, PASS:: invJac => invJ1DRad
      PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
      PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
  
  END TYPE meshVol1DRad

  ABSTRACT INTERFACE
    PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
      REAL(8), INTENT(in):: xi(1:3)
      REAL(8), ALLOCATABLE:: dPsi(:,:)

    END FUNCTION dPsi_interface

    PURE SUBROUTINE partialDer_interface(self, dPsi, dx)
      IMPORT meshVol1DRad

      CLASS(meshVol1DRad), INTENT(in):: self
      REAL(8), INTENT(in):: dPsi(1:,1:)
      REAL(8), INTENT(out), DIMENSION(1):: dx

    END SUBROUTINE partialDer_interface

  END INTERFACE

  TYPE, PUBLIC, EXTENDS(meshVol1DRad):: meshVol1DRadSegm
    !Element coordinates
    REAL(8):: r(1:2)
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
    !Connectivity to adjacent elements
    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
    REAL(8):: arNodes(1:2)
    CONTAINS
      PROCEDURE, PASS::   init        => initVol1DRadSegm
      PROCEDURE, PASS::   randPos     => randPos1DRadSeg
      PROCEDURE, PASS::   area        => areaRad
      PROCEDURE, NOPASS:: fPsi        => fPsiRad
      PROCEDURE, NOPASS:: dPsi        => dPsiRad
      PROCEDURE, PASS::   partialDer  => partialDerRad
      PROCEDURE, PASS::   elemK       => elemKRad
      PROCEDURE, PASS::   elemF       => elemFRad
      PROCEDURE, NOPASS:: inside      => insideRad
      PROCEDURE, PASS::   gatherEF    => gatherEFRad
      PROCEDURE, PASS::   gatherMF    => gatherMFRad
      PROCEDURE, PASS::   getNodes    => getNodesRad
      PROCEDURE, PASS::   phy2log     => phy2logRad
      PROCEDURE, PASS::   nextElement => nextElementRad

  END TYPE meshVol1DRadSegm

  CONTAINS
    !NODE FUNCTIONS
    !Init node element
    SUBROUTINE initNode1DRad(self, n, r)
      USE moduleSpecies
      USE moduleRefParam
      USE OMP_LIB
      IMPLICIT NONE

      CLASS(meshNode1DRad), INTENT(out):: self
      INTEGER, INTENT(in):: n
      REAL(8), INTENT(in):: r(1:3)

      self%n = n
      self%r = r(1)/L_ref
      !Node volume, to be determined in mesh
      self%v = 0.D0

      !Allocates output
      ALLOCATE(self%output(1:nSpecies))

      CALL OMP_INIT_LOCK(self%lock)

    END SUBROUTINE initNode1DRad

    PURE FUNCTION getCoord1DRad(self) RESULT(r)
      IMPLICIT NONE

      CLASS(meshNode1DRad), INTENT(in):: self
      REAL(8):: r(1:3)

      r = (/ self%r, 0.D0, 0.D0 /)

    END FUNCTION getCoord1DRad

    !EDGE FUNCTIONS
    !Inits edge element
    SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface)
      USE moduleSpecies
      USE moduleBoundary
      USE moduleErrors
      IMPLICIT NONE

      CLASS(meshEdge1DRad), INTENT(out):: self
      INTEGER, INTENT(in):: n
      INTEGER, INTENT(in):: p(:)
      INTEGER, INTENT(in):: bt
      INTEGER, INTENT(in):: physicalSurface
      REAL(8), DIMENSION(1:3):: r1
      INTEGER:: s

      self%n = n
      self%n1 => mesh%nodes(p(1))%obj
      !Get element coordinates
      r1 = self%n1%getCoordinates()

      self%r = r1(1)

      self%normal = (/ 1.D0, 0.D0, 0.D0 /)
      self%normal = self%normal/NORM2(self%normal)

      !Boundary index
      self%boundary => boundary(bt)
      ALLOCATE(self%fboundary(1:nSpecies))
      !Assign functions to boundary
      DO s = 1, nSpecies
        CALL pointBoundaryFunction(self, s)

      END DO

      !Physical Surface
      self%physicalSurface = physicalSurface

    END SUBROUTINE initEdge1DRad

    !Get nodes from edge
    PURE FUNCTION getNodes1DRad(self) RESULT(n)
      IMPLICIT NONE

      CLASS(meshEdge1DRad), INTENT(in):: self
      INTEGER, ALLOCATABLE:: n(:)

      ALLOCATE(n(1))
      n = (/ self%n1%n /)

    END FUNCTION getNodes1DRad

    PURE FUNCTION intersection1DRad(self, r0) RESULT(r)
      IMPLICIT NONE

      CLASS(meshEdge1DRad), INTENT(in):: self
      REAL(8), DIMENSION(1:3), INTENT(in):: r0
      REAL(8), DIMENSION(1:3):: r

      r = (/ self%r, 0.D0, 0.D0 /)

    END FUNCTION intersection1DRad

    !Calculates a 'random' position in edge
    FUNCTION randPos1DRad(self) RESULT(r)
      CLASS(meshEdge1DRad), INTENT(in):: self
      REAL(8):: r(1:3)

      r = (/ self%r, 0.D0, 0.D0 /)

    END FUNCTION randPos1DRad

    !VOLUME FUNCTIONS
    !SEGMENT FUNCTIONS
    !Init segment element
    SUBROUTINE initVol1DRadSegm(self, n, p, nodes)
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(out):: self
      INTEGER, INTENT(in):: n
      INTEGER, INTENT(in):: p(:)
      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
      REAL(8), DIMENSION(1:3):: r1, r2

      self%n = n
      self%n1 => nodes(p(1))%obj
      self%n2 => nodes(p(2))%obj
      !Get element coordinates
      r1 = self%n1%getCoordinates()
      r2 = self%n2%getCoordinates()
      self%r = (/ r1(1), r2(1) /)

      !Assign node volume
      CALL self%area()
      self%n1%v = self%n1%v + self%arNodes(1)
      self%n2%v = self%n2%v + self%arNodes(2)

      CALL OMP_INIT_LOCK(self%lock)

      ALLOCATE(self%listPart_in(1:nSpecies))
      ALLOCATE(self%totalWeight(1:nSpecies))

    END SUBROUTINE initVol1DRadSegm

    !Calculates a random position in 1D volume
    FUNCTION randPos1DRadSeg(self) RESULT(r)
      USE moduleRandom
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8):: r(1:3)
      REAL(8):: xii(1:3)
      REAL(8), ALLOCATABLE:: fPsi(:)

      xii(1)   = random(-1.D0, 1.D0)
      xii(2:3) = 0.D0

      fPsi = self%fPsi(xii)
      r(1) = DOT_PRODUCT(fPsi, self%r)

    END FUNCTION randPos1DRadSeg

    !Computes element area
    PURE SUBROUTINE areaRad(self)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(inout):: self
      REAL(8):: l !element length
      REAL(8):: fPsi(1:2)
      REAL(8):: r
      REAL(8):: detJ
      REAL(8):: Xii(1:3)

      self%volume  = 0.D0
      self%arNodes = 0.D0
      !1 point Gauss integral
      Xii = 0.D0
      fPsi = self%fPsi(Xii)
      detJ = self%detJac(Xii)
      !Computes total volume of the cell
      r    = DOT_PRODUCT(fPsi, self%r)
      l = 2.D0*detJ
      self%volume  = r*l
      !Computes volume per node
      Xii = (/-5.D-1,  0.D0, 0.D0/)
      r = DOT_PRODUCT(self%fPsi(Xii),self%r)
      self%arNodes(1) = fPsi(1)*r*l
      Xii = (/ 5.D-1,  0.D0, 0.D0/)
      r = DOT_PRODUCT(self%fPsi(Xii),self%r)
      self%arNodes(2) = fPsi(2)*r*l

    END SUBROUTINE areaRad

    !Computes element functions at point xii
    PURE FUNCTION fPsiRad(xi) RESULT(fPsi)
      IMPLICIT NONE

      REAL(8), INTENT(in):: xi(1:3)
      REAL(8), ALLOCATABLE:: fPsi(:)

      ALLOCATE(fPsi(1:2))

      fPsi(1) = 1.D0 - xi(1)
      fPsi(2) = 1.D0 + xi(1)
      fPsi    = fPsi * 5.D-1

    END FUNCTION fPsiRad

    !Computes element derivative shape function at Xii
    PURE FUNCTION dPsiRad(xi) RESULT(dPsi)
      IMPLICIT NONE

      REAL(8), INTENT(in):: xi(1:3)
      REAL(8), ALLOCATABLE:: dPsi(:,:)

      ALLOCATE(dPsi(1:1, 1:2))

      dPsi(1, 1) = -5.D-1
      dPsi(1, 2) =  5.D-1

    END FUNCTION dPsiRad

    !Computes partial derivatives of coordinates
    PURE SUBROUTINE partialDerRad(self, dPsi, dx)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: dPsi(1:,1:)
      REAL(8), INTENT(out), DIMENSION(1):: dx

      dx(1) = DOT_PRODUCT(dPsi(1,:), self%r)

    END SUBROUTINE partialDerRad

    !Computes local stiffness matrix
    PURE FUNCTION elemKRad(self) RESULT(localK)
      USE moduleConstParam, ONLY: PI2
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), ALLOCATABLE:: localK(:,:)
      REAL(8):: Xii(1:3)
      REAL(8):: dPsi(1:1, 1:2)
      REAL(8):: invJ(1), detJ
      REAL(8):: r, fPsi(1:2)
      INTEGER:: l

      ALLOCATE(localK(1:2, 1:2))
      localK  = 0.D0
      Xii     = 0.D0
      DO l = 1, 3
        xii(1) = corSeg(l)
        dPsi    = self%dPsi(Xii)
        detJ    = self%detJac(Xii, dPsi)
        invJ    = self%invJac(Xii, dPsi)
        fPsi    = self%fPsi(Xii)
        r       = DOT_PRODUCT(fPsi, self%r)
        localK  = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
                                  RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
                           r*wSeg(l)/detJ

      END DO

      localK = localK*PI2

    END FUNCTION elemKRad

    PURE FUNCTION elemFRad(self, source) RESULT(localF)
      USE moduleConstParam, ONLY: PI2
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: source(1:)
      REAL(8), ALLOCATABLE:: localF(:)
      REAL(8):: fPsi(1:2)
      REAL(8):: detJ, f, r
      REAL(8):: Xii(1:3)
      INTEGER:: l

      ALLOCATE(localF(1:2))
      localF = 0.D0
      Xii    = 0.D0

      DO l = 1, 3
        xii(1) = corSeg(l)
        detJ   = self%detJac(Xii)
        fPsi   = self%fPsi(Xii)
        r      = DOT_PRODUCT(fPsi, self%r)
        f      = DOT_PRODUCT(fPsi, source)
        localF = localF + f*fPsi*r*wSeg(l)*detJ

      END DO

    END FUNCTION elemFRad

    PURE FUNCTION insideRad(xi) RESULT(ins)
      IMPLICIT NONE

      REAL(8), INTENT(in):: xi(1:3)
      LOGICAL:: ins

      ins = xi(1) >=-1.D0 .AND. &
            xi(1) <= 1.D0

    END FUNCTION insideRad

    !Gathers EF at position Xii
    PURE FUNCTION gatherEFRad(self, xi) RESULT(EF)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      REAL(8):: dPsi(1, 1:2)
      REAL(8):: phi(1:2)
      REAL(8):: EF(1:3)
      REAL(8):: invJ

      phi = (/ self%n1%emData%phi, &
               self%n2%emData%phi /)

      dPsi  =  self%dPsi(xi)
      invJ  =  self%invJac(xi, dPsi)
      EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
      EF(2) =  0.D0
      EF(3) =  0.D0

    END FUNCTION gatherEFRad

    PURE FUNCTION gatherMFRad(self, xi) RESULT(MF)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      REAL(8):: fPsi(1:2)
      REAL(8):: MF_Nodes(1:2, 1:3)
      REAL(8):: MF(1:3)
      REAL(8):: invJ

      MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
                           self%n2%emData%B(1) /)
      MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
                           self%n2%emData%B(2) /)
      MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
                           self%n2%emData%B(3) /)

      fPsi  = self%fPsi(xi)
      MF = MATMUL(fPsi, MF_Nodes)

    END FUNCTION gatherMFRad

    !Get nodes from 1D volume
    PURE FUNCTION getNodesRad(self) RESULT(n)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      INTEGER, ALLOCATABLE:: n(:)

      ALLOCATE(n(1:2))
      n = (/ self%n1%n, self%n2%n /)

    END FUNCTION getNodesRad

    PURE FUNCTION phy2logRad(self, r) RESULT(xN)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: xN(1:3)

      xN = 0.D0
      xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0

    END FUNCTION phy2logRad

    !Get next element for a logical position xi
    SUBROUTINE nextElementRad(self, xi, nextElement)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      CLASS(meshElement), POINTER, INTENT(out):: nextElement
      
      NULLIFY(nextElement)
      IF     (xi(1) < -1.D0) THEN
        nextElement => self%e2

      ELSEIF (xi(1) >  1.D0) THEN
        nextElement => self%e1

      END IF

    END SUBROUTINE nextElementRad

    !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
    !Computes the element Jacobian determinant
    PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ)
      IMPLICIT NONE

      CLASS(meshVol1DRad), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
      REAL(8), ALLOCATABLE:: dPsi(:,:)
      REAL(8):: dJ
      REAL(8):: dx(1)

      IF (PRESENT(dPsi_in)) THEN
        dPsi = dPsi_in

      ELSE
        dPsi = self%dPsi(xi)

      END IF

      CALL self%partialDer(dPsi, dx)
      dJ = dx(1)

    END FUNCTION detJ1DRad

    !Computes the invers Jacobian
    PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ)
      IMPLICIT NONE

      CLASS(meshVol1DRad), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
      REAL(8), ALLOCATABLE:: dPsi(:,:)
      REAL(8):: dx(1)
      REAL(8):: invJ

      IF (PRESENT(dPsi_in)) THEN
        dPsi = dPsi_in

      ELSE
        dPsi = self%dPsi(xi)

      END IF

      CALL self%partialDer(dPsi, dx)
      invJ = 1.D0/dx(1)

    END FUNCTION invJ1DRad

    SUBROUTINE connectMesh1DRad(self)
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(inout):: self
      INTEGER:: e, et

      DO e = 1, self%numVols
        !Connect Vol-Vol
        DO et = 1, self%numVols
          IF (e /= et) THEN
            CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)

          END IF

        END DO

        SELECT TYPE(self)
        TYPE IS(meshParticles)
          !Connect Vol-Edge
          DO et = 1, self%numEdges
            CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)

          END DO

        END SELECT

      END DO

    END SUBROUTINE connectMesh1DRad

    SUBROUTINE connectVolVol(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVol), INTENT(inout):: elemA
      CLASS(meshVol), INTENT(inout):: elemB

        SELECT TYPE(elemA)
        TYPE IS(meshVol1DRadSegm)
          SELECT TYPE(elemB)
          TYPE IS(meshVol1DRadSegm)
            CALL connectSegmSegm(elemA, elemB)

          END SELECT

        END SELECT

    END SUBROUTINE connectVolVol

    SUBROUTINE connectSegmSegm(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB

      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
          elemA%n2%n == elemB%n1%n) THEN
        elemA%e1 => elemB
        elemB%e2 => elemA
        
      END IF

      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
          elemA%n1%n == elemB%n2%n) THEN

        elemA%e2 => elemB
        elemB%e1 => elemA

      END IF

    END SUBROUTINE connectSegmSegm

    SUBROUTINE connectVolEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVol),  INTENT(inout):: elemA
      CLASS(meshEdge), INTENT(inout):: elemB

      SELECT TYPE(elemA)
      TYPE IS (meshVol1DRadSegm)
        SELECT TYPE(elemB)
        CLASS IS(meshEdge1DRad)
          CALL connectSegmEdge(elemA, elemB)

        END SELECT

      END SELECT

    END SUBROUTINE connectVolEdge

    SUBROUTINE connectSegmEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
      CLASS(meshEdge1DRad),    INTENT(inout), TARGET:: elemB

      IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
          elemA%n2%n == elemB%n1%n) THEN

        elemA%e1 => elemB
        elemB%e2 => elemA

        !Revers the normal to point inside the domain
        elemB%normal = - elemB%normal

      END IF

      IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
          elemA%n1%n == elemB%n1%n) THEN

        elemA%e2 => elemB
        elemB%e1 => elemA

      END IF

    END SUBROUTINE connectSegmEdge

END MODULE moduleMesh1DRad