Boundaries for EM almost done

src/makefile

@@ -1,4 +1,4 @@
-OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o \
+OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshCommon.o $(OBJDIR)/moduleCompTime.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
 			  	$(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \

src/modules/init/moduleInput.f90

@@ -38,10 +38,15 @@ MODULE moduleInput
       CALL readSpecies(config)
       CALL checkStatus(config, "readSpecies")
 
-      !Read boundaries
+      !Read particle boundaries
-      CALL verboseError('Reading boundary conditions...')
+      CALL verboseError('Reading particle boundary conditions...')
-      CALL readBoundary(config)
+      CALL readBoundaryParticle(config)
-      CALL checkStatus(config, "readBoundary")
+      CALL checkStatus(config, "readBoundaryParticle")
+
+      ! read EM boundaries
+      CALL verboseError('Reading EM boundary conditions...')
+      CALL readBoundaryEM(config)
+      CALL checkStatus(config, "readBoundaryEM")
 
       !Read Geometry
       CALL verboseError('Reading Geometry...')
@@ -254,17 +259,7 @@ MODULE moduleInput
       IF (found) THEN
         CALL solver%initEM(EMType)
         SELECT CASE(EMType)
-        CASE("Electrostatic")
-          !Read BC
-          CALL readEMBoundary(config)
-
-        CASE("ElectrostaticBoltzmann")
-          !Read BC
-          CALL readEMBoundary(config)
-
         CASE("ConstantB")
-          !Read BC
-          CALL readEMBoundary(config)
           !Read constant magnetic field
           DO i = 1, 3
             WRITE(iString, '(i2)') i
@@ -791,7 +786,7 @@ MODULE moduleInput
     END SUBROUTINE readInteractions
 
     !Reads boundary conditions for the mesh
-    SUBROUTINE readBoundary(config)
+    SUBROUTINE readBoundaryParticle(config)
       use moduleMesh
       USE moduleErrors
       use moduleSpecies, only: nSpecies
@@ -837,7 +832,7 @@ MODULE moduleInput
       !Init the list of particles from surfaces
       CALL OMP_INIT_LOCK(partSurfaces%lock)
 
-    END SUBROUTINE readBoundary
+    END SUBROUTINE readBoundaryParticle
 
     !Read the geometry (mesh) for the case
     SUBROUTINE readGeometry(config)
@@ -1086,7 +1081,7 @@ MODULE moduleInput
 
     END SUBROUTINE readProbes
 
-    SUBROUTINE readEMBoundary(config)
+    SUBROUTINE readBoundaryEM(config)
       USE moduleMesh
       USE moduleOutput
       USE moduleErrors
@@ -1107,68 +1102,22 @@ MODULE moduleInput
       INTEGER:: info
       EXTERNAL:: dgetrf
 
-      CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
+      CALL config%info('boundaries.EM.models', found, n_children = nBoundariesEM)
       
       IF (found) THEN
-        ALLOCATE(boundaryEM(1:nBoundaryEM))
+        ALLOCATE(boundariesEM(1:nBoundariesEM))
 
       END IF
       
       DO b = 1, nBoundaryEM
         WRITE(bString, '(I2)') b
-        object = 'boundaryEM(' // TRIM(bString) // ')'
+        object = 'boundaries.EM.models(' // TRIM(bString) // ')'
 
-        CALL config%get(object // '.type', typeEM, found)
+        call boundariesEM(b)%init(config, object, b)
 
-        SELECT CASE(typeEM)
+      end do
-        CASE ("dirichlet")
-          CALL config%get(object // '.potential', potential, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
-
-          END IF
-
-          CALL config%get(object // '.physicalSurface', physicalSurface, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-
-          CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
-
-        CASE ("dirichletTime")
-          CALL config%get(object // '.potential', potential, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
-                               'readEMBoundary')
-
-          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', &
-                               'readEMBoundary')
-
-          END IF
-          temporalProfilePath = path // temporalProfile
-
-          CALL config%get(object // '.physicalSurface', physicalSurface, found)
-          IF (.NOT. found) THEN
-            CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
-                               'readEMBoundary')
-
-          END IF
-
-          CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
-
-        CASE DEFAULT
-          CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
-
-        END SELECT
-
-      END DO
 
+      ! TODO: Move this to the init of species
       ALLOCATE(qSpecies(1:nSpecies))
       DO s = 1, nSpecies
         SELECT TYPE(sp => species(s)%obj)
@@ -1182,9 +1131,10 @@ MODULE moduleInput
 
       END DO
 
+      ! TODO: Do this after the mesh has been read
       ! Modify K matrix due to boundary conditions
-      DO b = 1, nBoundaryEM
+      DO b = 1, nBoundariesEM
-        SELECT TYPE(boundary => boundaryEM(b)%obj)
+        SELECT TYPE(boundary => boundariesEM(b)%obj)
         TYPE IS(boundaryEMDirichlet)
           DO n = 1, boundary%nNodes
             ni = boundary%nodes(n)%obj%n
@@ -1208,11 +1158,11 @@ MODULE moduleInput
       !Compute the PLU factorization of K once boundary conditions have been read
       CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
       IF (info /= 0) THEN
-        CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
+        CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM')
 
       END IF
 
-    END SUBROUTINE readEMBoundary
+    END SUBROUTINE readBoundaryEM
 
     !Reads the injection of particles from the boundaries
     SUBROUTINE readInject(config)

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -199,7 +199,7 @@ MODULE moduleMeshInputGmsh2
           END SELECT
 
           !Associate boundary condition procedure.
-          bt = physicalSurface_to_id(boundaryType)
+          bt = physicalSurface_to_id(boundariesParticleLinking,boundaryType)
 
           CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex)
           DEALLOCATE(p)

src/modules/mesh/inout/text/moduleMeshInputText.f90

@@ -142,7 +142,7 @@ module moduleMeshInputText
             allocate(p(1))
             p(1) = n
             !Associate boundary condition procedure.
-            bt = physicalSurface_to_id(physicalID)
+            bt = physicalSurface_to_id(boundariesParticleLinking, physicalID)
             call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(bt)%speciesIndex)
 
             deallocate(p)

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -365,7 +365,7 @@ MODULE moduleMeshInputVTU
           END SELECT
 
           !Associate boundary condition procedure.
-          bt = physicalSurface_to_id(entitiesID(n))
+          bt = physicalSurface_to_id(boundariesParticleLinking,entitiesID(n))
 
           !Init edge
           CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex)

