!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
  USE moduleList
  USE moduleOutput, only: outputNode, emNode
  USE moduleCollisions
  use moduleSpecies, only: nSpecies

  ! MESH ELEMENTS
  !Generic mesh element
  TYPE, PUBLIC, ABSTRACT:: meshElement
    !Index
    INTEGER:: n = 0

  END TYPE meshElement

  !Generic Node element
  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshNode
    !Node volume
    REAL(8):: v = 0.D0
    !Output values
    TYPE(outputNode), ALLOCATABLE:: output(:)
    TYPE(emNode):: emData
    !Lock indicator for scattering
    INTEGER(KIND=OMP_LOCK_KIND):: lock
    CONTAINS
      !DEFERED PROCEDURES
      PROCEDURE(initNode_interface),    DEFERRED, PASS:: init
      PROCEDURE(getCoord_interface),    DEFERRED, PASS:: getCoordinates
      !GENERIC PROCEDURES
      PROCEDURE,                                  PASS:: resetOutput

  END TYPE meshNode

  interface
    pure module subroutine resetOutput(self)
      class(meshNode), intent(inout):: self

    end subroutine resetOutput

  end interface

  ABSTRACT INTERFACE
    !Interface of init a node (3D generic coordinates)
    SUBROUTINE initNode_interface(self, n, r)
      IMPORT:: meshNode
      CLASS(meshNode), INTENT(out):: self
      INTEGER, INTENT(in):: n
      REAL(8), INTENT(in):: r(1:3)

    END SUBROUTINE initNode_interface

    !Interface to get coordinates from node
    PURE FUNCTION getCoord_interface(self) RESULT(r)
      IMPORT:: meshNode
      CLASS(meshNode), INTENT(in):: self
      REAL(8):: r(1:3)

    END FUNCTION getCoord_interface

  END INTERFACE

  !Containers for nodes in the mesh
  TYPE:: meshNodeCont
    CLASS(meshNode), ALLOCATABLE:: obj

  END TYPE meshNodeCont

  ! Array of pointers to nodes.
  TYPE:: meshNodePointer
    CLASS(meshNode), POINTER:: obj
    contains
      procedure, pass:: meshNodePointer_equal_type_type
      procedure, pass:: meshNodePointer_equal_type_int
      generic:: operator(==) => meshNodePointer_equal_type_type, meshNodePointer_equal_type_int

  END TYPE meshNodePointer

  interface
    module subroutine meshNodePointer_add(self, node)
      class(meshNodePointer), allocatable, intent(inout):: self(:)
      integer, intent(in):: node

    end subroutine meshNodePointer_add

    module function meshNodePointer_equal_type_type(self, other) result(isEqual)
      class(meshNodePointer), intent(in):: self, other
      logical:: isEqual

    end function meshNodePointer_equal_type_type

    module function meshNodePointer_equal_type_int(self, other) result(isEqual)
      class(meshNodePointer), intent(in):: self
      integer, intent(in):: other
      logical:: isEqual

    end function meshNodePointer_equal_type_int

  end interface

  !Parent of Edge element
  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge
    !Nomber of nodes in the edge
    INTEGER:: nNodes
    !Connectivity to cells
    CLASS(meshCell), POINTER:: e1 => NULL(), e2 => NULL()
    !Connectivity to cells in meshColl
    CLASS(meshCell), POINTER:: eColl => NULL()
    !Normal vector
    REAL(8):: normal(1:3)
    ! Surface of edge
    REAL(8):: surface = 0.D0
    !Pointer to boundary per species
    TYPE(boundaryParticlePointer), allocatable:: boundariesParticle(:)
    !Physical surface for the edge
    CONTAINS
      !DEFERED PROCEDURES
      PROCEDURE(initEdge_interface),         DEFERRED, PASS:: init
      PROCEDURE(getNodesEdge_interface),     DEFERRED, PASS:: getNodes
      PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
      PROCEDURE(randPosEdge_interface),      DEFERRED, PASS:: randPos
      procedure(centerEdge_interface),       deferred, pass:: center
      procedure(centerXiEdge_interface),     deferred, nopass:: centerXi
      PROCEDURE(fPsi_interface),             DEFERRED, NOPASS:: fPsi
      !Gather value and spatial derivative on the nodes at position Xi
      PROCEDURE, PASS, PRIVATE:: gatherF_edge_scalar
      GENERIC:: gatherF  => gatherF_edge_scalar

  END TYPE meshEdge

  interface
    pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) result(f)
      class(meshEdge), intent(in):: self
      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8), intent(in):: valNodes(1:nNodes)
      real(8):: f

    end function gatherF_edge_scalar

  end interface

  ABSTRACT INTERFACE
    !Inits the edge parameters
    subroutine initEdge_interface(self, n, p, ps)
      use moduleSpecies, only:nSpecies
      import:: meshEdge

      class(meshEdge), intent(out):: self
      integer, intent(in):: n
      integer, intent(in):: p(:)
      integer, intent(in):: ps

    end subroutine initEdge_interface

    !Get nodes index from node
    PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n)
      IMPORT:: meshEdge
      CLASS(meshEdge), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      INTEGER:: n(1:nNodes)

    END FUNCTION getNodesEdge_interface

    !Returns the intersecction between an edge and a line defined by point r0
    PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r)
      IMPORT:: meshEdge
      CLASS(meshEdge), INTENT(in):: self
      REAL(8), INTENT(in), DIMENSION(1:3):: r0
      REAL(8):: r(1:3)

    END FUNCTION intersectionEdge_interface

    !Returns a random position in the edge
    FUNCTION randPosEdge_interface(self) RESULT(r)
      IMPORT:: meshEdge
      CLASS(meshEdge), INTENT(in):: self
      REAL(8):: r(1:3)

    END FUNCTION randPosEdge_interface

    ! Returns the center point of an edge
    FUNCTION centerEdge_interface(self) RESULT(r)
      IMPORT:: meshEdge
      CLASS(meshEdge), INTENT(in):: self
      REAL(8):: r(1:3)

    END FUNCTION centerEdge_interface

    ! Returns the center point of an edge in natural coordinates
    FUNCTION centerXiEdge_interface() RESULT(Xi)
      REAL(8):: Xi(1:3)

    END FUNCTION centerXiEdge_interface

  END INTERFACE

  !Containers for edges in the mesh
  TYPE:: meshEdgeCont
    CLASS(meshEdge), ALLOCATABLE:: obj

  END TYPE meshEdgeCont

  ! Array of pointers to edges.
  type:: meshEdgePointer
    class(meshEdge), pointer:: obj
    contains
      procedure, pass:: meshEdgePointer_equal_type_type
      procedure, pass:: meshEdgePointer_equal_type_int
      generic:: operator(==) => meshEdgePointer_equal_type_type, meshEdgePointer_equal_type_int

  end type meshEdgePointer

  interface
    module subroutine meshEdgePointer_add(self, edge)
      class(meshEdgePointer), allocatable, intent(inout):: self(:)
      integer, intent(in):: edge

    end subroutine meshEdgePointer_add

    module function meshEdgePointer_equal_type_type(self, other) result(isEqual)
      class(meshEdgePointer), intent(in):: self, other
      logical:: isEqual

    end function meshEdgePointer_equal_type_type

    module function meshEdgePointer_equal_type_int(self, other) result(isEqual)
      class(meshEdgePointer), intent(in):: self
      integer, intent(in):: other
      logical:: isEqual

    end function meshEdgePointer_equal_type_int

  end interface

  !Parent of cell element
  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshCell
    !Number of nodes in the cell
    INTEGER:: nNodes
    !Maximum collision rate
    REAL(8), ALLOCATABLE:: sigmaVrelMax(:)
    !Arrays for counting number of collisions
    TYPE(tallyCollisions), ALLOCATABLE:: tallyColl(:)
    !Volume
    REAL(8):: volume = 0.D0
    !List of particles inside the volume
    TYPE(listNode), ALLOCATABLE:: listPart_in(:)
    !Lock indicator for listPart_in
    INTEGER(KIND=OMP_LOCK_KIND):: lock
    !Total weight of particles inside cell
    REAL(8), ALLOCATABLE:: totalWeight(:)
    CONTAINS
      !DEFERRED PROCEDURES
      !Init the cell
      PROCEDURE(initCell_interface),         DEFERRED, PASS::   init
      !Get the index of the nodes
      PROCEDURE(getNodesCell_interface),     DEFERRED, PASS::   getNodes
      !Calculate random position on the cell
      PROCEDURE(randPosCell_interface),      DEFERRED, PASS::   randPos
      !Obtain functions and values of cell natural functions
      PROCEDURE(fPsi_interface),             DEFERRED, NOPASS:: fPsi
      PROCEDURE(dPsi_interface),             DEFERRED, NOPASS:: dPsi
      PROCEDURE(partialDer_interface),       DEFERRED, PASS::   partialDer
      PROCEDURE(detJac_interface),           DEFERRED, NOPASS:: detJac
      PROCEDURE(invJac_interface),           DEFERRED, NOPASS:: invJac
      !Procedures to get specific values in the node
      PROCEDURE(gatherCellArray_interface),  DEFERRED, PASS::   gatherElectricField
      PROCEDURE(gatherCellArray_interface),  DEFERRED, PASS::   gatherMagneticField
      !Compute K and F to solve PDE on the mesh
      PROCEDURE(elemK_interface),            DEFERRED, PASS::   elemK
      PROCEDURE(elemF_interface),            DEFERRED, PASS::   elemF
      !Check if particle is inside the cell
      PROCEDURE(inside_interface),           DEFERRED, NOPASS:: inside
      !Convert physical coordinates (r) into logical coordinates (Xi)
      PROCEDURE(phy2log_interface),          DEFERRED, PASS::   phy2log
      !Returns the neighbour element based on particle position outside the cell
      PROCEDURE(neighbourElement_interface), DEFERRED, PASS::   neighbourElement
      !Scatter properties of particles on cell nodes
      PROCEDURE, PASS:: scatter
      !Subroutine to find in which cell a particle is located
      PROCEDURE, PASS:: findCell
      !Gather value and spatial derivative on the nodes at position Xi
      procedure, pass, private:: gatherF_cell_scalar
      procedure, pass, private:: gatherF_cell_array
      procedure, pass, private:: gatherDF_cell_scalar
      generic:: gatherF  => gatherF_cell_scalar, gatherF_cell_array
      generic:: gatherDF => gatherDF_cell_scalar

  END TYPE meshCell

  interface
    recursive module subroutine findCell(self, part, oldCell)
      use moduleSpecies, only: particle

      class(meshCell), intent(inout):: self
      class(particle), intent(inout), target:: part
      class(meshCell), optional, intent(in):: oldCell

    end subroutine findCell

    module subroutine scatter(self, nNodes, part)
      use moduleSpecies, only: particle

      class(meshCell), intent(inout):: self
      integer, intent(in):: nNodes
      class(particle), intent(in):: part

    end subroutine scatter

    pure module function gatherF_cell_scalar(self, Xi, nNodes, valNodes) result(f)
      class(meshCell), intent(in):: self
      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8), intent(in):: valNodes(1:nNodes)
      real(8):: f

    end function gatherF_cell_scalar

    pure module function gatherF_cell_array(self, Xi, nNodes, valNodes) result(f)
      class(meshCell), intent(in):: self
      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8), intent(in):: valNodes(1:nNodes, 1:3)
      real(8):: f(1:3)

    end function gatherF_cell_array

    pure module function gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)

      class(meshCell), intent(in):: self
      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8), intent(in):: valNodes(1:nNodes)
      real(8):: df(1:3)

    end function gatherDF_cell_scalar

  end interface

  ABSTRACT INTERFACE
    SUBROUTINE initCell_interface(self, n, p, nodes)
      IMPORT:: meshCell, meshNodeCont
      CLASS(meshCell), INTENT(out):: self
      INTEGER, INTENT(in):: n
      INTEGER, INTENT(in):: p(:)
      TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)

    END SUBROUTINE initCell_interface

    PURE FUNCTION getNodesCell_interface(self, nNodes) RESULT(n)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      INTEGER:: n(1:nNodes)

    END FUNCTION getNodesCell_interface

    FUNCTION randPosCell_interface(self) RESULT(r)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      REAL(8):: r(1:3)

    END FUNCTION randPosCell_interface

    PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi)
      REAL(8), INTENT(in):: Xi(1:3)
      INTEGER, INTENT(in):: nNodes
      REAL(8):: fPsi(1:nNodes)

    END FUNCTION fPsi_interface

    PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi)
      REAL(8), INTENT(in):: Xi(1:3)
      INTEGER, INTENT(in):: nNodes
      REAL(8):: dPsi(1:3, 1:nNodes)

    END FUNCTION dPsi_interface

    PURE FUNCTION partialDer_interface(self, nNodes, dPsi) RESULT(pDer)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes)
      REAL(8):: pDer(1:3, 1:3)

    END FUNCTION partialDer_interface

    PURE FUNCTION detJac_interface(pDer) RESULT(dJ)
      REAL(8), INTENT(in):: pDer(1:3,1:3)
      REAL(8):: dJ

    END FUNCTION detJac_interface

    PURE FUNCTION invJac_interface(pDer) RESULT(invJ)
      REAL(8), INTENT(in):: pDer(1:3,1:3)
      REAL(8):: invJ(1:3,1:3)

    END FUNCTION invJac_interface

    PURE FUNCTION gatherCellArray_interface(self, Xi) RESULT(array)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      REAL(8):: array(1:3)

    END FUNCTION gatherCellArray_interface

    PURE FUNCTION elemK_interface(self, nNodes) RESULT(localK)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8):: localK(1:nNodes,1:nNodes)

    END FUNCTION elemK_interface

    PURE FUNCTION elemF_interface(self, nNodes, source) RESULT(localF)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: source(1:nNodes)
      REAL(8):: localF(1:nNodes)

    END FUNCTION elemF_interface

    PURE FUNCTION inside_interface(Xi) RESULT(ins)
      IMPORT:: meshCell
      REAL(8), INTENT(in):: Xi(1:3)
      LOGICAL:: ins

    END FUNCTION inside_interface

    PURE FUNCTION phy2log_interface(self,r) RESULT(Xi)
      IMPORT:: meshCell
      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)

    END FUNCTION phy2log_interface

    SUBROUTINE neighbourElement_interface(self, Xi, neighbourElement)
      IMPORT:: meshCell, meshElement
      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement

    END SUBROUTINE neighbourElement_interface

  END INTERFACE

  !Containers for cells in the mesh
  TYPE:: meshCellCont
    CLASS(meshCell), ALLOCATABLE:: obj

  END TYPE meshCellCont

  !Generic mesh type
  TYPE, ABSTRACT:: meshGeneric
    !Dimension of the mesh
    INTEGER:: dimen
    !Geometry of the mesh
    CHARACTER(:), ALLOCATABLE:: geometry
    !Number of elements
    INTEGER:: numNodes, numCells
    !Array of nodes
    TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
    !Array of cell elements
    TYPE(meshCellCont),  ALLOCATABLE:: cells(:)
    !PROCEDURES SPECIFIC OF FILE TYPE
    PROCEDURE(readMesh_interface),    POINTER, PASS::   readMesh     => NULL()
    PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial  => NULL()
    PROCEDURE(connectMesh_interface), POINTER, PASS::   connectMesh  => NULL()
    PROCEDURE(printColl_interface),   POINTER, PASS::   printColl    => NULL()
    CONTAINS
      !GENERIC PROCEDURES
      PROCEDURE, PASS:: doCollisions

  END TYPE meshGeneric

  interface
    module subroutine doCollisions(self)
      class(meshGeneric), intent(inout), target:: self

    end subroutine doCollisions

    module function findCellBrute(self, r) result(nVol)
      class(meshGeneric), intent(in):: self
      real(8), dimension(1:3), intent(in):: r
      integer:: nVol

    end function findCellBrute

  end interface

  ABSTRACT INTERFACE
    !Reads the mesh from a file
    SUBROUTINE readMesh_interface(self, filename)
      IMPORT meshGeneric
      CLASS(meshGeneric), INTENT(inout):: self
      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename

    END SUBROUTINE readMesh_interface

    SUBROUTINE readInitial_interface(filename, density, velocity, temperature)
      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature

    END SUBROUTINE readInitial_interface

    !Connects cell and edges to the mesh
    SUBROUTINE connectMesh_interface(self)
      IMPORT meshGeneric
      CLASS(meshGeneric), INTENT(inout):: self

    END SUBROUTINE connectMesh_interface

    !Prints number of collisions in each cell
    SUBROUTINE printColl_interface(self)
      IMPORT meshGeneric
      CLASS(meshGeneric), INTENT(in):: self

    END SUBROUTINE printColl_interface

  END INTERFACE

  !Particle mesh
  TYPE, EXTENDS(meshGeneric), PUBLIC:: meshParticles
    INTEGER:: numEdges
    !Array of boundary elements
    TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
    !Global stiffness matrix
    REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
    !Permutation matrix for P L U factorization
    INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
    !PROCEDURES SPECIFIC OF FILE TYPE
    PROCEDURE(printOutput_interface),  POINTER, PASS:: printOutput  => NULL()
    PROCEDURE(printEM_interface),      POINTER, PASS:: printEM      => NULL()
    PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
    CONTAINS
      !GENERIC PROCEDURES
      PROCEDURE, PASS:: constructGlobalK
      PROCEDURE, PASS:: doCoulomb

  END TYPE meshParticles

  interface
    module subroutine constructGlobalK(self)
      class(meshParticles), intent(inout):: self

    end subroutine constructGlobalK

    module subroutine doCoulomb(self)
      class(meshParticles), intent(in), target:: self

    end subroutine doCoulomb

  end interface

  ABSTRACT INTERFACE
    !Prints Species data
    SUBROUTINE printOutput_interface(self)
      IMPORT meshParticles
      CLASS(meshParticles), INTENT(in):: self

    END SUBROUTINE printOutput_interface

    !Prints EM info
    SUBROUTINE printEM_interface(self)
      IMPORT meshParticles
      CLASS(meshParticles), INTENT(in):: self

    END SUBROUTINE printEM_interface

    !Perform Coulomb Scattering
    !Prints average values
    SUBROUTINE printAverage_interface(self)
      IMPORT meshParticles
      CLASS(meshParticles), INTENT(in):: self

    END SUBROUTINE printAverage_interface

  END INTERFACE

  TYPE(meshParticles), TARGET:: mesh

  !Collision (MCC) mesh
  TYPE, EXTENDS(meshGeneric):: meshCollisions

  END TYPE meshCollisions

  TYPE(meshCollisions), TARGET:: meshColl

  ABSTRACT INTERFACE
    SUBROUTINE readMeshColl_interface(self, filename)
      IMPORT meshCollisions
      CLASS(meshCollisions), INTENT(inout):: self
      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename

    END SUBROUTINE readMeshColl_interface

    SUBROUTINE connectMeshColl_interface(self)
      IMPORT meshParticles

      CLASS(meshParticles), INTENT(inout):: self

    END SUBROUTINE connectMeshColl_interface

  END INTERFACE

  !Pointer to mesh used for MC collisions
  CLASS(meshGeneric), POINTER:: meshForMCC => NULL()

  !Procedure to find a cell for a particle in meshColl
  PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL()

  interface
    module subroutine findCellCollMesh(part)
      use moduleSpecies
      implicit none

      type(particle), intent(inout):: part

    end subroutine findCellCollMesh

    module subroutine findCellSameMesh(part)
      use moduleSpecies

      type(particle), intent(inout):: part

    end subroutine findCellSameMesh
  end interface

  ABSTRACT INTERFACE
    SUBROUTINE findCellColl_interface(part)
      USE moduleSpecies

      TYPE(particle), INTENT(inout):: part

    END SUBROUTINE findCellColl_interface

  END INTERFACE

  !Logical to indicate if an specific mesh for MC Collisions is used
  LOGICAL:: doubleMesh = .false.
  !Logical to indicate if MCC collisions are performed
  LOGICAL:: doMCCollisions = .FALSE.
  !Logical to indicate if Coulomb scattering is performed
  LOGICAL:: doCoulombScattering = .FALSE.
  !Logica to indicate if particles have to be listed in list inside the cells
  LOGICAL:: listInCells = .FALSE.
  !Complete path for the two meshes
  CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle

  ! BOUNDARY PARTICLE DEFINITIONS
  !Generic type for boundaries
  type, abstract, public:: boundaryParticleGeneric
    integer:: n
    character(:), allocatable:: name 
    procedure(updateParticle_interface), pointer, pass:: update => null()
    procedure(printParticle_interface),  pointer, pass:: print  => null()
    contains
      procedure(applyParticle_interface), deferred, pass:: apply

  end type boundaryParticleGeneric

  abstract interface
    ! Apply the effects of the boundary to the particle
    subroutine applyParticle_interface(self, edge, part)
      use moduleSpecies
      import boundaryParticleGeneric, meshEdge

      class(boundaryParticleGeneric), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine

    ! Update the values of the particle boundary model
    subroutine updateParticle_interface(self)
      import boundaryParticleGeneric

      class(boundaryParticleGeneric), intent(inout):: self

    end subroutine updateParticle_interface

    ! Write the values of the particle boundary model
    subroutine printParticle_interface(self, fileID)
      import boundaryParticleGeneric

      class(boundaryParticleGeneric), intent(inout):: self
      integer, intent(in):: fileID

    end subroutine printParticle_interface

  end interface

  !Reflecting boundary
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryReflection
    CONTAINS
      procedure, pass:: apply => reflection_apply

  END TYPE boundaryReflection

  interface
    module subroutine reflection_apply(self, edge, part)
      use moduleSpecies

      class(boundaryReflection), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine reflection_apply

  end interface

  !Absorption boundary
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption
    CONTAINS
      procedure, pass:: apply => absorption_apply

  END TYPE boundaryAbsorption

  interface
    module subroutine absorption_apply(self, edge, part)
      use moduleSpecies

      class(boundaryAbsorption), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine absorption_apply

  end interface

  !Transparent boundary
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent
    CONTAINS
      procedure, pass:: apply => transparent_apply

  END TYPE boundaryTransparent

  interface
    module subroutine transparent_apply(self, edge, part)
      use moduleSpecies

      class(boundaryTransparent), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine transparent_apply

  end interface

  !Symmetry axis
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis
    CONTAINS
      procedure, pass:: apply => symmetryAxis_apply

  END TYPE boundaryAxis

  interface
    module subroutine symmetryAxis_apply(self, edge, part)
      use moduleSpecies

      class(boundaryAxis), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine symmetryAxis_apply

  end interface

  !Wall Temperature boundary
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryWallTemperature
    !Thermal velocity of the wall: square root(Wall temperature X specific heat)
    REAL(8):: vTh
    CONTAINS
      procedure, pass:: apply => wallTemperature_apply

  END TYPE boundaryWallTemperature

  interface
    module subroutine wallTemperature_init(boundary, T, c)
      CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat

    end subroutine wallTemperature_init

    module subroutine wallTemperature_apply(self, edge, part)
      use moduleSpecies

      class(boundaryWallTemperature), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine wallTemperature_apply

  end interface

  !Ionization boundary
  TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryIonization
    REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
    CLASS(speciesGeneric), POINTER:: species !Ion species
    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
    TYPE(table1D):: crossSection
    REAL(8):: effectiveTime
    REAL(8):: eThreshold
    REAL(8):: deltaV
    integer:: tally
    integer(KIND=OMP_LOCK_KIND):: tally_lock
    CONTAINS
      procedure, pass:: apply => ionization_apply

  END TYPE boundaryIonization

  interface
    module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      real(8), intent(in):: mImpact
      real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
      integer, intent(in):: ion
      real(8), intent(in):: effTime
      character(:), allocatable, intent(in):: crossSection
      real(8), intent(in):: eThreshold
      integer, optional, intent(in):: electronSecondary

    end subroutine ionization_init

    module subroutine ionization_apply(self, edge, part)
      use moduleSpecies

      class(boundaryIonization), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine ionization_apply

  end interface

  ! Ensures quasi-neutrality by changing the reflection coefficient
  type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality
    real(8), allocatable:: alpha(:) ! Reflection parameter
    integer:: s_incident ! species index of the incident species
    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
    contains
      procedure, pass:: apply => quasiNeutrality_apply

  end type boundaryQuasiNeutrality

  interface
    module subroutine quasiNeutrality_init(boundary, s_incident)
      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      integer, intent(in):: s_incident

    end subroutine quasiNeutrality_init

    module subroutine quasiNeutrality_apply(self, edge, part)
      use moduleSpecies

      class(boundaryQuasiNeutrality), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine quasiNeutrality_apply

  end interface

  !Boundary for quasi-neutral outflow adjusting reflection coefficient
  type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive
    real(8), allocatable:: velocity_shift(:)
    integer:: s_incident ! species index of the incident species
    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
    contains
      procedure, pass:: apply => outflowAdaptive_apply

  end type boundaryOutflowAdaptive

  interface
    module subroutine outflowAdaptive_init(boundary, s_incident)
      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      integer, intent(in):: s_incident

    end subroutine outflowAdaptive_init

    module subroutine outflowAdaptive_apply(self, edge, part)
      use moduleSpecies

      class(boundaryOutflowAdaptive), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine outflowAdaptive_apply

  end interface

  ! Generic basic boundary conditions to use internally in the code
  interface
    module subroutine genericReflection(edge, part) 
      use moduleSpecies

      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine genericReflection

    module subroutine genericAbsorption(edge, part) 
      use moduleSpecies

      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine genericAbsorption

    module subroutine genericTransparent(edge, part) 
      use moduleSpecies

      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

    end subroutine genericTransparent

  end interface

  TYPE:: boundaryParticleCont
    CLASS(boundaryParticleGeneric), ALLOCATABLE:: obj

  END TYPE boundaryParticleCont

  type:: boundaryParticlePointer
    class(boundaryParticleGeneric), pointer:: obj

  end type boundaryParticlePointer

  !Number of boundaries
  INTEGER:: nBoundariesParticle = 0
  !Array for boundaries
  type(boundaryParticleCont), allocatable, target:: boundariesParticle(:)

  interface 
    ! Get the index of the species from its name
    module function boundaryParticleName_to_Index(boundaryName) result(bp)
      character(:), allocatable:: boundaryName
      integer:: bp

    end function boundaryParticleName_to_Index

    ! Update the particle boundary models
    module subroutine boundariesParticle_update()

    end subroutine boundariesParticle_update

    ! Write data about particle boundaries
    module subroutine boundariesParticle_write()

    end subroutine boundariesParticle_write
  end interface

  ! BOUNDARY ELECTROMAGNETIC DEFINITIONS
  ! Generic type for electromagnetic boundary conditions
  type, public, abstract:: boundaryEMGeneric
    integer:: n
    character(:), allocatable:: name 
    integer:: nNodes
    type(meshNodePointer), allocatable:: nodes(:)

    procedure(updateEM_interface),        pointer, pass:: update => null()
    procedure(printEMboundary_interface), pointer, pass:: print  => null()
    contains
      procedure(applyEM_interface), deferred, pass:: apply

  end type boundaryEMGeneric

  abstract interface
    ! Apply boundary condition to the load vector for the Poission equation
    subroutine applyEM_interface(self, vectorF)
      import boundaryEMGeneric

      class(boundaryEMGeneric), intent(in):: self
      real(8), intent(inout):: vectorF(:)

    end subroutine applyEM_interface

    ! Update the values of the EM boundary model
    subroutine updateEM_interface(self)
      import boundaryEMGeneric

      class(boundaryEMGeneric), intent(inout):: self

    end subroutine updateEM_interface

    ! Write the values of the particle boundary model
    subroutine printEMboundary_interface(self, fileID)
      import boundaryEMGeneric

      class(boundaryEMGeneric), intent(inout):: self
      integer, intent(in):: fileID

    end subroutine printEMboundary_interface

  end interface

  ! Extended types
  ! Dirichlet: Constant potential value
  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
    REAL(8):: potential
    CONTAINS
      ! boundaryEMGeneric DEFERRED PROCEDURES
      PROCEDURE, PASS:: apply => applyDirichlet

  END TYPE boundaryEMDirichlet

  interface
    module subroutine initDirichlet(self, config, object)
      use json_module

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object

    end subroutine initDirichlet

    module subroutine applyDirichlet(self, vectorF)
      class(boundaryEMDirichlet), intent(in):: self
      real(8), intent(inout):: vectorF(:)

    end subroutine applyDirichlet

  end interface

  ! Dirichlet time: time-dependent potential value
  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
    real(8):: potential
    real(8):: timeFactor
    type(table1D):: temporalProfile
    contains
      ! boundaryEMGeneric DEFERRED PROCEDURES
      procedure, pass:: apply => applyDirichletTime

  END TYPE boundaryEMDirichletTime

  interface
    module subroutine initDirichletTime(self, config, object)
      use json_module

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object

    end subroutine initDirichletTime

    module subroutine applyDirichletTime(self, vectorF)
      class(boundaryEMDirichletTime), intent(in):: self
      real(8), intent(inout):: vectorF(:)

    end subroutine applyDirichletTime

  end interface

  ! Floating: Floating surface, ie, zero net current
  type, extends(boundaryEMGeneric):: boundaryEMFloating
    real(8):: potential
    real(8):: invC ! Inverse of the capacitance of the surface
    real(8):: charge ! Charge accumulated on surface
    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
    contains
      ! boundaryEMGeneric DEFERRED PROCEDURES
      procedure, pass:: apply => applyFloating

  end type boundaryEMFloating

  interface
    module subroutine initFloating(self, config, object)
      use json_module

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object

    end subroutine initFloating

    module subroutine applyFloating(self, vectorF)
      class(boundaryEMFloating), intent(in):: self
      real(8), intent(inout):: vectorF(:)

    end subroutine applyFloating

  end interface

  ! Container for boundary conditions
  TYPE:: boundaryEMCont
    CLASS(boundaryEMGeneric), ALLOCATABLE:: obj

  END TYPE boundaryEMCont

  INTEGER:: nBoundariesEM
  TYPE(boundaryEMCont), ALLOCATABLE, target:: boundariesEM(:)

  !Information of charge and reference parameters for rho vector
  REAL(8), ALLOCATABLE:: qSpecies(:)

  ! Update the EM boundary models
  interface 
    module function boundaryEMName_to_Index(boundaryName) result(bp)
      character(:), allocatable:: boundaryName
      integer:: bp

    end function boundaryEMName_to_Index

    module subroutine boundariesEM_update()

    end subroutine boundariesEM_update

    module subroutine boundariesEM_write()

    end subroutine boundariesEM_write

  end interface

  ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS
  ! Link physical surface to edges
  type:: physicalSurface
    integer:: index
    class(meshNodePointer), allocatable:: nodes(:)
    class(meshEdgePointer), allocatable:: edges(:)
    class(boundaryParticlePointer), allocatable:: particles(:)
    class(boundaryEMGeneric), pointer:: EM => null()

  end type

  integer:: nPhysicalSurfaces
  type(physicalSurface), allocatable:: physicalSurfaces(:)

  ! Get ID from array based on physical surface
  interface
    module function physicalSurface_to_index(ps) result(index)
      integer:: ps
      integer:: index

    end function physicalSurface_to_index

  end interface

END MODULE moduleMesh