!moduleMesh1D: 1D cartesian module
!              x == x
!              y == unused
!              z == unused
MODULE moduleMesh1D
  USE moduleMesh
  IMPLICIT NONE

  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1D
    !Element coordinates
    REAL(8):: x = 0.D0
    CONTAINS
      PROCEDURE, PASS:: init => initNode1D
      PROCEDURE, PASS:: getCoordinates => getCoord1D

  END TYPE meshNode1D

  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1D
    !Element coordinates
    REAL(8):: x = 0.D0
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL()
    CONTAINS
      PROCEDURE, PASS:: init     => initEdge1D
      PROCEDURE, PASS:: getNodes => getNodes1D
      PROCEDURE, PASS:: randPos  => randPos1D

  END TYPE meshEdge1D

  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1D
    CONTAINS
      PROCEDURE, PASS:: detJac => detJ1D
      PROCEDURE, PASS:: invJac => invJ1D
      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
      PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
      PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
  
  END TYPE meshVol1D

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

    END FUNCTION fPsi_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 meshVol1D
      CLASS(meshVol1D), 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(meshVol1D):: meshVol1DSegm
    !Element coordinates
    REAL(8):: x(1:2)
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
    !Connectivity to adjacent elements
    CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
    REAL(8):: arNodes(1:2)
    CONTAINS
      PROCEDURE, PASS::   init        => initVol1DSegm
      PROCEDURE, PASS::   area        => areaSegm
      PROCEDURE, NOPASS:: fPsi        => fPsiSegm
      PROCEDURE, NOPASS:: dPsi        => dPsiSegm
      PROCEDURE, PASS::   partialDer  => partialDerSegm
      PROCEDURE, PASS::   elemK       => elemKSegm
      PROCEDURE, PASS::   elemF       => elemFSegm
      PROCEDURE, NOPASS:: weight      => weightSegm
      PROCEDURE, NOPASS:: inside      => insideSegm
      PROCEDURE, PASS::   scatter     => scatterSegm
      PROCEDURE, PASS::   gatherEF    => gatherEFSegm
      PROCEDURE, PASS::   getNodes    => getNodesSegm
      PROCEDURE, PASS::   phy2log     => phy2logSegm
      PROCEDURE, PASS::   nextElement => nextElementSegm
      PROCEDURE, PASS::   resetOutput => resetOutputSegm

  END TYPE meshVol1DSegm

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

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

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

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

    END SUBROUTINE initNode1D

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

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

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

    END FUNCTION getCoord1D

    !EDGE FUNCTIONS
    !Inits edge element
    SUBROUTINE initEdge1D(self, n, p, bt, physicalSurface)
      IMPLICIT NONE

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

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

      self%x = r1(1)

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

      !Boundary index
      self%bt = bt
      !Physical Surface
      self%physicalSurface = physicalSurface

    END SUBROUTINE initEdge1D

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

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

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

    END FUNCTION getNodes1D

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

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

    END FUNCTION randPos1D

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

      CLASS(meshVol1DSegm), INTENT(out):: self
      INTEGER, INTENT(in):: n
      INTEGER, INTENT(in):: p(:)
      REAL(8), DIMENSION(1:3):: r1, r2

      self%n = n
      self%n1 => mesh%nodes(p(1))%obj
      self%n2 => mesh%nodes(p(2))%obj
      !Get element coordinates
      r1 = self%n1%getCoordinates()
      r2 = self%n2%getCoordinates()
      self%x = (/ 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)

      self%sigmaVrelMax = sigma_ref/L_ref**2

      CALL OMP_INIT_LOCK(self%lock)

    END SUBROUTINE initVol1DSegm

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

      CLASS(meshVol1DSegm), INTENT(inout):: self
      REAL(8):: l !element length
      REAL(8):: fPsi(1:2)
      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)
      l = 2.D0*detJ
      self%volume  =      l
      self%arNodes = fPsi*l

    END SUBROUTINE areaSegm

    !Computes element functions at point xii
    PURE FUNCTION fPsiSegm(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 fPsiSegm

    !Computes element derivative shape function at Xii
    PURE FUNCTION dPsiSegm(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 dPsiSegm

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

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

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

    END SUBROUTINE partialDerSegm

    !Computes local stiffness matrix
    PURE FUNCTION elemKSegm(self) RESULT(ke)
      IMPLICIT NONE

      CLASS(meshVol1DSegm), INTENT(in):: self
      REAL(8):: ke(1:2,1:2)
      REAL(8):: Xii(1:3)
      REAL(8):: dPsi(1:1, 1:2)
      REAL(8):: invJ

      ke      = 0.D0
      Xii     = (/ 0.D0, 0.D0, 0.D0 /)
      dPsi    = self%dPsi(Xii)
      invJ    = self%invJac(Xii, dPsi)
      ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) 
      ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) 
      ke      = 2.D0*ke*invJ

    END FUNCTION elemKSegm

    PURE FUNCTION elemFSegm(self, source) RESULT(localF)
      IMPLICIT NONE

      CLASS(meshVol1DSegm), INTENT(in):: self
      REAL(8), INTENT(in):: source(1:)
      REAL(8), ALLOCATABLE:: localF(:)
      REAL(8):: fPsi(1:2)
      REAL(8):: detJ
      REAL(8):: Xii(1:3)

      Xii = 0.D0
      fPsi = self%fPsi(Xii)
      detJ = self%detJac(Xii)
      ALLOCATE(localF(1:2))
      localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ

    END FUNCTION elemFSegm

    PURE FUNCTION weightSegm(xi) RESULT(w)
      IMPLICIT NONE

      REAL(8), INTENT(in):: xi(1:3)
      REAL(8):: w(1:3)

      w = fPsiSegm(xi)

    END FUNCTION weightSegm

    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

    SUBROUTINE scatterSegm(self, part)
      USE moduleOutput
      USE moduleSpecies
      IMPLICIT NONE

      CLASS(meshVol1DSegm), INTENT(in):: self
      CLASS(particle), INTENT(in):: part
      TYPE(outputNode), POINTER:: vertex
      REAL(8):: w_p(1:2)
      REAL(8):: tensorS(1:3,1:3)

      w_p = self%weight(part%xi)
      tensorS = outerProduct(part%v, part%v)

      vertex => self%n1%output(part%sp)
      vertex%den          = vertex%den          + part%weight*w_p(1)
      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(1)*part%v(:)
      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS

      vertex => self%n2%output(part%sp)
      vertex%den          = vertex%den          + part%weight*w_p(2)
      vertex%mom(:)       = vertex%mom(:)       + part%weight*w_p(2)*part%v(:)
      vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS

    END SUBROUTINE scatterSegm

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

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

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

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

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

    END FUNCTION getNodesSegm

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

      CLASS(meshVol1DSegm), 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%x(1))/(self%x(2) - self%x(1)) - 1.D0

    END FUNCTION phy2logSegm

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

      CLASS(meshVol1DSegm), INTENT(in):: self
      REAL(8), INTENT(in):: xi(1:3)
      CLASS(*), 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 nextElementSegm

    !Reset the output of nodes in element
    PURE SUBROUTINE resetOutputSegm(self)
      USE moduleSpecies
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshVol1DSegm), INTENT(inout):: self
      INTEGER:: k

      DO k = 1, nSpecies
        self%n1%output(k)%den     = 0.D0
        self%n1%output(k)%mom     = 0.D0
        self%n1%output(k)%tensorS = 0.D0

        self%n2%output(k)%den     = 0.D0
        self%n2%output(k)%mom     = 0.D0
        self%n2%output(k)%tensorS = 0.D0

      END DO

    END SUBROUTINE resetOutputSegm


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

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

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

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


END MODULE moduleMesh1D