src/modules/mesh/moduleMesh.f90

@@ -72,9 +72,20 @@ MODULE moduleMesh
   TYPE:: meshNodePointer
     CLASS(meshNode), POINTER:: obj
     CONTAINS
+      procedure, pass:: meshNodePointer_equal
+      generic:: operator(==) => meshNodePointer_equal
 
   END TYPE meshNodePointer
 
+  interface
+    module function meshNodePointer_equal(self, other) result(isEqual)
+      class(meshNodePointer), intent(in):: self, other
+      logical:: isEqual
+
+    end function meshNodePointer_equal
+
+  end interface
+
   !Parent of Edge element
   TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge
     !Nomber of nodes in the edge
@@ -89,7 +100,6 @@ MODULE moduleMesh
     REAL(8):: surface = 0.D0
     !Pointer to boundary per species
     TYPE(boundaryParticlePointer), allocatable:: boundariesParticle(:)
-    class(boundaryEMGeneric), pointer:: boundaryEM => null()
     !Physical surface for the edge
     CONTAINS
       !DEFERED PROCEDURES
@@ -118,16 +128,17 @@ MODULE moduleMesh
 
   ABSTRACT INTERFACE
     !Inits the edge parameters
-    SUBROUTINE initEdge_interface(self, n, p, bt)
+    subroutine initEdge_interface(self, n, p, btPart, btEM)
       use moduleSpecies, only:nSpecies
-      IMPORT:: meshEdge
+      import:: meshEdge
 
