diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 117b77d..7b685ab 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -29,6 +29,7 @@ MODULE moduleMesh1DCart PROCEDURE, PASS:: intersection => intersection1DCart PROCEDURE, PASS:: randPos => randPosEdge procedure, pass:: center => centerEdgePoint + procedure, nopass:: centerXi => centerXiPoint procedure, nopass:: fPsi => fPsiPoint END TYPE meshEdge1DCart diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 3f6df2b..99aef3e 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -29,6 +29,7 @@ MODULE moduleMesh1DRad PROCEDURE, PASS:: intersection => intersection1DRad PROCEDURE, PASS:: randPos => randPos1DRad procedure, pass:: center => centerEdgePoint + procedure, nopass:: centerXi => centerXiPoint procedure, nopass:: fPsi => fPsiPoint END TYPE meshEdge1DRad diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 301cc6b..46fe1f1 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -29,6 +29,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: intersection => intersection2DCartEdge PROCEDURE, PASS:: randPos => randPosEdge procedure, pass:: center => centerEdgeSegm + procedure, nopass:: centerXi => centerXiSegm procedure, nopass:: fPsi => fPsiSegm END TYPE meshEdge2DCart diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 0610b87..bc30e0e 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -29,6 +29,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: intersection => intersection2DCylEdge PROCEDURE, PASS:: randPos => randPosEdge procedure, pass:: center => centerEdgeSegm + procedure, nopass:: centerXi => centerXiSegm procedure, nopass:: fPsi => fPsiSegm END TYPE meshEdge2DCyl diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index a68d1f7..03f7f18 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -31,7 +31,8 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: randPos => randPosEdgeTria procedure, pass:: center => centerEdgeTria !PARTICULAR PROCEDURES - PROCEDURE, NOPASS:: fPsi => fPsiTria + procedure, nopass:: centerXi => centerXiTria + PROCEDURE, NOPASS:: fPsi => fPsiTria END TYPE meshEdge3DCartTria diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index d1a6920..5c050e9 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -119,6 +119,7 @@ MODULE moduleMesh PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos procedure(centerEdge_interface), deferred, pass:: center + procedure(centerXiEdge_interface), deferred, nopass:: centerXi PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi !Gather value and spatial derivative on the nodes at position Xi PROCEDURE, PASS, PRIVATE:: gatherF_edge_scalar @@ -185,6 +186,12 @@ MODULE moduleMesh END FUNCTION centerEdge_interface + ! Returns the center point of an edge in natural coordinates + FUNCTION centerXiEdge_interface() RESULT(Xi) + REAL(8):: Xi(1:3) + + END FUNCTION centerXiEdge_interface + END INTERFACE !Containers for edges in the mesh diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index a940dab..2e2a743 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -346,7 +346,6 @@ submodule(moduleMesh) boundaryParticle class(boundaryParticleGeneric), intent(inout):: self integer, intent(in):: fileID - integer:: e write(fileID, '(A)') self%name select type(self) @@ -432,17 +431,25 @@ submodule(moduleMesh) boundaryParticle end do if (s == self%s_incident) then - density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) + density_incident = edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes) else - density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) + density_rest = density_rest + edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes) end if end do ! Correction for this time step - alpha = 1.d0 - density_incident/density_rest + if (density_rest > 1.0d-10) then + ! If there is a rest population, correct + alpha = 1.d0 - density_incident/density_rest + + else + ! If not, alpha is assumed unchaged + alpha = 0.d0 + + end if ! Apply correction with a factor of 0.1 to avoid fast changes self%alpha(edge%n) = self%alpha(edge%n) + 1.0d-2 * alpha @@ -472,7 +479,7 @@ submodule(moduleMesh) boundaryParticle type is(boundaryQuasiNeutrality) write(fileID, '(A,",",A)') '"Edge id"', '"alpha"' do e = 1, size(self%edges) - write(fileID, '('//fmtInt//',",",'//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) + write(fileID, '('//fmtColInt//','//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n) end do diff --git a/src/modules/mesh/moduleMeshCommon.f90 b/src/modules/mesh/moduleMeshCommon.f90 index e4413aa..98b8a82 100644 --- a/src/modules/mesh/moduleMeshCommon.f90 +++ b/src/modules/mesh/moduleMeshCommon.f90 @@ -15,14 +15,69 @@ module moduleMeshCommon real(8), parameter:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) ! Center point in natural coordinates + ! Point + real(8), parameter:: cenPoint(1:3) = (/ 0.0d0, 0.0d0, 0.0d0 /) ! Segment - real(8), parameter:: cenSeg(1:3) = (/ 0.5d0, 0.0d0, 0.0d0 /) + real(8), parameter:: cenSeg(1:3) = (/ 0.0d0, 0.0d0, 0.0d0 /) ! Tria - real(8), parameter:: cenTria(1:3) = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 0.0d0 /) + real(8), parameter:: cenTria(1:3) = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 0.0d0 /) ! Quad - real(8), parameter:: cenQuad(1:3) = (/ 0.0d0, 0.0d0, 0.0d0 /) + real(8), parameter:: cenQuad(1:3) = (/ 0.0d0, 0.0d0, 0.0d0 /) + ! Tetra + real(8), parameter:: cenTetra(1:3) = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 1.d0/3.d0 /) contains + ! RETURN CENTER IN NATURAL COORDINATES + ! Point + pure function centerXiPoint() result(Xi) + implicit none + + real(8):: Xi(1:3) + + Xi = cenPoint + + end function centerXiPoint + + ! Segment + pure function centerXiSegm() result(Xi) + implicit none + + real(8):: Xi(1:3) + + Xi = cenSeg + + end function centerXiSegm + + ! Tria + pure function centerXiTria() result(Xi) + implicit none + + real(8):: Xi(1:3) + + Xi = cenTria + + end function centerXiTria + + ! Quad + pure function centerXiQuad() result(Xi) + implicit none + + real(8):: Xi(1:3) + + Xi = cenQuad + + end function centerXiQuad + + ! Quad + pure function centerXiTetra() result(Xi) + implicit none + + real(8):: Xi(1:3) + + Xi = cenTetra + + end function centerXiTetra + ! ELEMENT FUNCTIONS ! Point pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)