From fbbb0d5d135caf3f83b77c4bd1e9ce9cce30b893 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 5 Feb 2026 15:30:50 +0100 Subject: [PATCH] fPsi functions for edges I need to make a common module for mesh, many functions for elements are shared. Also, try to reduce the 'select type' statements, but I don't know enough Fortran for it. --- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 25 ++++++++++++----- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 25 ++++++++++++----- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 24 ++++++++++++++--- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 28 +++++++++++++++----- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 11 ++++---- src/modules/mesh/moduleMesh.f90 | 20 ++++++++++++++ src/modules/mesh/moduleMeshBoundary.f90 | 9 ++++++- 7 files changed, 114 insertions(+), 28 deletions(-) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index a304c5e..dd86109 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -27,10 +27,11 @@ MODULE moduleMesh1DCart CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS !meshEdge DEFERRED PROCEDURES - PROCEDURE, PASS:: init => initEdge1DCart - PROCEDURE, PASS:: getNodes => getNodes1DCart - PROCEDURE, PASS:: intersection => intersection1DCart - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: init => initEdge1DCart + PROCEDURE, PASS:: getNodes => getNodes1DCart + PROCEDURE, PASS:: intersection => intersection1DCart + PROCEDURE, PASS:: randPos => randPosEdge + procedure, nopass:: fPsi => fPsiPoint END TYPE meshEdge1DCart @@ -172,7 +173,19 @@ MODULE moduleMesh1DCart END FUNCTION randPosEdge - !VOLUME FUNCTIONS + !Compute element functions at point Xi + pure function fPsiPoint(Xi, nNodes) RESULT(fPsi) + implicit none + + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8):: fPsi(1:nNodes) + + fPsi = (/ 1.D0 /) + + end function fPsiPoint + + !CELL FUNCTIONS !SEGMENT FUNCTIONS !Init element SUBROUTINE initCell1DCartSegm(self, n, p, nodes) @@ -446,7 +459,7 @@ MODULE moduleMesh1DCart END SUBROUTINE volumeSegm - !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS + !COMMON FUNCTIONS FOR 1D CELL ELEMENTS !Compute element Jacobian determinant PURE FUNCTION detJ1DCart(pDer) RESULT(dJ) IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index a46c6a2..3e36711 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -27,10 +27,11 @@ MODULE moduleMesh1DRad CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS !meshEdge DEFERRED PROCEDURES - PROCEDURE, PASS:: init => initEdge1DRad - PROCEDURE, PASS:: getNodes => getNodes1DRad - PROCEDURE, PASS:: intersection => intersection1DRad - PROCEDURE, PASS:: randPos => randPos1DRad + PROCEDURE, PASS:: init => initEdge1DRad + PROCEDURE, PASS:: getNodes => getNodes1DRad + PROCEDURE, PASS:: intersection => intersection1DRad + PROCEDURE, PASS:: randPos => randPos1DRad + procedure, nopass:: fPsi => fPsiPoint END TYPE meshEdge1DRad @@ -172,7 +173,19 @@ MODULE moduleMesh1DRad END FUNCTION randPos1DRad - !VOLUME FUNCTIONS + !Compute element functions at point Xi + pure function fPsiPoint(Xi, nNodes) RESULT(fPsi) + implicit none + + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8):: fPsi(1:nNodes) + + fPsi = (/ 1.D0 /) + + end function fPsiPoint + + !CELL FUNCTIONS !SEGMENT FUNCTIONS !Init element SUBROUTINE initCell1DRadSegm(self, n, p, nodes) @@ -461,7 +474,7 @@ MODULE moduleMesh1DRad END SUBROUTINE volumeSegm - !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS + !COMMON FUNCTIONS FOR 1D CELL ELEMENTS !Compute element Jacobian determinant PURE FUNCTION detJ1DRad(pDer) RESULT(dJ) IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index db1ef4f..092c0d6 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -32,10 +32,11 @@ MODULE moduleMesh2DCart CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS !meshEdge DEFERRED PROCEDURES - PROCEDURE, PASS:: init => initEdge2DCart - PROCEDURE, PASS:: getNodes => getNodes2DCart - PROCEDURE, PASS:: intersection => intersection2DCartEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: init => initEdge2DCart + PROCEDURE, PASS:: getNodes => getNodes2DCart + PROCEDURE, PASS:: intersection => intersection2DCartEdge + PROCEDURE, PASS:: randPos => randPosEdge + procedure, nopass:: fPsi => fPsiSegm END TYPE meshEdge2DCart @@ -234,6 +235,21 @@ MODULE moduleMesh2DCart END FUNCTION randPosEdge + !Compute element functions at point Xi + pure function fPsiSegm(Xi, nNodes) RESULT(fPsi) + implicit none + + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8):: fPsi(1:nNodes) + + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) + + fPsi = fPsi * 0.50D0 + + end function fPsiSegm + !VOLUME FUNCTIONS !QUAD FUNCTIONS !Init element diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 1f5f33c..109512e 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -32,10 +32,11 @@ MODULE moduleMesh2DCyl CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS !meshEdge DEFERRED PROCEDURES - PROCEDURE, PASS:: init => initEdge2DCyl - PROCEDURE, PASS:: getNodes => getNodes2DCyl - PROCEDURE, PASS:: intersection => intersection2DCylEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: init => initEdge2DCyl + PROCEDURE, PASS:: getNodes => getNodes2DCyl + PROCEDURE, PASS:: intersection => intersection2DCylEdge + PROCEDURE, PASS:: randPos => randPosEdge + procedure, nopass:: fPsi => fPsiSegm END TYPE meshEdge2DCyl @@ -243,7 +244,22 @@ MODULE moduleMesh2DCyl END FUNCTION randPosEdge - !VOLUME FUNCTIONS + !Compute element functions at point Xi + pure function fPsiSegm(Xi, nNodes) RESULT(fPsi) + implicit none + + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8):: fPsi(1:nNodes) + + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) + + fPsi = fPsi * 0.50D0 + + end function fPsiSegm + + !CELL FUNCTIONS !QUAD FUNCTIONS !Init element SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) @@ -929,7 +945,7 @@ MODULE moduleMesh2DCyl END SUBROUTINE volumeTria - !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS + !COMMON FUNCTIONS FOR CYLINDRICAL CELL ELEMENTS !Compute element Jacobian determinant PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ) IMPLICIT NONE diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 8170df7..0a4ba79 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -30,7 +30,7 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: intersection => intersection3DCartTria PROCEDURE, PASS:: randPos => randPosEdgeTria !PARTICULAR PROCEDURES - PROCEDURE, NOPASS, PRIVATE:: fPsi => fPsiEdgeTria + PROCEDURE, NOPASS:: fPsi => fPsiTria END TYPE meshEdge3DCartTria @@ -204,7 +204,7 @@ MODULE moduleMesh3DCart Xi(2) = random( 0.D0, 1.D0 - Xi(1)) Xi(3) = 0.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 3) r = (/DOT_PRODUCT(fPsi, self%x), & DOT_PRODUCT(fPsi, self%y), & DOT_PRODUCT(fPsi, self%z)/) @@ -212,17 +212,18 @@ MODULE moduleMesh3DCart END FUNCTION randPosEdgeTria !Shape functions for triangular surface - PURE FUNCTION fPsiEdgeTria(Xi) RESULT(fPsi) + PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: fPsi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) fPsi(1) = 1.D0 - Xi(1) - Xi(2) fPsi(2) = Xi(1) fPsi(3) = Xi(2) - END FUNCTION fPsiEdgeTria + END FUNCTION fPsiTria !VOLUME FUNCTIONS !TETRA FUNCTIONS diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index c4ac680..8a88e78 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -97,6 +97,10 @@ MODULE moduleMesh PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + !Gather value and spatial derivative on the nodes at position Xi + PROCEDURE, PASS, PRIVATE:: gatherF_edge_scalar + GENERIC:: gatherF => gatherF_edge_scalar END TYPE meshEdge @@ -561,6 +565,22 @@ MODULE moduleMesh END FUNCTION gatherF_cell_scalar + ! Gather the value of valNodes at position Xi of an edge + pure function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f) + implicit none + + class(meshEdge), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes) + real(8):: f + real(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = dot_product(fPsi, valNodes) + + end function gatherF_edge_scalar + !Gather the value of valNodes (array) at position Xi PURE FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f) IMPLICIT NONE diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index e81e796..6b8be87 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -223,12 +223,19 @@ MODULE moduleMeshBoundary class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part + real(8), allocatable:: density(:) + integer:: nNodes + integer, allocatable:: nodes(:) + select type(bound => edge%boundary%bTypes(part%species%n)%obj) type is(boundaryQuasiNeutrality) + nNodes = edge%nNodes + nodes = edge%getNodes(nNodes) + bound%alpha = 1.0d-1 - if (random() < bound%alpha) then + if (random() <= bound%alpha) then call reflection(edge, part) else