!moduleMesh2DCyl: 2D aXial symmetric extension of generic mesh from GMSH format.
!               x == z
!               y == r
!               z == theta (unused)
MODULE moduleMesh2DCyl
  USE moduleMesh
  use moduleMeshCommon
  IMPLICIT NONE

  TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl
    !Element coordinates
    REAL(8):: r = 0.D0, z = 0.D0
    CONTAINS
      !meshNode DEFERRED PROCEDURES
      PROCEDURE, PASS:: init           => initNode2DCyl
      PROCEDURE, PASS:: getCoordinates => getCoord2DCyl

  END TYPE meshNode2DCyl

  TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCyl
    !Element coordinates
    REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
    CONTAINS
      !meshEdge DEFERRED PROCEDURES
      PROCEDURE, PASS::   init         => initEdge2DCyl
      PROCEDURE, PASS::   getNodes     => getNodes2DCyl
      PROCEDURE, PASS::   intersection => intersection2DCylEdge
      PROCEDURE, PASS::   randPos      => randPosEdge
      procedure, pass::   center       => centerEdgeSegm
      procedure, nopass:: centerXi     => centerXiSegm
      procedure, nopass:: fPsi         => fPsiSegm

  END TYPE meshEdge2DCyl

  !Quadrilateral volume element
  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylQuad
    !Element coordinates
    REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
    !Connectivity to adjacent elements
    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()

    CONTAINS
      !meshCell DEFERRED PROCEDURES
      PROCEDURE, PASS::   init                => initCellQuad2DCyl
      PROCEDURE, PASS::   getNodes            => getNodesQuad
      PROCEDURE, PASS::   randPos             => randPosCellQuad
      PROCEDURE, NOPASS:: fPsi                => fPsiQuad
      PROCEDURE, NOPASS:: dPsi                => dPsiQuad
      PROCEDURE, PASS::   partialDer          => partialDerQuad
      PROCEDURE, NOPASS:: detJac              => detJ2DCyl
      PROCEDURE, NOPASS:: invJac              => invJ2DCyl
      PROCEDURE, PASS::   gatherElectricField => gatherEFQuad
      PROCEDURE, PASS::   gatherMagneticField => gatherMFQuad
      PROCEDURE, PASS::   elemK               => elemKQuad
      PROCEDURE, PASS::   elemF               => elemFQuad
      PROCEDURE, NOPASS:: inside              => insideQuad
      PROCEDURE, PASS::   phy2log             => phy2logQuad
      PROCEDURE, PASS::   neighbourElement    => neighbourElementQuad
      !PARTICLUAR PROCEDURES
      PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad

  END TYPE meshCell2DCylQuad

  !Triangular volume element
  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylTria
    !Element coordinates
    REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0
    !Connectivity to nodes
    CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
    !Connectivity to adjacent elements
    CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()

    CONTAINS
      !meshCell DEFERRED PROCEDURES
      PROCEDURE, PASS::   init                => initCellTria2DCyl
      PROCEDURE, PASS::   getNodes            => getNodesTria
      PROCEDURE, PASS::   randPos             => randPosCellTria
      PROCEDURE, NOPASS:: fPsi                => fPsiTria
      PROCEDURE, NOPASS:: dPsi                => dPsiTria
      PROCEDURE, PASS::   partialDer          => partialDerTria
      PROCEDURE, NOPASS:: detJac              => detJ2DCyl
      PROCEDURE, NOPASS:: invJac              => invJ2DCyl
      PROCEDURE, PASS::   gatherElectricField => gatherEFTria
      PROCEDURE, PASS::   gatherMagneticField => gatherMFTria
      PROCEDURE, PASS::   elemK               => elemKTria
      PROCEDURE, PASS::   elemF               => elemFTria
      PROCEDURE, NOPASS:: inside              => insideTria
      PROCEDURE, PASS::   phy2log             => phy2logTria
      PROCEDURE, PASS::   neighbourElement    => neighbourElementTria
      !PARTICULAR PROCEDURES
      PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria

  END TYPE meshCell2DCylTria

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

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

      self%n = n
      self%z = r(1)/L_ref
      self%r = r(2)/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 initNode2DCyl

    !Get coordinates from node
    PURE FUNCTION getCoord2DCyl(self) RESULT(r)
      IMPLICIT NONE

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

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

    END FUNCTION getCoord2DCyl

    !EDGE FUNCTIONS
    !Init edge element
    SUBROUTINE initEdge2DCyl(self, n, p, ps)
      USE moduleSpecies
      USE moduleErrors
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE

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

      self%n  = n
      self%nNodes = SIZE(p)
      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%z  = (/r1(1), r2(1)/)
      self%r  = (/r1(2), r2(2)/)
      !Edge surface
      IF (self%z(2) /= self%z(1)) THEN
        self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1))

      ELSE
        self%surface = ABS(self%r(2)**2 - self%r(1)**2)

      END IF
      self%surface = self%surface * PI
      !Normal vector
      self%normal = (/ -(self%r(2)-self%r(1)), &
                         self%z(2)-self%z(1) , &
                         0.D0 /)
      self%normal = self%normal/NORM2(self%normal)

      ! 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)
      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n2%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 initEdge2DCyl

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

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

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

    END FUNCTION getNodes2DCyl

    !Calculate intersection between position and edge
    PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r)
      IMPLICIT NONE

      CLASS(meshEdge2DCyl), INTENT(in):: self
      REAL(8), DIMENSION(1:3), INTENT(in):: r0
      REAL(8), DIMENSION(1:3):: r
      REAL(8), DIMENSION(1:3):: edge0, edgeV
      REAL(8):: tI

      edge0 = (/self%z(1), self%r(1), 0.D0 /)
      edgeV = (/self%z(2), self%r(2), 0.D0 /) - edge0

      tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)

      r = edge0 + tI*edgeV

    END FUNCTION intersection2DCylEdge

    !Calculate a random position in edge
    FUNCTION randPosEdge(self) RESULT(r)
      USE moduleRandom
      IMPLICIT NONE

      CLASS(meshEdge2DCyl), INTENT(in):: self
      real(8):: Xi(1:3)
      real(8):: r(1:3)
      real(8):: fPsi(1:2)

      Xi    = 0.d0
      Xi(1) = random()

      fPsi = self%fPsi(Xi, 2)

      r = (/dot_product(fPsi, self%z), &
            dot_product(fPsi, self%r), &
            0.d0/)

    END FUNCTION randPosEdge

    function centerEdgeSegm(self) result(r)
      use moduleMeshCommon, only: cenSeg
      implicit none

      class(meshEdge2DCyl), intent(in):: self
      real(8):: r(1:3)
      real(8):: fPsi(1:2)

      fPsi = self%fPsi(cenSeg, 2)

      r = (/dot_product(fPsi, self%z), &
            dot_product(fPsi, self%r), &
            0.d0/)

    end function centerEdgeSegm

    !CELL FUNCTIONS
    !QUAD FUNCTIONS
    !Init element
    SUBROUTINE initCellQuad2DCyl(self, n, p, nodes)
      IMPLICIT NONE

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

      !Assign node index
      self%n  = n

      !Assign number of nodes of cell
      self%nNodes = SIZE(p)

      !Assign nodes to element
      self%n1 => nodes(p(1))%obj
      self%n2 => nodes(p(2))%obj
      self%n3 => nodes(p(3))%obj
      self%n4 => nodes(p(4))%obj
      !Get element coordinates
      r1 = self%n1%getCoordinates()
      r2 = self%n2%getCoordinates()
      r3 = self%n3%getCoordinates()
      r4 = self%n4%getCoordinates()
      self%z  = (/r1(1), r2(1), r3(1), r4(1)/)
      self%r  = (/r1(2), r2(2), r3(2), r4(2)/)

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

    !Get nodes from quadrilateral element
    PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n)
      IMPLICIT NONE

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

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

    END FUNCTION getNodesQuad

    !Random position in quadrilateral volume
    FUNCTION randPosCellQuad(self) RESULT(r)
      USE moduleRandom
      IMPLICIT NONE

      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8):: r(1:3)
      REAL(8):: Xi(1:3)
      REAL(8):: fPsi(1:4)

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

      fPsi = self%fPsi(Xi, 4)

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

    END FUNCTION randPosCellQuad

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

      CLASS(meshCell2DCylQuad), 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:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), &
                        DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /)
      pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), &
                        DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /)
      pDer(3,3) = 1.D0

    END FUNCTION partialDerQuad

    PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: phi(1:4)

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

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

    END FUNCTION gatherEFQuad

    PURE FUNCTION gatherMFQuad(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: B(1:4,1:3)

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

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

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

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

    END FUNCTION gatherMFQuad

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

      CLASS(meshCell2DCylQuad), 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:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: r
      REAL(8):: invJ(1:3,1:3), detJ
      INTEGER:: l, m

      localK = 0.D0

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

        END DO
      END DO
      localK = localK*PI2

    END FUNCTION elemKQuad

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

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

      localF = 0.D0

      Xi = 0.D0
      DO l = 1, 3
        Xi(1) = corQuad(l)
        DO m = 1, 3
          Xi(2)  = corQuad(m)
          dPsi   = self%dPsi(Xi, 4)
          pDer   = self%partialDer(4, dPsi)
          detJ   = self%detJac(pDer)
          fPsi = self%fPsi(Xi, 4)
          r      = DOT_PRODUCT(fPsi, self%r)
          f      = DOT_PRODUCT(fPsi, source)
          localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ

         END DO
      END DO
      localF = localF*PI2

    END FUNCTION elemFQuad

    !Check if Xi is inside the element
    PURE FUNCTION insideQuad(Xi) RESULT(ins)
      IMPLICIT NONE

      REAL(8), INTENT(in):: Xi(1:3)
      LOGICAL:: ins
      real(8):: eps

      eps = 1.0d-10

      ins = (Xi(1) >= -1.D0-eps .AND. Xi(1) <= 1.D0+eps) .AND. &
            (Xi(2) >= -1.D0-eps .AND. Xi(2) <= 1.D0+eps)

    END FUNCTION insideQuad

    !Transform physical coordinates to element coordinates with a Taylor series
    PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
      IMPLICIT NONE

      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
      REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: conv

      !Iterative newton method to transform coordinates
      conv  = 1.D0
      Xi0   = 0.D0
      Xi(3) = 0.D0

      DO WHILE(conv > 1.D-4)
        fPsi           = self%fPsi(Xi0, 4)
        x0(1)          = dot_product(fPsi, self%z)
        x0(2)          = dot_product(fPsi, self%r)
        deltaR         = r(1:2) - x0
        dPsi           = self%dPsi(Xi0, 4)
        pDer           = self%partialDer(4, dPsi)
        detJ           = self%detJac(pDer)
        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
        Xi0(1:2)       = Xi(1:2)

      END DO

    END FUNCTION phy2logQuad

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

      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
      REAL(8):: XiArray(1:4)
      INTEGER:: nextInt

      XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
      nextInt = MAXLOC(XiArray,1)
      !Select the higher value of directions and searches in that direction
      NULLIFY(neighbourElement)
      SELECT CASE (nextInt)
      CASE (1)
        neighbourElement => self%e1
      CASE (2)
        neighbourElement => self%e2
      CASE (3)
        neighbourElement => self%e3
      CASE (4)
        neighbourElement => self%e4
      END SELECT

    END SUBROUTINE neighbourElementQuad

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

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

      self%volume  = 0.D0

      !2D 1 point Gauss Quad Integral
      Xi = 0.D0
      dPsi = self%dPsi(Xi, 4)
      pDer = self%partialDer(4, dPsi)
      detJ = self%detJac(pDer)
      fPsi = self%fPsi(Xi, 4)
      r = DOT_PRODUCT(fPsi,self%r)
      !Computes total volume of the cell
      self%volume  = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral)
      !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis.
      Xi = (/-5.D-1, -0.25D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
      Xi = (/ 5.D-1, -0.25D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
      Xi = (/ 5.D-1,  0.75D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
      Xi = (/-5.D-1,  0.75D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8

    END SUBROUTINE volumeQuad

    !TRIA FUNCTIONS
    !Init element
    SUBROUTINE initCellTria2DCyl(self, n, p, nodes)
      USE moduleRefParam
      IMPLICIT NONE

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

      !Assign node index
      self%n = n

      !Assign number of nodes of cell
      self%nNodes = SIZE(p)

      !Assign nodes to element
      self%n1 => nodes(p(1))%obj
      self%n2 => nodes(p(2))%obj
      self%n3 => nodes(p(3))%obj
      !Get element coordinates
      r1 = self%n1%getCoordinates()
      r2 = self%n2%getCoordinates()
      r3 = self%n3%getCoordinates()
      self%z  = (/r1(1), r2(1), r3(1)/)
      self%r  = (/r1(2), r2(2), r3(2)/)
      !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 initCellTria2DCyl

    !Random position in cell
    PURE FUNCTION getNodesTria(self, nNodes) RESULT(n)
      IMPLICIT NONE

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

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

    END FUNCTION getNodesTria

    !Random position in cell
    FUNCTION randPosCellTria(self) RESULT(r)
      USE moduleRandom
      IMPLICIT NONE

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

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

      fPsi = self%fPsi(Xi, 4)

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

    END FUNCTION randPosCellTria

    !Compute the derivatives in global coordinates
    PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer)
      IMPLICIT NONE

      CLASS(meshCell2DCylTria), 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:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%z(1:3)), &
                        DOT_PRODUCT(dPsi(2,1:3),self%z(1:3)) /)
      pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%r(1:3)), &
                        DOT_PRODUCT(dPsi(2,1:3),self%r(1:3)) /)

    END FUNCTION partialDerTria

    !Gather electric field at position Xi
    PURE FUNCTION gatherEFTria(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell2DCylTria), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: phi(1:3)

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

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

    END FUNCTION gatherEFTria

    !Gather magnetic field at position Xi
    PURE FUNCTION gatherMFTria(self, Xi) RESULT(array)
      IMPLICIT NONE
      CLASS(meshCell2DCylTria), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)
      REAL(8):: B(1:3,1:3)

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

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

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

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

    END FUNCTION gatherMFTria

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

      CLASS(meshCell2DCylTria), 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:3)
      REAL(8):: r
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: invJ(1:3,1:3), detJ
      INTEGER:: l

      localK = 0.D0

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

      END DO
      localK = localK*PI2

    END FUNCTION elemKTria

    !Compute element local source vector
    PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF)
      USE moduleConstParam, ONLY: PI2
      IMPLICIT NONE

      CLASS(meshCell2DCylTria), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: source(1:nNodes)
      REAL(8):: localF(1:nNodes)
      REAL(8):: fPsi(1:3)
      REAL(8):: dPsi(1:3, 1:3), 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 2D Gauss Quad Integral
      DO l = 1, 4
        Xi(1)  = Xi1Tria(l)
        Xi(2)  = Xi2Tria(l)
        dPsi   = self%dPsi(Xi, 3)
        pDer   = self%partialDer(3, dPsi)
        detJ   = self%detJac(pDer)
        fPsi   = self%fPsi(Xi, 3)
        r      = DOT_PRODUCT(fPsi, self%r)
        f      = DOT_PRODUCT(fPsi, source)
        localF = localF + r*f*fPsi*wTria(l)*detJ

      END DO
      localF = localF*PI2

    END FUNCTION elemFTria

    !Check if Xi is inside the element
    PURE FUNCTION insideTria(Xi) RESULT(ins)
      IMPLICIT NONE

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

      ins =        Xi(1)         >= 0.D0 .AND. &
                           Xi(2) >= 0.D0 .AND. &
            1.D0 - Xi(1) - Xi(2) >= 0.D0

    END FUNCTION insideTria

    !Transform physical coordinates to element coordinates
    PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
      IMPLICIT NONE

      CLASS(meshCell2DCylTria), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
      REAL(8):: dPsi(1:3,1:4)
      REAL(8):: pDer(1:3, 1:3)

      !Direct method to convert coordinates
      Xi(3)          = 0.D0
      deltaR         = (/ r(1) - self%z(1), r(2) - self%r(1) /)
      dPsi           = self%dPsi(Xi, 3)
      pDer           = self%partialDer(3, dPsi)
      detJ           = self%detJac(pDer)
      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ

    END FUNCTION phy2logTria

    !Get the neighbour cell for a logical position Xi
    SUBROUTINE neighbourElementTria(self, Xi, neighbourElement)
      IMPLICIT NONE

      CLASS(meshCell2DCylTria), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
      REAL(8):: XiArray(1:3)
      INTEGER:: nextInt

      XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /)
      nextInt = MINLOC(XiArray,1)
      NULLIFY(neighbourElement)
      SELECT CASE (nextInt)
      CASE (1)
        neighbourElement => self%e1
      CASE (2)
        neighbourElement => self%e2
      CASE (3)
        neighbourElement => self%e3
      END SELECT

    END SUBROUTINE neighbourElementTria

    !Calculate volume for triangular element
    PURE SUBROUTINE volumeTria(self)
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE

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

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

    END SUBROUTINE volumeTria

    !COMMON FUNCTIONS FOR CYLINDRICAL CELL ELEMENTS
    !Compute element Jacobian determinant
    PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ)
      IMPLICIT NONE

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

      dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1)

    END FUNCTION detJ2DCyl

    !Compute element Jacobian inverse matrix (without determinant)
    PURE FUNCTION invJ2DCyl(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:2) = (/  pDer(2,2), -pDer(2,1) /)
      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
      invJ(3, 3)   = 1.D0

    END FUNCTION invJ2DCyl

    SUBROUTINE connectMesh2DCyl(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 connectMesh2DCyl

    !Selects type of elements to build connection
    SUBROUTINE connectCellCell(elemA, elemB)
      IMPLICIT NONE

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

      SELECT TYPE(elemA)
      TYPE IS(meshCell2DCylQuad)
        !Element A is a quadrilateral
        SELECT TYPE(elemB)
        TYPE IS(meshCell2DCylQuad)
          !Element B is a quadrilateral
          CALL connectQuadQuad(elemA, elemB)

        TYPE IS(meshCell2DCylTria)
          !Element B is a triangle
          CALL connectQuadTria(elemA, elemB)

        END SELECT

      TYPE IS(meshCell2DCylTria)
        !Element A is a Triangle
        SELECT TYPE(elemB)
        TYPE IS(meshCell2DCylQuad)
          !Element B is a quadrilateral
          CALL connectQuadTria(elemB, elemA)

        TYPE IS(meshCell2DCylTria)
          !Element B is a triangle
          CALL connectTriaTria(elemA, elemB)

        END SELECT

      END SELECT

    END SUBROUTINE connectCellCell

    SUBROUTINE connectCellEdge(elemA, elemB)
      IMPLICIT NONE

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

      SELECT TYPE(elemB)
      CLASS IS(meshEdge2DCyl)
        SELECT TYPE(elemA)
        TYPE IS(meshCell2DCylQuad)
          !Element A is a quadrilateral
          CALL connectQuadEdge(elemA, elemB)

        TYPE IS(meshCell2DCylTria)
          !Element A is a triangle
          CALL connectTriaEdge(elemA, elemB)

        END SELECT

      END SELECT

    END SUBROUTINE connectCellEdge

    SUBROUTINE connectQuadQuad(elemA, elemB)
      IMPLICIT NONE

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

      !Check direction 1
      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
        IF     (elemA%n1%n == elemB%n4%n .AND. &
                elemA%n2%n == elemB%n3%n) THEN
          elemA%e1 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
                elemA%n2%n == elemB%n2%n) THEN
          elemA%e1 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                elemA%n2%n == elemB%n1%n) THEN
          elemA%e1 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n1%n == elemB%n1%n .AND. &
                elemA%n2%n == elemB%n4%n) THEN
          elemA%e1 => elemB
          elemB%e4 => elemA

        END IF

      END IF

      !Check direction 2
      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
        IF     (elemA%n2%n == elemB%n1%n .AND. &
                elemA%n3%n == elemB%n4%n) THEN
          elemA%e2 => elemB
          elemB%e4 => elemA

        ELSEIF (elemA%n2%n == elemB%n4%n .AND. &
                elemA%n3%n == elemB%n3%n) THEN
          elemA%e2 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
                elemA%n3%n == elemB%n2%n) THEN
          elemA%e2 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                elemA%n3%n == elemB%n1%n) THEN
          elemA%e2 => elemB
          elemB%e1 => elemA

        END IF

      END IF

      !Check direction 3
      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
        IF     (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n4%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n3%n == elemB%n1%n .AND. &
                elemA%n4%n == elemB%n4%n) THEN
          elemA%e3 => elemB
          elemB%e4 => elemA

        ELSEIF (elemA%n3%n == elemB%n4%n .AND. &
                elemA%n4%n == elemB%n3%n) THEN
          elemA%e3 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
                elemA%n4%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

        END IF

      END IF

      !Check direction 4
      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
        IF     (elemA%n4%n == elemB%n3%n .AND. &
                elemA%n1%n == elemB%n2%n) THEN
          elemA%e4 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e4 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n4%n == elemB%n1%n .AND. &
                elemA%n1%n == elemB%n4%n) THEN
          elemA%e4 => elemB
          elemB%e4 => elemA

        ELSEIF (elemA%n4%n == elemB%n4%n .AND. &
                elemA%n1%n == elemB%n3%n) THEN
          elemA%e4 => elemB
          elemB%e3 => elemA

        END IF

      END IF

    END SUBROUTINE connectQuadQuad

    SUBROUTINE connectQuadTria(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA
      CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemB

      !Check direction 1
      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
        IF     (elemA%n1%n == elemB%n1%n .AND. &
                elemA%n2%n == elemB%n3%n) THEN
          elemA%e1 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
                elemA%n2%n == elemB%n2%n) THEN
          elemA%e1 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                elemA%n2%n == elemB%n1%n) THEN
          elemA%e1 => elemB
          elemB%e1 => elemA

        END IF

      END IF

      !Check direction 2
      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
        IF     (elemA%n2%n == elemB%n1%n .AND. &
                elemA%n3%n == elemB%n3%n) THEN
          elemA%e2 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
                elemA%n3%n == elemB%n2%n) THEN
          elemA%e2 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                elemA%n3%n == elemB%n1%n) THEN
          elemA%e2 => elemB
          elemB%e1 => elemA

        END IF

      END IF

      !Check direction 3
      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
        IF     (elemA%n3%n == elemB%n1%n .AND. &
                elemA%n4%n == elemB%n3%n) THEN
          elemA%e3 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
                elemA%n4%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n4%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        END IF

      END IF

      !Check direction 4
      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
        IF     (elemA%n4%n == elemB%n1%n .AND. &
                elemA%n1%n == elemB%n3%n) THEN
          elemA%e4 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
                elemA%n1%n == elemB%n2%n) THEN
          elemA%e4 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e4 => elemB
          elemB%e1 => elemA

        END IF

      END IF

    END SUBROUTINE connectQuadTria

    SUBROUTINE connectTriaTria(elemA, elemB)
      IMPLICIT NONE

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

      !Check direction 1
      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
        IF     (elemA%n1%n == elemB%n1%n .AND. &
                elemA%n2%n == elemB%n3%n) THEN
          elemA%e1 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                elemA%n2%n == elemB%n1%n) THEN
          elemA%e1 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
                elemA%n2%n == elemB%n2%n) THEN
          elemA%e1 => elemB
          elemB%e2 => elemA

        END IF

      END IF

      !Check direction 2
      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
        IF     (elemA%n2%n == elemB%n1%n .AND. &
                elemA%n3%n == elemB%n3%n) THEN
          elemA%e2 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                elemA%n3%n == elemB%n1%n) THEN
          elemA%e2 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
                elemA%n3%n == elemB%n2%n) THEN
          elemA%e2 => elemB
          elemB%e2 => elemA

        END IF


      END IF

      !Check direction 3
      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
        IF     (elemA%n3%n == elemB%n1%n .AND. &
                elemA%n1%n == elemB%n3%n) THEN
          elemA%e3 => elemB
          elemB%e3 => elemA

        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
                elemA%n1%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

        END IF


      END IF

    END SUBROUTINE connectTriaTria

    SUBROUTINE connectQuadEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA
      CLASS(meshEdge2DCyl),    INTENT(inout), TARGET:: elemB

      !Check direction 1
      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
        IF     (elemA%n1%n == elemB%n1%n .AND. &
                elemA%n2%n == elemB%n2%n) THEN
          elemA%e1 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n1%n == elemB%n2%n .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

      END IF

      !Check direction 2
      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
        IF     (elemA%n2%n == elemB%n1%n .AND. &
                elemA%n3%n == elemB%n2%n) THEN
          elemA%e2 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                elemA%n3%n == elemB%n1%n) THEN
          elemA%e2 => elemB
          elemB%e2 => elemA

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

        END IF

      END IF

      !Check direction 3
      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
        IF     (elemA%n3%n == elemB%n1%n .AND. &
                elemA%n4%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n4%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

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

        END IF

      END IF

      !Check direction 4
      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
        IF     (elemA%n4%n == elemB%n1%n .AND. &
                elemA%n1%n == elemB%n2%n) THEN
          elemA%e4 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e4 => elemB
          elemB%e2 => elemA

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

        END IF

      END IF

    END SUBROUTINE connectQuadEdge

    SUBROUTINE connectTriaEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemA
      CLASS(meshEdge2DCyl),    INTENT(inout), TARGET:: elemB

      !Check direction 1
      IF (.NOT. ASSOCIATED(elemA%e1)) THEN
        IF     (elemA%n1%n == elemB%n1%n .AND. &
                elemA%n2%n == elemB%n2%n) THEN
          elemA%e1 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n1%n == elemB%n2%n .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

      END IF

      !Check direction 2
      IF (.NOT. ASSOCIATED(elemA%e2)) THEN
        IF     (elemA%n2%n == elemB%n1%n .AND. &
                elemA%n3%n == elemB%n2%n) THEN
          elemA%e2 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                elemA%n3%n == elemB%n1%n) THEN
          elemA%e2 => elemB
          elemB%e2 => elemA

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

        END IF

      END IF

      !Check direction 3
      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
        IF     (elemA%n3%n == elemB%n1%n .AND. &
                elemA%n1%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

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

        END IF

      END IF

    END SUBROUTINE connectTriaEdge

END MODULE moduleMesh2DCyl