It seems to work, but I don't think this is quite useful for the manuscript

src/modules/init/moduleInput.f90

@@ -1367,6 +1367,8 @@ MODULE moduleInput
               ! If the boundary for the species is linked to the one analysing, add the edges
               if (associated(physicalSurfaces(ps)%EM, bound)) then
                 bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes]
+
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
               end if
 
             end do

src/modules/mesh/moduleMesh.f90

@@ -964,6 +964,7 @@ MODULE moduleMesh
     type(meshNodePointer), allocatable:: nodes(:)
 
     procedure(updateEM_interface), pointer, pass:: update => null()
+    procedure(printEM_interface),  pointer, pass:: print  => null()
     contains
       procedure(applyEM_interface), deferred, pass:: apply
 
@@ -987,6 +988,15 @@ MODULE moduleMesh
 
     end subroutine updateEM_interface
 
+    ! Write the values of the particle boundary model
+    subroutine printEM_interface(self, fileID)
+      import boundaryEMGeneric
+
+      class(boundaryEMGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+
+    end subroutine printEM_interface
+
   end interface
 
   ! Extended types
@@ -1049,6 +1059,8 @@ MODULE moduleMesh
   ! Floating: Floating surface, ie, zero net current
   type, extends(boundaryEMGeneric):: boundaryEMFloating
     real(8):: potential
+    real(8):: invC ! Inverse of the capacitance of the surface
+    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
       ! boundaryEMGeneric DEFERRED PROCEDURES
       procedure, pass:: apply => applyFloating
@@ -1097,6 +1109,10 @@ MODULE moduleMesh
 
     end subroutine boundariesEM_update
 
+    module subroutine boundariesEM_write()
+
+    end subroutine boundariesEM_write
+
   end interface
 
   ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -127,14 +127,15 @@ submodule(moduleMesh) boundaryEM
     ! Init
     module subroutine initFloating(self, config, object)
       use json_module
-      use moduleRefParam, only: Volt_ref
+      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
+      real(8):: potential, capacitance
       logical:: found
 
       select type(self)
@@ -148,8 +149,20 @@ submodule(moduleMesh) boundaryEM
 
         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%update => updateFloating
 
+        allocate(self%edges(0))
+
       end select
 
     end subroutine initFloating
@@ -174,40 +187,60 @@ submodule(moduleMesh) boundaryEM
     ! Update
     subroutine updateFloating(self)
       use moduleMesh, only: qSpecies, meshNode
+      use moduleCaseParam, only: tauMin
       implicit none
 
       class(boundaryEMGeneric), intent(inout):: self
-      integer:: n, s
+      integer:: e, n, s
       integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: mom_nodes(:)
       class(meshNode), pointer:: node
-      real(8):: momNode, momNodeTotal, totalCurrent
+      real(8):: mom_center, edgeDensityCurrent, totalCurrent
 
       totalCurrent = 0.0d0
       select type(self)
       type is(boundaryEMFloating)
-        do n = 1, self%nNodes
+        do e = 1, size(self%edges)
-          momNodeTotal = 0.0d0
+          edge => self%edges(e)%obj
 
-          node => self%nodes(n)%obj
+          edgeDensityCurrent = 0.0d0
 
+          nodes = edge%getNodes(edge%nNodes)
+          allocate(mom_nodes(1:edge%nNodes))
           do s = 1, nSpecies
-            print *, s, node%output(s)%mom
+            mom_center = 0.0d0
-            momNode = norm2(node%output(s)%mom)
+            do n = 1, self%nNodes
+              node => mesh%nodes(nodes(n))%obj
 
-            momNodeTotal = momNodeTotal + qSpecies(s) * momNode
+              ! Minus sign to get the flux exiting the edge
+              mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)
+
+            end do
+
+            mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+            
+            edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center
 
           end do
 
-          totalCurrent = totalCurrent + momNodeTotal
+          totalCurrent = totalCurrent + edgeDensityCurrent * edge%surface
 
         end do
 
-        print *, totalCurrent
+        self%potential = self%potential +  self%invC * totalCurrent * tauMin
+
+        ! print *, totalCurrent, self%potential
 
       end select
 
     end subroutine updateFloating
 
+    ! Write
+    subroutine writeFloating()
+
+    end subroutine writeFloating
+
     ! Get the index of the boundary based on the name
     module function boundaryEMName_to_Index(boundaryName) result(bp)
       use moduleErrors

src/modules/mesh/moduleMesh@elements.f90

@@ -176,6 +176,14 @@ submodule(moduleMesh) elements
 
           END DO
 
+        TYPE IS(boundaryEMFloating)
+          DO n = 1, boundary%nNodes
+            ni = boundary%nodes(n)%obj%n
+            self%K(ni, :)  = 0.D0
+            self%K(ni, ni) = 1.D0
+
+          END DO
+
         END SELECT
 
       END DO