-      CLASS(meshEdge), INTENT(out):: self
+      class(meshEdge), intent(out):: self
-      INTEGER, INTENT(in):: n
+      integer, intent(in):: n
-      INTEGER, INTENT(in):: p(:)
+      integer, intent(in):: p(:)
-      INTEGER, INTENT(in):: bt(1:nSpecies)
+      integer, intent(in):: btPart(1:nSpecies)
+      integer, intent(in):: btEM
 
-    END SUBROUTINE initEdge_interface
+    end subroutine initEdge_interface
 
     !Get nodes index from node
     PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n)
@@ -799,18 +810,11 @@ MODULE moduleMesh
   !Array for linking boundaries
   type(boundaryParticleLinking), allocatable, dimension(:):: boundariesParticleLinking
 
-  interface
-    module function physicalSurface_to_id(physicalSurface) result(b)
-      integer:: physicalSurface
-      integer:: b
-
-    end function physicalSurface_to_id
-
-  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(:)
 
@@ -818,6 +822,7 @@ MODULE moduleMesh
     contains
       procedure, pass::                              init   => initBoundaryEM
       procedure(applyEM_interface), deferred, pass:: apply
+      procedure, pass::                              addNodes
 
   end type boundaryEMGeneric
 
@@ -832,6 +837,12 @@ MODULE moduleMesh
 
     end subroutine initBoundaryEM
 
+    module subroutine addNodes(self, nodes)
+      class(boundaryEMGeneric), intent(inout):: self
+      type(meshNodePointer), intent(in):: nodes(:)
+
+    end subroutine addNodes
+
   end interface
 
   abstract interface
@@ -850,7 +861,7 @@ MODULE moduleMesh
 
       class(boundaryEMGeneric), intent(in):: self
 
-  end subroutine updateEM_interface
+    end subroutine updateEM_interface
 
   end interface
 
@@ -894,10 +905,34 @@ MODULE moduleMesh
 
   END TYPE boundaryEMCont
 
-  INTEGER:: nBoundaryEM
+  type:: boundaryEMLinking
-  TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
+    integer:: physicalSurface
+    class(boundaryEMGeneric), pointer:: boundary => null()
+
+  end type boundaryEMLinking
+
+  INTEGER:: nBoundariesEM
+  TYPE(boundaryEMCont), ALLOCATABLE, target:: boundariesEM(:)
 
   !Information of charge and reference parameters for rho vector
   REAL(8), ALLOCATABLE:: qSpecies(:)
 
+  ! Get ID from array based on physical surface
+  interface physicalSurface_to_id
+    module function physicalSurface_to_boundaryParticleId(boundaryArray, physicalSurface) result(b)
+      type(boundaryParticleLinking):: boundaryArray(:)
+      integer:: physicalSurface
+      integer:: b
+
+    end function physicalSurface_to_boundaryParticleId
+
+    module function physicalSurface_to_boundaryEMId(boundaryArray, physicalSurface) result(b)
+      type(boundaryEMLinking):: boundaryArray(:)
+      integer:: physicalSurface
+      integer:: b
+
+    end function physicalSurface_to_boundaryEMId
+
+  end interface physicalSurface_to_id
+
 END MODULE moduleMesh

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -1,102 +1,152 @@
 submodule(moduleMesh) boundaryEM
   CONTAINS
-    SUBROUTINE findNodes(self, physicalSurface)
+    subroutine addNodes(self, nodes)
-      USE moduleMesh
+      implicit none
-      IMPLICIT NONE
 
-      CLASS(boundaryEMGeneric), INTENT(inout):: self
+      class(boundaryEMGeneric), intent(inout):: self
-      INTEGER, INTENT(in):: physicalSurface
+      type(meshNodePointer), intent(in)::       nodes(:)
-      CLASS(meshEdge), POINTER:: edge
+      INTEGER:: n, nn, nNodes
-      INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
+      logical:: inArray
-      INTEGER:: nNodes, nodesNew
-      INTEGER:: e, n
 
-      !Temporal array to hold nodes
-      ALLOCATE(nodes(0))
 
