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

  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
    !Element coordinates
    REAL(8):: r = 0.D0
    CONTAINS
      !meshNode DEFERRED PROCEDURES
      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
      !meshEdge DEFERRED PROCEDURES
      PROCEDURE, PASS::   init         => initEdge1DRad
      PROCEDURE, PASS::   getNodes     => getNodes1DRad
      PROCEDURE, PASS::   intersection => intersection1DRad
      PROCEDURE, PASS::   randPos      => randPos1DRad
      procedure, nopass:: fPsi         => fPsiPoint

  END TYPE meshEdge1DRad

  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DRadSegm
    !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()
    CONTAINS
      !meshCell DEFERRED PROCEDURES
      PROCEDURE, PASS::   init                => initCell1DRadSegm
      PROCEDURE, PASS::   getNodes            => getNodesSegm
      PROCEDURE, PASS::   randPos             => randPos1DRadSegm
      PROCEDURE, NOPASS:: fPsi                => fPsiSegm
      PROCEDURE, NOPASS:: dPsi                => dPsiSegm
      PROCEDURE, PASS::   partialDer          => partialDerSegm
      PROCEDURE, NOPASS:: detJac              => detJ1DRad
      PROCEDURE, NOPASS:: invJac              => invJ1DRad
      PROCEDURE, PASS::   gatherElectricField => gatherEFSegm
      PROCEDURE, PASS::   gatherMagneticField => gatherMFSegm
      PROCEDURE, PASS::   elemK               => elemKSegm
      PROCEDURE, PASS::   elemF               => elemFSegm
      PROCEDURE, NOPASS:: inside              => insideSegm
      PROCEDURE, PASS::   phy2log             => phy2logSegm
      PROCEDURE, PASS::   neighbourElement    => neighbourElementSegm
      !PARTICLUAR PROCEDURES
      PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm

  END TYPE meshCell1DRadSegm

  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

      !Allocate 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
    !Init edge element
    SUBROUTINE initEdge1DRad(self, n, p, ps)
      USE moduleSpecies, only:nSpecies
      USE moduleErrors
      IMPLICIT NONE

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

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

      self%r = r1(1)

      self%surface = 1.D0

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

      ! Get Physical Surface index based on input numering
      ps_index = physicalSurface_to_index(ps)

      ! Add elements to physical surface
      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)

      ! Link edge to particle boundaries
      allocate(self%boundariesParticle(1:nSpecies))
      do s = 1, nSpecies
        self%boundariesParticle(s)%obj => physicalSurfaces(ps_index)%particles(s)%obj

      end do

    end subroutine initEdge1DRad

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

      CLASS(meshEdge1DRad), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      INTEGER:: n(1:nNodes)

      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

    !Calculate 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

    !CELL FUNCTIONS
    !SEGMENT FUNCTIONS
    !Init element
    SUBROUTINE initCell1DRadSegm(self, n, p, nodes)
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), 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%nNodes = SIZE(p)
      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%calculateVolume()

      CALL OMP_INIT_LOCK(self%lock)

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

    END SUBROUTINE initCell1DRadSegm

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

      CLASS(meshCell1DRadSegm), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      INTEGER:: n(1:nNodes)

      n = (/ self%n1%n, self%n2%n /)

    END FUNCTION getNodesSegm

    !Random position in 1D volume
    FUNCTION randPos1DRadSegm(self) RESULT(r)
      USE moduleRandom
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(in):: self
      REAL(8):: r(1:3)
      REAL(8):: Xi(1:3)
      REAL(8):: fPsi(1:2)

      Xi    = 0.D0
      Xi(1) = random(-1.D0, 1.D0)

      fPsi = self%fPsi(Xi, 2)

      r    = 0.D0
      r(1) = DOT_PRODUCT(fPsi, self%r)

    END FUNCTION randPos1DRadSegm

    !Partial derivative in global coordinates
    PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
      REAL(8):: pDer(1:3, 1:3)

      pDer = 0.D0

      pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%r(1:2))
      pDer(2,2) = 1.D0
      pDer(3,3) = 1.D0

    END FUNCTION partialDerSegm

    PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: phi(1:2)

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

      array = -self%gatherDF(Xi, 2, phi)

    END FUNCTION gatherEFSegm

    PURE FUNCTION gatherMFSegm(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell1DRadSegm), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: B(1:2,1:3)

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

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

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

      array = self%gatherF(Xi, 2, B)

    END FUNCTION gatherMFSegm

    !Compute element local stiffness matrix
    PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK)
      USE moduleConstParam, ONLY: PI2
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8):: localK(1:nNodes,1:nNodes)
      REAL(8):: Xi(1:3)
      REAL(8):: dPsi(1:3, 1:2)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: r
      REAL(8):: invJ(1:3, 1:3), detJ
      INTEGER:: l

      localK = 0.D0

      Xi = 0.D0
      !Start 1D Gauss Quad Integral
      DO l = 1, 3
        Xi(1)  = corSeg(l)
        dPsi   = self%dPsi(Xi, 2)
        pDer   = self%partialDer(2, dPsi)
        detJ   = self%detJac(pDer)
        invJ   = self%invJac(pDer)
        r      = self%gatherF(Xi, 2, self%r)
        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
                                           MATMUL(invJ,dPsi))* &
                          r*wSeg(l)/detJ

      END DO
      localK = localK*PI2

    END FUNCTION elemKSegm

    !Compute the local source vector for a force f
    PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF)
      USE moduleConstParam, ONLY: PI2
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: source(1:nNodes)
      REAL(8):: localF(1:nNodes)
      REAL(8):: fPsi(1:2)
      REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
      REAL(8):: Xi(1:3)
      REAL(8):: detJ, f
      REAL(8):: r
      INTEGER:: l

      localF = 0.D0

      Xi = 0.D0
      !Start 1D Gauss Quad Integral
      DO l = 1, 3
        Xi(1) = corSeg(l)
        dPsi   = self%dPsi(Xi, 2)
        pDer   = self%partialDer(2, dPsi)
        detJ   = self%detJac(pDer)
        fPsi   = self%fPsi(Xi, 2)
        r      = DOT_PRODUCT(fPsi, self%r)
        f      = DOT_PRODUCT(fPsi, source)
        localF = localF + r*f*fPsi*wSeg(l)*detJ

      END DO
      localF = localF*PI2

    END FUNCTION elemFSegm

    PURE FUNCTION insideSegm(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 insideSegm

    PURE FUNCTION phy2logSegm(self, r) RESULT(Xi)
      IMPLICIT NONE

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

      Xi = 0.D0

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

    END FUNCTION phy2logSegm

    !Get the next element for a logical position Xi
    SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement)
      IMPLICIT NONE

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

      ELSEIF (Xi(1) >  1.D0) THEN
        neighbourElement => self%e1

      END IF

    END SUBROUTINE neighbourElementSegm

    !Compute element volume
    PURE SUBROUTINE volumeSegm(self)
      USE moduleConstParam, ONLY: PI4
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(inout):: self
      REAL(8):: Xi(1:3)
      REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
      REAL(8):: detJ
      REAL(8):: fPsi(1:2)
      REAL(8):: r

      self%volume  = 0.D0
      !1D 1 point Gauss Quad Integral
      Xi = 0.D0
      dPsi = self%dPsi(Xi, 2)
      pDer = self%partialDer(2, dPsi)
      detJ = self%detJac(pDer)
      fPsi = self%fPsi(Xi, 2)
      r = DOT_PRODUCT(fPsi, self%r)
      !Compute total volume of the cell
      self%volume  = r*detJ*PI4 !2*2PI
      !Compute volume per node
      Xi = (/-5.D-1,  0.D0, 0.D0/)
      r = self%gatherF(Xi, 2, self%r)
      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI4
      Xi = (/ 5.D-1,  0.D0, 0.D0/)
      r = self%gatherF(Xi, 2, self%r)
      self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI4

    END SUBROUTINE volumeSegm

    !COMMON FUNCTIONS FOR 1D CELL ELEMENTS
    !Compute element Jacobian determinant
    PURE FUNCTION detJ1DRad(pDer) RESULT(dJ)
      IMPLICIT NONE

      REAL(8), INTENT(in):: pDer(1:3, 1:3)
      REAL(8):: dJ

      dJ = pDer(1, 1)

    END FUNCTION detJ1DRad

    !Compute element Jacobian inverse matrix (without determinant)
    PURE FUNCTION invJ1DRad(pDer) RESULT(invJ)
      IMPLICIT NONE

      REAL(8), INTENT(in):: pDer(1:3, 1:3)
      REAL(8):: invJ(1:3,1:3)

      invJ = 0.D0

      invJ(1, 1) = 1.D0/pDer(1, 1)
      invJ(2, 2) = 1.D0
      invJ(3, 3) = 1.D0

    END FUNCTION invJ1DRad

    SUBROUTINE connectMesh1DRad(self)
      IMPLICIT NONE

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

      DO e = 1, self%numCells
        !Connect Cell-Cell
        DO et = 1, self%numCells
          IF (e /= et) THEN
            CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)

          END IF

        END DO

        SELECT TYPE(self)
        TYPE IS(meshParticles)
          !Connect Cell-Edge
          DO et = 1, self%numEdges
            CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)

          END DO

        END SELECT

      END DO

    END SUBROUTINE connectMesh1DRad

    SUBROUTINE connectCellCell(elemA, elemB)
      IMPLICIT NONE

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

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

          END SELECT

        END SELECT

    END SUBROUTINE connectCellCell

    SUBROUTINE connectSegmSegm(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemA
      CLASS(meshCell1DRadSegm), 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 connectCellEdge(elemA, elemB)
      IMPLICIT NONE

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

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

        END SELECT

      END SELECT

    END SUBROUTINE connectCellEdge

    SUBROUTINE connectSegmEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshCell1DRadSegm), 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

        !Rever 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