From d211ed9a716a7db5898839af46a3b9f65d0e305e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 18 Feb 2026 15:14:59 +0100 Subject: [PATCH] Almost done with EM, not to modify the EM module itself --- src/modules/init/moduleInput.f90 | 65 ++++++----------- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 12 ++-- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 12 ++-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 10 ++- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 10 ++- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 12 ++-- .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 8 ++- .../mesh/inout/text/moduleMeshInputText.f90 | 7 +- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 8 ++- src/modules/mesh/moduleMesh.f90 | 11 ++- src/modules/mesh/moduleMesh@boundaryEM.f90 | 33 +++++++-- .../mesh/moduleMesh@boundaryParticle.f90 | 5 +- src/modules/mesh/moduleMesh@elements.f90 | 69 ++++++++++++++----- 13 files changed, 165 insertions(+), 97 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index ed34f3f..8cfcee2 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -812,15 +812,15 @@ MODULE moduleInput END DO - ! Read linking - CALL config%info('boundaries.particles.linking', found, n_children = nPhysicalSurfaces) + ! Link physical surfaces to boundaries + call config%info('boundaries.particles.linking', found, n_children = nPhysicalSurfaces) allocate(boundaryParticlesLinking(1:nPhysicalSurfaces)) allocate(speciesNames(1:nSpecies)) do b = 1, nPhysicalSurfaces write(iString, '(i2)') b object = 'boundary.particles.linking(' // trim(iString) // ')' - call config%get(object // '.physicalSurfaces', boundaryParticlesLinking(b)%physicalSurface, found) - call config%get(object // '.models', speciesNames, found) + call config%get(object // '.physicalSurface', boundaryParticlesLinking(b)%physicalSurface, found) + call config%get(object // '.model', speciesNames, found) allocate(boundaryParticlesLinking(b)%speciesIndex(1:nSpecies)) do s = 1, nSpecies boundaryParticlesLinking(b)%speciesIndex(s) = boundaryParticlesName_to_Index(speciesNames(s)) @@ -1093,14 +1093,9 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: typeEM - REAL(8):: potential - INTEGER:: physicalSurface - CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath - INTEGER:: b, s, n, ni + INTEGER:: b CHARACTER(2):: bString - INTEGER:: info - EXTERNAL:: dgetrf + character(len=100), allocatable:: modelName(:) CALL config%info('boundaries.EM.models', found, n_children = nBoundariesEM) @@ -1109,14 +1104,27 @@ MODULE moduleInput END IF - DO b = 1, nBoundaryEM - WRITE(bString, '(I2)') b + do b = 1, nBoundaryEM + write(bString, '(I2)') b object = 'boundaries.EM.models(' // TRIM(bString) // ')' call boundariesEM(b)%init(config, object, b) end do + ! Link physical surfaces to boundaries + call config%info('boundaries.EM.linking', found, n_children = nPhysicalSurfaces) + allocate(boundaryEMLinking(1:nPhysicalSurfaces)) + allocate(speciesNames(1:nSpecies)) + do b = 1, nPhysicalSurfaces + write(iString, '(i2)') b + object = 'boundary.EM.linking(' // trim(iString) // ')' + call config%get(object // '.physicalSurface', boundaryEMLinking(b)%physicalSurface, found) + call config%get(object // '.model', modelName, found) + boundaryEMLinking(b)%model => boundariesEM(boundaryEMName_to_Index(modelName(s)))%obj + + end do + ! TODO: Move this to the init of species ALLOCATE(qSpecies(1:nSpecies)) DO s = 1, nSpecies @@ -1131,37 +1139,6 @@ MODULE moduleInput END DO - ! TODO: Do this after the mesh has been read - ! Modify K matrix due to boundary conditions - DO b = 1, nBoundariesEM - SELECT TYPE(boundary => boundariesEM(b)%obj) - TYPE IS(boundaryEMDirichlet) - DO n = 1, boundary%nNodes - ni = boundary%nodes(n)%obj%n - mesh%K(ni, :) = 0.D0 - mesh%K(ni, ni) = 1.D0 - - END DO - - TYPE IS(boundaryEMDirichletTime) - DO n = 1, boundary%nNodes - ni = boundary%nodes(n)%obj%n - mesh%K(ni, :) = 0.D0 - mesh%K(ni, ni) = 1.D0 - - END DO - - END SELECT - - END DO - - !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', 'readBoundaryEM') - - END IF - END SUBROUTINE readBoundaryEM !Reads the injection of particles from the boundaries diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 7a77ddc..e4ac1d7 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -98,7 +98,7 @@ MODULE moduleMesh1DCart !EDGE FUNCTIONS !Init edge element - SUBROUTINE initEdge1DCart(self, n, p, bt) + SUBROUTINE initEdge1DCart(self, n, p, btPart, btEM) USE moduleSpecies, only: nSpecies USE moduleErrors IMPLICIT NONE @@ -106,7 +106,8 @@ MODULE moduleMesh1DCart CLASS(meshEdge1DCart), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nSpecies) + INTEGER, INTENT(in):: btPart(1:nSpecies) + integer, intent(in):: btEM REAL(8), DIMENSION(1:3):: r1 INTEGER:: s @@ -126,11 +127,14 @@ MODULE moduleMesh1DCart allocate(self%boundariesParticle(1:nSpecies)) !Assign functions to boundary do s = 1, nSpecies - self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj + self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj end do - END SUBROUTINE initEdge1DCart + ! Add nodes to EM boundary + call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes)) + + end subroutine initEdge1DCart !Get nodes from edge PURE FUNCTION getNodes1DCart(self, nNodes) RESULT(n) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index fcc438a..edfece5 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -98,7 +98,7 @@ MODULE moduleMesh1DRad !EDGE FUNCTIONS !Init edge element - SUBROUTINE initEdge1DRad(self, n, p, bt) + SUBROUTINE initEdge1DRad(self, n, p, btPart, btEM) USE moduleSpecies, only:nSpecies USE moduleErrors IMPLICIT NONE @@ -106,7 +106,8 @@ MODULE moduleMesh1DRad CLASS(meshEdge1DRad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nSpecies) + INTEGER, INTENT(in):: btPart(1:nSpecies) + integer, intent(in):: btEM REAL(8), DIMENSION(1:3):: r1 INTEGER:: s @@ -126,11 +127,14 @@ MODULE moduleMesh1DRad allocate(self%boundariesParticle(1:nSpecies)) !Assign functions to boundary do s = 1, nSpecies - self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj + self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj end do - END SUBROUTINE initEdge1DRad + ! Add nodes to EM boundary + call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes)) + + end subroutine initEdge1DRad !Get nodes from edge PURE FUNCTION getNodes1DRad(self, nNodes) RESULT(n) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 5bfb1e4..f928f7f 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -133,7 +133,7 @@ MODULE moduleMesh2DCart !EDGE FUNCTIONS !Init edge element - SUBROUTINE initEdge2DCart(self, n, p, bt) + SUBROUTINE initEdge2DCart(self, n, p, btPart, btEm) USE moduleSpecies, only:nSpecies USE moduleErrors IMPLICIT NONE @@ -141,7 +141,8 @@ MODULE moduleMesh2DCart CLASS(meshEdge2DCart), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nSpecies) + INTEGER, INTENT(in):: btPart(1:nSpecies) + integer, intent(in):: btEM REAL(8), DIMENSION(1:3):: r1, r2 INTEGER:: s @@ -165,10 +166,13 @@ MODULE moduleMesh2DCart allocate(self%boundariesParticle(1:nSpecies)) !Assign functions to boundary do s = 1, nSpecies - self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj + self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj end do + ! Add nodes to EM boundary + call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes)) + END SUBROUTINE initEdge2DCart !Get nodes from edge diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index e606641..137fd68 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -133,7 +133,7 @@ MODULE moduleMesh2DCyl !EDGE FUNCTIONS !Init edge element - SUBROUTINE initEdge2DCyl(self, n, p, bt) + SUBROUTINE initEdge2DCyl(self, n, p, btPart, btEM) USE moduleSpecies USE moduleErrors USE moduleConstParam, ONLY: PI @@ -142,7 +142,8 @@ MODULE moduleMesh2DCyl CLASS(meshEdge2DCyl), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nsPecies) + INTEGER, INTENT(in):: btPart(1:nsPecies) + integer, intent(in):: btEM REAL(8), DIMENSION(1:3):: r1, r2 INTEGER:: s @@ -174,10 +175,13 @@ MODULE moduleMesh2DCyl allocate(self%boundariesParticle(1:nSpecies)) !Assign functions to boundary do s = 1, nSpecies - self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj + self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj end do + ! Add nodes to EM boundary + call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes)) + END SUBROUTINE initEdge2DCyl !Get nodes from edge diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 1655ab8..78f19c7 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -104,7 +104,7 @@ MODULE moduleMesh3DCart !EDGE FUNCTIONS !Init surface element - SUBROUTINE initEdge3DCartTria(self, n, p, bt) + SUBROUTINE initEdge3DCartTria(self, n, p, btPart, btEM) USE moduleSpecies, only: nSpecies USE moduleErrors USE moduleMath @@ -114,7 +114,8 @@ MODULE moduleMesh3DCart CLASS(meshEdge3DCartTria), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nSpecies) + INTEGER, INTENT(in):: btPart(1:nSpecies) + integer, intent(in):: btEM REAL(8), DIMENSION(1:3):: r1, r2, r3 REAL(8), DIMENSION(1:3):: vec1, vec2 INTEGER:: s @@ -145,12 +146,15 @@ MODULE moduleMesh3DCart !Boundary particle linking allocate(self%boundariesParticle(1:nSpecies)) - !Assign functions to boundary + !Assign functions to particle boundary do s = 1, nSpecies - self%boundariesParticle(s)%obj => boundariesParticle(bt(s))%obj + self%boundariesParticle(s)%obj => boundariesParticle(btPart(s))%obj end do + ! Add nodes to EM boundary + call boundariesEMLinking(btEM)%model%addNodes(self%getNodes(self%nNodes)) + END SUBROUTINE initEdge3DCartTria !Get nodes from surface diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index a495ab9..cbf1dd7 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -36,7 +36,8 @@ MODULE moduleMeshInputGmsh2 CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8):: r(1:3) !3 generic coordinates INTEGER, ALLOCATABLE:: p(:) !Array for nodes - INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0, bt = 0 + INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0 + integer:: btPart = 0, btEM = 0 INTEGER:: totalNumElem INTEGER:: numEdges INTEGER:: boundaryType @@ -199,9 +200,10 @@ MODULE moduleMeshInputGmsh2 END SELECT !Associate boundary condition procedure. - bt = physicalSurface_to_id(boundariesParticleLinking,boundaryType) + btPart = physicalSurface_to_id(boundariesParticleLinking,boundaryType) + btEM = physicalSurface_to_id(boundariesEMLinking, boundaryType) - CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex) + CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(btPart)%speciesIndex, btEM) DEALLOCATE(p) END DO diff --git a/src/modules/mesh/inout/text/moduleMeshInputText.f90 b/src/modules/mesh/inout/text/moduleMeshInputText.f90 index 3ada5b9..2f685ea 100644 --- a/src/modules/mesh/inout/text/moduleMeshInputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshInputText.f90 @@ -47,7 +47,7 @@ module moduleMeshInputText integer:: physicalID integer:: n, c integer, allocatable:: p(:) - integer:: bt + integer:: btPart, btEM fileID = 10 @@ -142,8 +142,9 @@ module moduleMeshInputText allocate(p(1)) p(1) = n !Associate boundary condition procedure. - bt = physicalSurface_to_id(boundariesParticleLinking, physicalID) - call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(bt)%speciesIndex) + btPart = physicalSurface_to_id(boundariesParticleLinking, physicalID) + btEM = physicalSurface_to_id(boundariesEMLinking, physicalID) + call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(btPart)%speciesIndex, btEM) deallocate(p) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index c775db8..4b7542d 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -173,7 +173,7 @@ MODULE moduleMeshInputVTU REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates INTEGER:: n, e, c INTEGER, ALLOCATABLE:: p(:) - INTEGER:: bt + INTEGER:: btPart, btEM fileID = 10 @@ -365,10 +365,12 @@ MODULE moduleMeshInputVTU END SELECT !Associate boundary condition procedure. - bt = physicalSurface_to_id(boundariesParticleLinking,entitiesID(n)) + btPart = physicalSurface_to_id(boundariesParticleLinking, entitiesID(n)) + btEM = physicalSurface_to_id(boundariesEMLinking, entitiesID(n)) !Init edge - CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex) + CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(btPart)%speciesIndex, btEM) + DEALLOCATE(p) END DO diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 3f7e139..b78b61c 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -495,7 +495,7 @@ MODULE moduleMesh END TYPE meshParticles interface - pure module subroutine constructGlobalK(self) + module subroutine constructGlobalK(self) class(meshParticles), intent(inout):: self end subroutine constructGlobalK @@ -626,6 +626,12 @@ MODULE moduleMesh end subroutine initBoundaryParticle + module function boundaryParticleName_to_Index(boundaryName) result(bp) + character(:), allocatable:: boundaryName + integer:: bp + + end function boundaryParticleName_to_Index + end interface abstract interface @@ -915,12 +921,13 @@ MODULE moduleMesh type:: boundaryEMLinking integer:: physicalSurface - class(boundaryEMGeneric), pointer:: boundary => null() + class(boundaryEMGeneric), pointer:: model => null() end type boundaryEMLinking INTEGER:: nBoundariesEM TYPE(boundaryEMCont), ALLOCATABLE, target:: boundariesEM(:) + type(boundaryEMLinking), allocatable, dimension(:):: boundariesEMLinking !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index a319bdd..a9b20ae 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -1,5 +1,29 @@ submodule(moduleMesh) boundaryEM CONTAINS + function boundaryEMName_to_Index(boundaryName) result(bp) + use moduleErrors + implicit none + + character(:), allocatable:: boundaryName + integer:: bp + integer:: b + + bp = 0 + do b = 1, nBoundaries + 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 + subroutine addNodes(self, nodes) implicit none @@ -159,7 +183,7 @@ submodule(moduleMesh) boundaryEM CLASS(boundaryEMDirichlet), INTENT(in):: self REAL(8), INTENT(inout):: vectorF(:) - INTEGER:: n, ni + integer:: n DO n = 1, self%nNodes self%nodes(n)%obj%emData%phi = self%potential @@ -178,7 +202,7 @@ submodule(moduleMesh) boundaryEM CLASS(boundaryEMDirichletTime), INTENT(in):: self REAL(8), INTENT(inout):: vectorF(:) REAL(8):: timeFactor - INTEGER:: n, ni + integer:: n timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) @@ -190,7 +214,7 @@ submodule(moduleMesh) boundaryEM END SUBROUTINE applyDirichletTime - function physicalSurface_to_boundaryEMId(boundaryArray,physicalSurface) result(b) + function physicalSurface_to_boundaryEMId(boundaryArray, physicalSurface) result(b) implicit none type(boundaryEMLinking):: boundaryArray(:) @@ -198,9 +222,10 @@ submodule(moduleMesh) boundaryEM integer:: b integer:: bt + b = 0 do bt=1, nPhysicalSurfaces if (boundaryArray(bt)%physicalSurface == physicalSurface) then - b = boundaryArray(bt)%boundary%n + b = boundaryArray(bt)%model%n exit diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 8ed7ac8..fbc55eb 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -1,7 +1,7 @@ !moduleMeshBoundary: Boundary functions for the mesh edges submodule(moduleMesh) boundaryParticle contains - module function boundaryParticleName_to_Index(boundaryName) result(bp) + function boundaryParticleName_to_Index(boundaryName) result(bp) use moduleErrors implicit none @@ -19,7 +19,7 @@ submodule(moduleMesh) boundaryParticle end do if (bp == 0) then - call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryName2Index') + call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryParticleName_to_Index') end if @@ -503,6 +503,7 @@ submodule(moduleMesh) boundaryParticle integer:: b integer:: bt + b = 0 do bt=1, nPhysicalSurfaces if (boundaryArray(bt)%physicalSurface == physicalSurface) then b = bt diff --git a/src/modules/mesh/moduleMesh@elements.f90 b/src/modules/mesh/moduleMesh@elements.f90 index 2291919..81c03e3 100644 --- a/src/modules/mesh/moduleMesh@elements.f90 +++ b/src/modules/mesh/moduleMesh@elements.f90 @@ -1,7 +1,7 @@ submodule(moduleMesh) elements CONTAINS !Reset the output of node - PURE module SUBROUTINE resetOutput(self) + PURE SUBROUTINE resetOutput(self) USE moduleSpecies USE moduleOutput IMPLICIT NONE @@ -40,40 +40,73 @@ submodule(moduleMesh) elements end function meshNodePointer_equal_type_int !Constructs the global K matrix - PURE module SUBROUTINE constructGlobalK(self) + SUBROUTINE constructGlobalK(self) IMPLICIT NONE CLASS(meshParticles), INTENT(inout):: self INTEGER:: c - INTEGER, ALLOCATABLE:: n(:) + INTEGER, ALLOCATABLE:: nodes(:) REAL(8), ALLOCATABLE:: localK(:,:) INTEGER:: i, j + integer:: n, b, ni + INTEGER:: info + EXTERNAL:: dgetrf DO c = 1, self%numCells associate(nNodes => self%cells(c)%obj%nNodes) - ALLOCATE(n(1:nNodes)) + ALLOCATE(nodes(1:nNodes)) ALLOCATE(localK(1:nNodes, 1:nNodes)) - n = self%cells(c)%obj%getNodes(nNodes) + nodes = self%cells(c)%obj%getNodes(nNodes) localK = self%cells(c)%obj%elemK(nNodes) DO i = 1, nNodes DO j = 1, nNodes - self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) + self%K(nodes(i), nodes(j)) = self%K(nodes(i), nodes(j)) + localK(i, j) END DO END DO - DEALLOCATE(n, localK) + DEALLOCATE(nodes, localK) end associate END DO + ! Modify K matrix due to EM boundary conditions + DO b = 1, nBoundariesEM + SELECT TYPE(boundary => boundariesEM(b)%obj) + TYPE IS(boundaryEMDirichlet) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + self%K(ni, :) = 0.D0 + self%K(ni, ni) = 1.D0 + + END DO + + TYPE IS(boundaryEMDirichletTime) + 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 + + !Compute the PLU factorization of K once boundary conditions have been read + CALL dgetrf(self%numNodes, self%numNodes, self%K, self%numNodes, self%IPIV, info) + IF (info /= 0) THEN + CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM') + + END IF + END SUBROUTINE constructGlobalK ! Gather the value of valNodes at position Xi of an edge - pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f) + pure function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f) implicit none class(meshEdge), intent(in):: self @@ -89,7 +122,7 @@ submodule(moduleMesh) elements end function gatherF_edge_scalar !Gather the value of valNodes (scalar) at position Xi - PURE module FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f) + PURE FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self @@ -105,7 +138,7 @@ submodule(moduleMesh) elements END FUNCTION gatherF_cell_scalar !Gather the value of valNodes (array) at position Xi - PURE module FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f) + PURE FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self @@ -121,7 +154,7 @@ submodule(moduleMesh) elements END FUNCTION gatherF_cell_array !Gather the spatial derivative of valNodes (scalar) at position Xi - PURE module FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df) + PURE FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self @@ -146,7 +179,7 @@ submodule(moduleMesh) elements END FUNCTION gatherDF_cell_scalar !Scatters particle properties into cell nodes - module SUBROUTINE scatter(self, nNodes, part) + SUBROUTINE scatter(self, nNodes, part) USE moduleMath USE moduleSpecies USE OMP_LIB @@ -184,7 +217,7 @@ submodule(moduleMesh) elements END SUBROUTINE scatter !Find next cell for particle - RECURSIVE module SUBROUTINE findCell(self, part, oldCell) + RECURSIVE SUBROUTINE findCell(self, part, oldCell) USE moduleSpecies USE moduleErrors USE OMP_LIB @@ -248,7 +281,7 @@ submodule(moduleMesh) elements END SUBROUTINE findCell !If Coll and Particle are the same, simply copy the part%cell into part%cellColl - module SUBROUTINE findCellSameMesh(part) + SUBROUTINE findCellSameMesh(part) USE moduleSpecies IMPLICIT NONE @@ -261,7 +294,7 @@ submodule(moduleMesh) elements !TODO: try to combine this with the findCell for a regular mesh !Find the volume in which particle reside in the mesh for collisions !No boundary interaction taken into account - module SUBROUTINE findCellCollMesh(part) + SUBROUTINE findCellCollMesh(part) USE moduleSpecies IMPLICIT NONE @@ -311,7 +344,7 @@ submodule(moduleMesh) elements !Returns index of volume associated to a position (if any) !If no voulme is found, returns 0 !WARNING: This function is slow and should only be used in initialization phase - module FUNCTION findCellBrute(self, r) RESULT(nVol) + FUNCTION findCellBrute(self, r) RESULT(nVol) USE moduleSpecies IMPLICIT NONE @@ -337,7 +370,7 @@ submodule(moduleMesh) elements END FUNCTION findCellBrute !Computes collisions in element - module SUBROUTINE doCollisions(self) + SUBROUTINE doCollisions(self) USE moduleCollisions USE moduleSpecies USE moduleList @@ -501,7 +534,7 @@ submodule(moduleMesh) elements END SUBROUTINE doCollisions - module SUBROUTINE doCoulomb(self) + SUBROUTINE doCoulomb(self) USE moduleCoulomb USE moduleRandom USE moduleOutput