!Module to solve the electromagnetic field
MODULE moduleEM
  IMPLICIT NONE

  TYPE:: boundaryEM
    CHARACTER(:), ALLOCATABLE:: typeEM
    INTEGER:: physicalSurface
    REAL(8):: potential

    CONTAINS
      PROCEDURE, PASS:: apply

  END TYPE boundaryEM

  INTEGER:: nBoundaryEM
  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)

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

  CONTAINS
    !Apply boundary conditions to the K matrix for Poisson's equation
    SUBROUTINE apply(self, edge)
      USE moduleMesh
      IMPLICIT NONE

      CLASS(boundaryEM), INTENT(in):: self
      CLASS(meshEdge):: edge
      INTEGER, ALLOCATABLE:: nodes(:)
      INTEGER:: n

      nodes = edge%getNodes()

      DO n = 1, SIZE(nodes)
        SELECT CASE(self%typeEM)
        CASE ("dirichlet")
          mesh%K(nodes(n), :)        = 0.D0
          mesh%K(nodes(n), nodes(n)) = 1.D0
        
          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential

        END SELECT

      END DO

    END SUBROUTINE

    !Assemble the source vector based on the charge density to solve Poisson's equation
    SUBROUTINE assembleSourceVector(vectorF)
      USE moduleMesh
      USE moduleRefParam
      IMPLICIT NONE

      REAL(8), INTENT(out):: vectorF(1:mesh%numNodes)
      REAL(8), ALLOCATABLE:: localF(:)
      INTEGER, ALLOCATABLE:: nodes(:)
      REAL(8), ALLOCATABLE:: rho(:)
      INTEGER:: nNodes
      INTEGER:: e, i, ni
      CLASS(meshNode), POINTER:: node

      !$OMP SINGLE
      vectorF = 0.D0
      !$OMP END SINGLE

      !$OMP DO REDUCTION(+:vectorF)
      DO e = 1, mesh%numCells
        nodes  = mesh%cells(e)%obj%getNodes()
        nNodes = SIZE(nodes)
        !Calculates charge density (rho) in element nodes
        ALLOCATE(rho(1:nNodes))
        rho = 0.D0
        DO i = 1, nNodes
          ni = nodes(i)
          node => mesh%nodes(ni)%obj
          rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))

        END DO

        !Calculates local F vector
        localF = mesh%cells(e)%obj%elemF(rho)

        !Assign local F to global F
        DO i = 1, nNodes
          ni = nodes(i)
          vectorF(ni) = vectorF(ni) + localF(i)

        END DO

        DEALLOCATE(localF)
        DEALLOCATE(nodes, rho)

      END DO
      !$OMP END DO

      !Apply boundary conditions
      !$OMP DO
      DO i = 1, mesh%numNodes
        node => mesh%nodes(i)%obj

        SELECT CASE(node%emData%type)
        CASE ("dirichlet")
          vectorF(i) = node%emData%phi

        END SELECT

      END DO
      !$OMP END DO

    END SUBROUTINE assembleSourceVector

    !Solving the Poisson equation for electrostatic potential
    SUBROUTINE solveElecField()
      USE moduleMesh
      USE moduleErrors
      IMPLICIT NONE

      INTEGER, SAVE:: INFO
      INTEGER:: n
      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
      EXTERNAL:: dgetrs

      !$OMP SINGLE
      ALLOCATE(tempF(1:mesh%numNodes))
      !$OMP END SINGLE

      CALL assembleSourceVector(tempF)

      !$OMP SINGLE
      CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
                              mesh%IPIV, tempF,  mesh%numNodes, info)
      !$OMP END SINGLE

      IF (info == 0) THEN
        !Suscessful resolution of Poission equation
        !$OMP DO
        DO n = 1, mesh%numNodes
          mesh%nodes(n)%obj%emData%phi = tempF(n)

        END DO
        !$OMP END DO

      ELSE
        !$OMP SINGLE
        CALL criticalError('Poisson equation failed', 'solveElecField')
        !$OMP END SINGLE

      END IF

      !$OMP SINGLE
      DEALLOCATE(tempF)
      !$OMP END SINGLE

    END SUBROUTINE solveElecField

END MODULE moduleEM