-      ! Loop thorugh the edges and identify those that are part of the boundary
+      nNodes = size(nodes)
-      DO e = 1, mesh%numEdges
+      ! Collect all nodes that are not already in the boundary
-        edge => mesh%edges(e)%obj
-        IF (edge%physicalSurface == physicalSurface) THEN
-          ! Edge is of the right boundary index
-          ! Get nodes in the edge
-          nNodes = edge%nNodes
-          nodesEdge = edge%getNodes(nNodes)
-          ! Collect all nodes that are not already in the temporal array
-          DO n = 1, nNodes
-            IF (ANY(nodes == nodesEdge(n))) THEN
-              ! Node already in array, skip
-              CYCLE
-
-            ELSE
-              ! If not, add element to array of nodes
-              nodes = [nodes, nodesEdge(n)]
-
-            END IF
-
-          END DO
-
-        END IF
-
-      END DO
-
-      ! Point boundary to nodes
-      nNodes = SIZE(nodes)
-      ALLOCATE(self%nodes(nNodes))
-      self%nNodes = nNodes
       DO n = 1, nNodes
-        self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
+        inArray = .false.
+        ! I cannot use the procedure 'any' for this
+        do nn = 1 , size(self%nodes)
+          IF (self%nodes(nn) == nodes(n)) THEN
+            ! Node already in array
+            inArray = .true.
+            exit
 
-      END DO
+          end if
 
