submodule(moduleMesh) boundaryEM
  CONTAINS
    ! Dirichlet boundary condition
    ! Init
    module SUBROUTINE initDirichlet(self, config, object)
      use json_module
      USE moduleRefParam, ONLY: Volt_ref
      use moduleErrors
      IMPLICIT NONE

      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object
      REAL(8):: potential
      logical:: found

      select type(self)
      type is(boundaryEMDirichlet)
        call config%get(object // '.potential', potential, found)
        if (.not. found) then
          call criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readBoundaryEM')

        end if

        self%potential = potential / Volt_ref

      end select

    end subroutine initDirichlet

    ! Apply
    module SUBROUTINE applyDirichlet(self, vectorF)
      IMPLICIT NONE

      CLASS(boundaryEMDirichlet), INTENT(in):: self
      REAL(8), INTENT(inout):: vectorF(:)
      integer:: n

      DO n = 1, self%nNodes
        associate(node => self%nodes(n)%obj)
          node%emData%phi  = self%potential
          vectorF(node%n)  = node%emData%phi

        end associate

      END DO

    END SUBROUTINE applyDirichlet

    ! Dirichlet time-dependent boundary condition
    ! Init
    module subroutine initDirichletTime(self, config, object)
      use json_module
      use moduleRefParam, ONLY: Volt_ref, ti_ref
      use moduleErrors
      use moduleOutput, only: path
      implicit none

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object
      real(8):: potential
      character(:), allocatable:: temporalProfile
      logical:: found
      character(:), allocatable:: temporalProfilePath

      select type(self)
      type is(boundaryEMDirichletTime)
        call config%get(object // '.potential', potential, found)
        if (.not. found) then
          call criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
                             'readBoundaryEM')

        end if

        call config%get(object // '.temporalProfile', temporalProfile, found)
        if (.not. found) then
          call criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
                             'readBoundaryEM')

        end if
        temporalProfilePath = path // temporalProfile

        self%potential = potential / Volt_ref

        call self%temporalProfile%init(temporalProfile)

        call self%temporalProfile%convert(1.D0/ti_ref, 1.D0)

        self%update => updateDirichletTime

      end select

    END SUBROUTINE initDirichletTime

    ! Apply
    module SUBROUTINE applyDirichletTime(self, vectorF)
      USE moduleMesh
      USE moduleCaseParam, ONLY: timeStep, tauMin
      IMPLICIT NONE

      CLASS(boundaryEMDirichletTime), INTENT(in):: self
      REAL(8), INTENT(inout):: vectorF(:)
      integer:: n

      DO n = 1, self%nNodes
        self%nodes(n)%obj%emData%phi  = self%potential * self%timeFactor
        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi

      END DO

    END SUBROUTINE applyDirichletTime

    ! Update
    subroutine updateDirichletTime(self)
      implicit none

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

      select type(self)
      type is(boundaryEMDirichletTime)
        self%timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)

      end select

    end subroutine updateDirichletTime

    ! Neumann constant value
    ! Init
    module subroutine initNeumann(self, config, object)
      use json_module
      use moduleRefParam, only: EF_ref
      use moduleErrors, only: criticalError
      implicit none

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object
      real(8):: electricField
      logical:: found

      select type(self)
      type is(boundaryEMNeumann)
        call config%get(object // '.electricField', electricField, found)
        if (.not. found) then
          call criticalError('Required parameter "electricField" for Neumann boundary condition not found', 'readBoundaryEM')

        end if

        self%electricField = electricField / EF_ref

      end select

    end subroutine initNeumann

    ! Apply
    module subroutine applyNeumann(self, vectorF)
      implicit none

      class(boundaryEMNeumann), intent(in):: self
      real(8), intent(inout):: vectorF(:)
      integer:: n

      do n = 1, self%nNodes
        associate(node => self%nodes(n)%obj)
          vectorF(node%n) = vectorF(node%n) - self%electricField/node%v

        end associate

      end do

    end subroutine applyNeumann

    ! Floating surface
    ! Init
    module subroutine initFloating(self, config, object)
      use json_module
      use moduleRefParam, only: Volt_ref, L_ref
      use moduleConstParam, only: eps_0
      use moduleErrors
      implicit none

      class(boundaryEMGeneric), allocatable, intent(inout):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object
      real(8):: potential, capacitance
      logical:: found

      select type(self)
      type is(boundaryEMFloating)
        call config%get(object // '.potential', potential, found)
        if (.not. found) then
          call criticalError('Required parameter "potential" for Floating boundary condition not found', &
                             'readBoundaryEM')

        end if

        self%potential = potential / Volt_ref

        call config%get(object // '.capacitance', capacitance, found)
        if (.not. found) then
          call warningError('Parameter "capacitance" for floating boundary condition not found. Using the default value of 1 nF.')
          capacitance = 1.0d-9

        end if
    
        ! Inverse of non-dimensional capacitance
        self%invC = 1.0d0 / (capacitance / (eps_0 * L_ref))

        self%charge = 0.0d0

        self%update => updateFloating

        self%print  => writeFloating

        allocate(self%edges(0))

      end select

    end subroutine initFloating
    
    ! Apply
    module subroutine applyFloating(self, vectorF)
      use moduleMesh
      implicit none

      class(boundaryEMFloating), intent(in):: self
      real(8), intent(inout):: vectorF(:)
      integer:: n

      do n = 1, self%nNodes
        self%nodes(n)%obj%emData%phi  = self%potential
        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi

      end do

    end subroutine applyFloating

    ! Update
    subroutine updateFloating(self)
      use moduleCaseParam, only: tauMin
      use moduleRefParam, only: Vol_ref, n_ref
      implicit none

      class(boundaryEMGeneric), intent(inout):: self
      integer:: e, n, s
      integer, allocatable:: nodes(:)
      class(meshEdge), pointer:: edge
      real(8), allocatable:: mom_nodes(:)
      class(meshNode), pointer:: node
      real(8):: mom_center, edgeDensityCurrent, chargeTau

      select type(self)
      type is(boundaryEMFloating)
        do e = 1, self%nEdges
          edge => self%edges(e)%obj

          chargeTau = 0.0d0
          edgeDensityCurrent = 0.0d0

          nodes = edge%getNodes(edge%nNodes)
          allocate(mom_nodes(1:edge%nNodes))
          do s = 1, nSpecies
            mom_center = 0.0d0
            do n = 1, self%nNodes
              node => mesh%nodes(nodes(n))%obj

                ! Minus sign to get the flux exiting the surface (opposite to the normal) 
                ! Momentum is reescaled as value at the node has no data about the node volume
              mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref)

            end do

            mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)

            edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center

          end do

          chargeTau = edgeDensityCurrent * edge%surface * tauMin

          ! Accumulate charge of speceis on surface
          self%charge = self%charge + chargeTau

        end do

        self%potential = self%potential + chargeTau * self%invC

      end select

    end subroutine updateFloating

    ! Write
    subroutine writeFloating(self, fileID)
      use moduleOutput, only: fmtColReal
      use moduleConstParam, only: qe
      use moduleRefParam, only: Volt_ref, v_ref, n_ref, L_ref, ti_ref
      implicit none

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

      write(fileID, '(A)') self%name
      select type(self)
      type is(boundaryEMFloating)
        write(fileID, '(A,",",A)') 'Total charge (C)', 'Potential (V)'
        write(fileID, '(*('//fmtColReal//'))') self%charge*qe*v_ref*n_ref*L_ref**2*ti_ref, self%potential*Volt_ref

      end select

    end subroutine writeFloating

    ! Free current
    ! Init
    module subroutine initFreeCurrent(self, config, object)
      use json_module
      use moduleErrors, only: warningError
      implicit none

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

      select type(self)
      type is(boundaryEMFreeCurrent)
        allocate(self%edges(0))

        self%counter = 0

        call config%get(object // '.every', every, found)
        if (found) then
          self%every = every

        else
          self%every = 50
          call warningError('Using "every" = 50 for Free Current boundary')

        end if

        self%update => updateFreeCurrent
        self%print  => writeFreeCurrent

      end select

    end subroutine initFreeCurrent

    ! Apply
    module subroutine applyFreeCurrent(self, vectorF)
      implicit none

      class(boundaryEMFreeCurrent), intent(in):: self
      real(8), intent(inout):: vectorF(:)
      integer:: e, n
      integer, allocatable:: nodes(:)
      real(8), allocatable:: fPsi(:)

      do e = 1, self%nEdges
        associate(edge => self%edges(e)%obj)
          fPsi = edge%fPsi(edge%centerXi(), edge%nNodes)

          nodes = edge%getNodes(edge%nNodes)

          do n = 1, edge%nNodes
            associate(node => mesh%nodes(nodes(n))%obj)
              ! Assign to each node the corresponding weight of the Neumann BC
              vectorF(node%n) = vectorF(node%n) - fPsi(n) * self%electricField(e)/node%v

            end associate

          end do

        end associate

        deallocate(fPsi, nodes)

      end do

    end subroutine applyFreeCurrent

    ! Update
    subroutine updateFreeCurrent(self)
      use moduleCaseParam, only: tauMin
      use moduleRefParam, only: Vol_ref, n_ref
      implicit none

      class(boundaryEMGeneric), intent(inout):: self
      integer:: e, n, s
      integer, allocatable:: nodes(:)
      real(8), allocatable:: mom_nodes(:)
      class(meshNode), pointer:: node
      real(8):: currentDensity_center

      select type(self)
      type is(boundaryEMFreeCurrent)
        self%counter = self%counter + 1

        do e=1, self%nEdges
          currentDensity_center    = 0.0d0

          associate(edge => self%edges(e)%obj)

            nodes = edge%getNodes(edge%nNodes)
            allocate(mom_nodes(1:edge%nNodes))

            do s = 1, nSpecies
              do n = 1, self%nNodes
                node => mesh%nodes(nodes(n))%obj

                ! Minus sign to get the flux exiting the surface (opposite to the normal) 
                ! Momentum is reescaled as value at the node has no data about the node volume
                mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref)

              end do

              currentDensity_center = currentDensity_center + qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)

            end do

            self%deltaElectricField(e) = self%deltaElectricField(e) + currentDensity_center * tauMin

            deallocate(nodes, mom_nodes)

          end associate


        end do

        if (self%counter == self%every) then
          self%electricField = self%electricField - self%deltaElectricField

          self%counter = 0
          self%deltaElectricField = 0.0d0

        end if

      end select

    end subroutine updateFreeCurrent

    ! Write
    subroutine writeFreeCurrent(self, fileID)
      ! use moduleOutput, only: fmtColReal
      use moduleRefParam, only: EF_ref
      use moduleOutput, only: fmtColInt, fmtColReal
      implicit none

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

      write(fileID, '(A)') self%name
      select type(self)
      type is(boundaryEMFreeCurrent)
        do e = 1, self%nEdges
          write(fileID, '(A,",",A)') 'Edge id', 'Electric Field (V m^-1)'
          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%electricField(e)*EF_ref

        end do

      end select

    end subroutine writeFreeCurrent

    ! Get the index of the boundary based on the name
    module function boundaryEMName_to_Index(boundaryName) result(bp)
      use moduleErrors
      implicit none

      character(:), allocatable:: boundaryName
      integer:: bp
      integer:: b

      bp = 0
      do b = 1, nBoundariesEM
        if (boundaryName == boundariesEM(b)%obj%name) then
          bp = boundariesEM(b)%obj%n

        end if

      end do

      if (bp == 0) then
        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryEMName_to_Index')
      
      end if
    
    end function boundaryEMName_to_Index

    ! Update all EM boundaries
    module subroutine boundariesEM_update()
      implicit none

      integer:: b

      do b = 1, nBoundariesEM
        if (associated(boundariesEM(b)%obj%update)) then
          call boundariesEM(b)%obj%update()

        end if

      end do

    end subroutine boundariesEM_update

    ! Write all EM boundaries
    module subroutine boundariesEM_write()
      use moduleCaseparam, only: timeStep
      use moduleOutput, only: fileID_boundaryEM, &
                              formatFileName, &
                              informFileCreation,&
                              boundaryEMOutput,&
                              generateFilePath, &
                              prefix
      implicit none

      integer:: b
      character(:), allocatable:: fileName

      if (boundaryEMOutput) then
        fileName = formatFileName(prefix, 'boundariesEM', 'csv', timeStep)
        call informFileCreation(fileName)
        open(fileID_boundaryEM, file = generateFilePath(fileName))

        do b = 1, nBoundariesEM
          if (associated(boundariesEM(b)%obj%update)) then
            call boundariesEM(b)%obj%print(fileID_boundaryEM)

          end if

        end do

        close(fileID_boundaryEM)

      end if

    end subroutine boundariesEM_write

end submodule boundaryEM