-    END SUBROUTINE findNodes
+        end do
+
+        if (.not. inArray) then
+          ! If not, add element to array of nodes
+          self%nodes = [self%nodes, nodes(n)]
+
+        end if
+
+      end do
+
+    end subroutine addNodes
+
+    module subroutine initBoundaryEM(self, config, object, b)
+      use json_module
+      use moduleErrors
+      implicit none
+
+      class(boundaryEMGeneric), allocatable, intent(out):: self
+      type(json_file), intent(inout):: config
+      character(:), allocatable, intent(in):: object
+      integer, intent(in):: b
+      character(:), allocatable:: bType
+      logical:: found
+
+      self%n = b
+      call config%get(object // '.name', self%name, found)
+      if (.not. found) then
+        call criticalError('Required parameter "name" for EM boundary condition not found', &
+                           'initBoundaryEM')
+
+      end if
+
+      call config%get(object // '.type', bType,     found)
+      if (.not. found) then
+        call criticalError('Required parameter "type" for EM boundary condition not found', &
+                           'initBoundaryEM')
+
+      end if
+
+      select case(bType)
+      case ("dirichlet")
+        ! Allocate boundary edge
+        ALLOCATE(boundaryEMDirichlet:: self)
+
+        CALL initDirichlet(self, config, object)
+
+      case ("dirichlettime")
+        ! Allocate boundary edge
+        ALLOCATE(boundaryEMDirichletTime:: self)
+
+        call initDirichletTime(self, config, object)
+
+      case default
+        call criticalError('Boundary type ' // bType // ' not supported', 'readBoundaryEM')
+
+      end select
+
+    end subroutine initBoundaryEM
 
     ! Initialize Dirichlet boundary condition
-    SUBROUTINE initDirichlet(self, physicalSurface, potential)
+    SUBROUTINE initDirichlet(self, config, object)
+      use json_module
       USE moduleRefParam, ONLY: Volt_ref
+      use moduleErrors
       IMPLICIT NONE
 
       CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
-      INTEGER, INTENT(in):: physicalSurface
+      type(json_file), intent(inout):: config
-      REAL(8), INTENT(in):: potential
+      character(:), allocatable, intent(in):: object
+      REAL(8):: potential
+      logical:: found
 
-      ! Allocate boundary edge
+      select type(self)
-      ALLOCATE(boundaryEMDirichlet:: 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
 
-      SELECT TYPE(self)
-      TYPE IS(boundaryEMDirichlet)
         self%potential = potential / Volt_ref
 
-        CALL findNodes(self, physicalSurface)
+      end select
 
-      END SELECT
+    end subroutine initDirichlet
-
-    END SUBROUTINE initDirichlet
 
     ! Initialize Dirichlet boundary condition
-    SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile)
+    subroutine initDirichletTime(self, config, object)
-      USE moduleRefParam, ONLY: Volt_ref, ti_ref
+      use json_module
-      IMPLICIT NONE
+      use moduleRefParam, ONLY: Volt_ref, ti_ref
+      use moduleErrors
+      implicit none
 
-      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
+      class(boundaryEMGeneric), allocatable, intent(out):: self
-      INTEGER, INTENT(in):: physicalSurface
+      type(json_file), intent(inout):: config
-      REAL(8), INTENT(in):: potential
+      character(:), allocatable, intent(in):: object
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile
+      real(8):: potential
+      character(:), allocatable:: temporalProfile
+      logical:: found
+      character(:), allocatable:: temporalProfilePath
 
-      ! Allocate boundary edge
+      select type(self)
-      ALLOCATE(boundaryEMDirichletTime:: 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
 
-      SELECT TYPE(self)
-      TYPE IS(boundaryEMDirichletTime)
         self%potential = potential / Volt_ref
 
-        CALL findNodes(self, physicalSurface)
+        call self%temporalProfile%init(temporalProfile)
 
-        CALL self%temporalProfile%init(temporalProfile)
+        call self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
 
-        CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
+      end select
-
-      END SELECT
 
     END SUBROUTINE initDirichletTime
 
@@ -138,4 +188,24 @@ submodule(moduleMesh) boundaryEM
 
     END SUBROUTINE applyDirichletTime
 
+    function physicalSurface_to_boundaryEMId(boundaryArray,physicalSurface) result(b)
+      implicit none
+
+      type(boundaryEMLinking):: boundaryArray(:)
+      integer:: physicalSurface
+      integer:: b
+      integer:: bt
+
+      do bt=1, nPhysicalSurfaces
+        if (boundaryArray(bt)%physicalSurface == physicalSurface) then
+          b = boundaryArray(bt)%boundary%n
+          
+          exit
+
+        end if
+
+      end do
+
+    end function physicalSurface_to_boundaryEMId
+
 end submodule boundaryEM

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -48,7 +48,19 @@ submodule(moduleMesh) boundaryParticle
 
       self%n = b
       CALL config%get(object // '.name', self%name, found)
+      if (.not. found) then
+        call criticalError('Required parameter "name" for EM boundary condition not found', &
+                           'initBoundaryParticle')
+
+      end if
+
       CALL config%get(object // '.type', bType,     found)
+      if (.not. found) then
+        call criticalError('Required parameter "type" for EM boundary condition not found', &
+                           'initBoundaryParticle')
+
+      end if
+
       SELECT CASE(bType)
       CASE('reflection')
         ALLOCATE(boundaryReflection:: self)
@@ -483,15 +495,16 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine genericTransparent
 
-    function physicalSurface_to_id(physicalSurface) result(b)
+    function physicalSurface_to_boundaryParticleId(boundaryArray, physicalSurface) result(b)
       implicit none
 
+      type(boundaryParticleLinking):: boundaryArray(:)
       integer:: physicalSurface
       integer:: b
       integer:: bt
 
       do bt=1, nPhysicalSurfaces
-        if (boundariesParticleLinking(bt)%physicalSurface == physicalSurface) then
+        if (boundaryArray(bt)%physicalSurface == physicalSurface) then
           b = bt
           
           exit
@@ -500,6 +513,6 @@ submodule(moduleMesh) boundaryParticle
 
       end do
 
-    end function physicalSurface_to_id
+    end function physicalSurface_to_boundaryParticleId
 
 end submodule boundaryParticle

src/modules/mesh/moduleMesh@elements.f90

@@ -18,6 +18,16 @@ submodule(moduleMesh) elements
 
     END SUBROUTINE resetOutput
 
+    module function meshNodePointer_equal(self, other) result(isEqual)
+      implicit none
+
+      class(meshNodePointer), intent(in):: self, other
+      logical:: isEqual
+
+      isEqual = self%obj%n == other%obj%n
+
+    end function meshNodePointer_equal
+
     !Constructs the global K matrix
     PURE module SUBROUTINE constructGlobalK(self)
       IMPLICIT NONE