From db6b0a2c03128352af55a0abf9b102304a57125d Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 22 Mar 2021 12:39:34 +0100 Subject: [PATCH] Fixed an issue with reflection of particles in all geometries and also assigning the normal vector in 2D and 3D. 3D Cartesian geometry is working properly, although it needs testing. Still issue with ionization boundary. --- src/makefile | 2 +- src/modules/makefile | 6 +- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 7 +- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 56 +++-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 25 ++- .../mesh/2DCart/moduleMesh2DCartRead.f90 | 70 +++--- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 26 ++- .../mesh/2DCyl/moduleMesh2DCylRead.f90 | 90 +++++--- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 41 ++-- .../mesh/3DCart/moduleMesh3DCartRead.f90 | 211 +++++++++--------- src/modules/mesh/moduleMesh.f90 | 4 +- src/modules/mesh/moduleMeshBoundary.f90 | 20 +- src/modules/moduleInject.f90 | 3 +- src/modules/moduleMath.f90 | 42 ++++ src/modules/moduleOutput.f90 | 11 +- 15 files changed, 349 insertions(+), 265 deletions(-) create mode 100644 src/modules/moduleMath.f90 diff --git a/src/makefile b/src/makefile index 7bb277a..b398134 100644 --- a/src/makefile +++ b/src/makefile @@ -3,7 +3,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ - $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ + $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ $(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o \ $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \ $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \ diff --git a/src/modules/makefile b/src/modules/makefile index 28fc4ec..cea6d73 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -2,11 +2,11 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \ moduleOutput.o moduleInput.o moduleSolver.o \ moduleCollisions.o moduleTable.o moduleParallel.o \ - moduleEM.o moduleRandom.o + moduleEM.o moduleRandom.o moduleMath.o all: $(OBJS) mesh.o -mesh.o: moduleCollisions.o moduleBoundary.o +mesh.o: moduleCollisions.o moduleBoundary.o moduleMath.o $(MAKE) -C mesh all moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90 @@ -21,7 +21,7 @@ moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90 +moduleOutput.o: moduleMath.o moduleSpecies.o moduleRefParam.o moduleOutput.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ moduleRandom.o: moduleConstParam.o moduleRandom.f90 diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 366fa84..020b971 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -148,7 +148,6 @@ MODULE moduleMesh1DCart self%x = r1(1) self%normal = (/ 1.D0, 0.D0, 0.D0 /) - self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -176,11 +175,11 @@ MODULE moduleMesh1DCart END FUNCTION getNodes1DCart - PURE FUNCTION intersection1DCart(self, r0, v0) RESULT(r) + PURE FUNCTION intersection1DCart(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r r = (/ self%x, 0.D0, 0.D0 /) @@ -383,7 +382,7 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm SUBROUTINE scatterSegm(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index d0f8311..08609cd 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -7,6 +7,9 @@ MODULE moduleMesh1DRad USE moduleMeshBoundary IMPLICIT NONE + REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) + REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /) + TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad !Element coordinates REAL(8):: r = 0.D0 @@ -174,11 +177,11 @@ MODULE moduleMesh1DRad END FUNCTION getNodes1DRad - PURE FUNCTION intersection1DRad(self, r0, v0) RESULT(r) + PURE FUNCTION intersection1DRad(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r r = (/ self%r, 0.D0, 0.D0 /) @@ -317,19 +320,26 @@ MODULE moduleMesh1DRad REAL(8):: ke(1:2,1:2) REAL(8):: Xii(1:3) REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ + REAL(8):: invJ(1), detJ REAL(8):: r, fPsi(1:2) + INTEGER:: l ke = 0.D0 Xii = 0.D0 - dPsi = self%dPsi(Xii) - invJ = self%invJac(Xii, dPsi) - fPsi = self%fPsi(Xii) - r = DOT_PRODUCT(fPsi, self%r) - ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) - ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) - ke = 2.D0*ke*invJ - ke = ke*r*PI2 + DO l = 1, 3 + xii(1) = corSeg(l) + dPsi = self%dPsi(Xii) + detJ = self%detJac(Xii, dPsi) + invJ = self%invJac(Xii, dPsi) + fPsi = self%fPsi(Xii) + r = DOT_PRODUCT(fPsi, self%r) + ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r*wSeg(l)/detJ + + END DO + + ke = ke*PI2 END FUNCTION elemKRad @@ -341,17 +351,23 @@ MODULE moduleMesh1DRad REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:2) - REAL(8):: r - REAL(8):: detJ + REAL(8):: detJ, f, r REAL(8):: Xii(1:3) + INTEGER:: l - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - r = DOT_PRODUCT(fPsi,self%r) ALLOCATE(localF(1:2)) - localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ - localF = localF*r*PI2 + localF = 0.D0 + Xii = 0.D0 + + DO l = 1, 3 + xii(1) = corSeg(l) + detJ = self%detJac(Xii) + fPsi = self%fPsi(Xii) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) + localF = localF + f*fPsi*r*wSeg(l)*detJ + + END DO END FUNCTION elemFRad @@ -377,7 +393,7 @@ MODULE moduleMesh1DRad END FUNCTION insideRad SUBROUTINE scatterRad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 40c218c..10b12fe 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -190,9 +190,9 @@ MODULE moduleMesh2DCart self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) !Normal vector - self%normal = (/ self%y(2)-self%y(1), & - self%x(2)-self%x(1), & - 0.D0 /) + self%normal = (/ -(self%y(2)-self%y(1)), & + self%x(2)-self%x(1) , & + 0.D0 /) self%normal = self%normal/NORM2(self%normal) !Boundary index @@ -243,20 +243,21 @@ MODULE moduleMesh2DCart END FUNCTION getNodes2DCart - PURE FUNCTION intersection2DCartEdge(self, r0, v0) RESULT(r) + PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%x(1), self%y(1), 0.D0 /) + edge0 = (/self%x(1), self%y(1), 0.D0 /) + edgeV = (/self%x(2), self%y(2), 0.D0 /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection2DCartEdge @@ -496,7 +497,7 @@ MODULE moduleMesh2DCart !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -855,7 +856,7 @@ MODULE moduleMesh2DCart !Scatter properties of particles into element SUBROUTINE scatterTria(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 index b0f2ee3..96664df 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -437,15 +437,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -456,15 +456,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -475,15 +475,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n4%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -494,15 +494,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n4%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e4 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -521,16 +521,16 @@ MODULE moduleMesh2DCartRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + END IF END IF @@ -540,15 +540,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -559,15 +559,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index e0f3a27..d48af83 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -191,9 +191,9 @@ MODULE moduleMesh2DCyl self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) !Normal vector - self%normal = (/ self%r(2)-self%r(1), & - self%z(2)-self%z(1), & - 0.D0 /) + self%normal = (/ -(self%r(2)-self%r(1)), & + self%z(2)-self%z(1) , & + 0.D0 /) self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -221,20 +221,22 @@ MODULE moduleMesh2DCyl END FUNCTION getNodes2DCyl - PURE FUNCTION intersection2DCylEdge(self, r0, v0) RESULT(r) + PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r) + USE moduleMath IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%z(1), self%r(1), 0.D0 /) + edge0 = (/self%z(1), self%r(1), 0.D0 /) + edgeV = (/self%z(2), self%r(2), 0.D0 /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection2DCylEdge @@ -507,7 +509,7 @@ MODULE moduleMesh2DCyl !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -874,7 +876,7 @@ MODULE moduleMesh2DCyl !Scatter properties of particles into element SUBROUTINE scatterTria(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 index 379a0fb..cbe3faf 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 @@ -217,6 +217,29 @@ MODULE moduleMesh2DCylRead END SUBROUTINE connectedVolEdge + PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) + IMPLICIT NONE + + INTEGER, DIMENSION(1:2), INTENT(in):: nodesA, nodesB + LOGICAL:: coincident + INTEGER:: i + + coincident = .FALSE. + DO i = 1, 2 + IF (ANY(nodesA(i) == nodesB)) THEN + coincident = .TRUE. + + ELSE + coincident = .FALSE. + EXIT + + END IF + + END DO + + END FUNCTION coincidentNodes + + SUBROUTINE connectedQuadQuad(elemA, elemB) IMPLICIT NONE @@ -437,15 +460,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -456,15 +479,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -475,15 +498,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n4%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -494,12 +517,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n4%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e4 => elemB - elemB%e2 => elemA + elemB%e1 => elemA ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -518,15 +544,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -537,15 +563,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -556,15 +582,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 315cb42..ec7f2b1 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -134,6 +134,7 @@ MODULE moduleMesh3DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleMath IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(out):: self @@ -142,6 +143,7 @@ MODULE moduleMesh3DCart INTEGER, INTENT(in):: bt INTEGER, INTENT(in):: physicalSurface REAL(8), DIMENSION(1:3):: r1, r2, r3 + REAL(8), DIMENSION(1:3):: vec1, vec2 INTEGER:: s self%n = n @@ -156,10 +158,14 @@ MODULE moduleMesh3DCart self%y = (/r1(2), r2(2), r3(2)/) self%z = (/r1(3), r2(3), r3(3)/) !Normal vector - self%normal = (/ (self%y(2)-self%y(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)), & - (self%x(2)-self%x(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%x(3)-self%x(1)), & - (self%x(2)-self%x(1))*(self%y(3)-self%y(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)) /) - self%normal = self%normal/NORM2(self%normal) + vec1 = (/ self%x(2) - self%x(1), & + self%y(2) - self%y(1), & + self%z(2) - self%z(1) /) + vec2 = (/ self%x(3) - self%x(1), & + self%y(3) - self%y(1), & + self%z(3) - self%z(1) /) + self%normal = crossProduct(vec1, vec2) + self%normal = normalize(self%normal) !Boundary index self%boundary => boundary(bt) @@ -187,20 +193,21 @@ MODULE moduleMesh3DCart END FUNCTION getNodes3DCartTria - PURE FUNCTION intersection3DCartTria(self, r0, v0) RESULT(r) + PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), INTENT(in):: r0(1:3) REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%x(1), self%y(1), self%z(1) /) + edge0 = (/self%x(1), self%y(1), self%z(1) /) + edgeV = (/self%x(2), self%y(2), self%z(2) /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection3DCartTria @@ -476,7 +483,7 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra SUBROUTINE scatterTetra(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -577,7 +584,7 @@ MODULE moduleMesh3DCart INTEGER:: nextInt !TODO: Review when connectivity - xiArray = (/ xi(3), xi(2), 1.D0 - xi(1) - xi(2) - xi(3), xi(1) /) + xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /) nextInt = MINLOC(xiArray, 1) NULLIFY(nextElement) SELECT CASE(nextInt) @@ -647,9 +654,11 @@ MODULE moduleMesh3DCart invJ(2,2) = (dx(1)*dz(3) - dx(3)*dz(1)) invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1)) - invJ(3,1) = -(dx(2)*dy(3) - dx(3)*dy(2)) - invJ(3,2) = (dx(1)*dy(3) - dx(3)*dy(1)) - invJ(3,3) = -(dx(1)*dy(2) - dx(2)*dy(1)) + invJ(3,1) = (dx(2)*dy(3) - dx(3)*dy(2)) + invJ(3,2) = -(dx(1)*dy(3) - dx(3)*dy(1)) + invJ(3,3) = (dx(1)*dy(2) - dx(2)*dy(1)) + + invJ = TRANSPOSE(invJ) END FUNCTION invJ3DCart diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 index c6b021a..01aef8e 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 @@ -269,25 +269,25 @@ MODULE moduleMesh3DCartRead !Check surface 2 IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e2 => elemB elemB%e1 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e2 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e3 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB @@ -299,25 +299,25 @@ MODULE moduleMesh3DCartRead !Check surface 3 IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e3 => elemB elemB%e1 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e2 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e3 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB @@ -328,7 +328,7 @@ MODULE moduleMesh3DCartRead END IF !Check surface 4 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (.NOT. ASSOCIATED(elemA%e4)) THEN IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN @@ -360,77 +360,78 @@ MODULE moduleMesh3DCartRead END SUBROUTINE connectedTetraTetra SUBROUTINE connectedTetraEdge(elemA, elemB) + USE moduleMath IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB + INTEGER:: nodesEdge(1:3) + REAL(8), DIMENSION(1:3):: vec1, vec2 + REAL(8):: normVol(1:3) + + nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /) !Check surface 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), & + nodesEdge)) THEN - elemA%e1 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - ELSEIF ((elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) .OR. & - (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n)) THEN + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e1 => elemB - elemB%e2 => elemA + elemA%e1 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF - END IF !Check surface 2 IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + nodesEdge)) THEN - elemA%e2 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(3) - elemA%x(2), & + elemA%y(3) - elemA%y(2), & + elemA%z(3) - elemA%z(2) /) + vec2 = (/ elemA%x(4) - elemA%x(2), & + elemA%y(4) - elemA%y(2), & + elemA%z(4) - elemA%z(2) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - ELSEIF ((elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) .OR. & - (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n)) THEN + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e2 => elemB - elemB%e2 => elemA + elemA%e2 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF @@ -438,66 +439,66 @@ MODULE moduleMesh3DCartRead !Check surface 3 IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), & + nodesEdge)) THEN - elemA%e3 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - ELSEIF ((elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) .OR. & - (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n)) THEN + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e3 => elemB - elemB%e2 => elemA + elemA%e3 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF END IF !Check surface 4 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), & + nodesEdge)) THEN - elemA%e4 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - ELSEIF ((elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) .OR. & - (elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n)) THEN + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e4 => elemB - elemB%e2 => elemA + elemA%e4 => elemB + elemB%e1 => elemA + + + ELSE + + elemA%e4 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index f5f7b47..137b256 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -97,10 +97,10 @@ MODULE moduleMesh END FUNCTION getNodesEdge_interface !Returns the intersecction between an edge and a line defined by point r0 and vector v0 - PURE FUNCTION intersectionEdge_interface(self, r0, v0) RESULT(r) + PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self - REAL(8), INTENT(in), DIMENSION(1:3):: r0, v0 + REAL(8), INTENT(in), DIMENSION(1:3):: r0 REAL(8):: r(1:3) END FUNCTION intersectionEdge_interface diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 2252b84..7f5fb88 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -13,30 +13,26 @@ MODULE moduleMeshBoundary !rp = intersection between particle and edge !rpp = final position of particle !vpp = final velocity of particle - REAL(8), DIMENSION(1:3):: rp, rpp, vpp + REAL(8), DIMENSION(1:3):: rp, vpp + REAL(8):: tI REAL(8):: taup !time step for reflecting process !Reflect particle velocity vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal + part%v = vpp - !Computes the intersection between particle and surface - rp = edge%intersection(part%r, part%v) + rp = edge%intersection(part%r) - !Computes the reflection time step - taup = NORM2(part%r - rp)*tau(part%sp) + part%r = 2.D0*(rp - part%r) + part%r - !New position of particle - rpp = rp + vpp*taup - - !assign new parameters to particle - part%r = rpp - part%v = vpp + !particle is assumed to be inside part%n_in = .TRUE. END SUBROUTINE reflection !Absoption in a surface SUBROUTINE absorption(edge, part) + USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -45,7 +41,7 @@ MODULE moduleMeshBoundary REAL(8):: rpp(1:3) !Position of particle projected to the edge REAL(8):: d !Distance from particle to edge - rpp = edge%intersection(part%r, part%v) + rpp = edge%intersection(part%r) d = NORM2(rpp - part%r) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index df9df48..9d50abc 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -200,6 +200,7 @@ MODULE moduleInject USE moduleSolver USE moduleMesh USE moduleRandom + USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self @@ -241,7 +242,7 @@ MODULE moduleInject partInj(n)%vol = randomEdge%e2%n ELSE - PRINT *, "ERROR NO VOL ASSOCIATED TO EDGE" + CALL criticalError("No Volume associated to edge", 'addParticles') END IF diff --git a/src/modules/moduleMath.f90 b/src/modules/moduleMath.f90 new file mode 100644 index 0000000..213afa5 --- /dev/null +++ b/src/modules/moduleMath.f90 @@ -0,0 +1,42 @@ +MODULE moduleMath + IMPLICIT NONE + + CONTAINS + !Outer product of two tensors + PURE FUNCTION outerProduct(a,b) RESULT(s) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a, b + REAL(8):: s(1:3,1:3) + + s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) + + END FUNCTION outerProduct + + !Cross product of two 3D vectors + PURE FUNCTION crossProduct(a, b) RESULT(c) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a, b + REAL(8), DIMENSION(1:3):: c + + c = 0.D0 + + c(1) = a(2)*b(3) - a(3)*b(2) + c(2) = -(a(1)*b(3) - a(3)*b(1)) + c(3) = a(1)*b(2) - a(2)*b(1) + + END FUNCTION crossProduct + + !Normalizes a 3D vector + PURE FUNCTION normalize(a) RESULT(b) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a + REAL(8), DIMENSION(1:3):: b + + b = a / NORM2(a) + + END FUNCTION normalize + +END MODULE moduleMath diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 index b4b8317..963302d 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -31,16 +31,6 @@ MODULE moduleOutput LOGICAL:: emOutput = .FALSE. CONTAINS - FUNCTION outerProduct(a,b) RESULT(s) - IMPLICIT NONE - - REAL(8), DIMENSION(1:3):: a, b - REAL(8):: s(1:3,1:3) - - s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) - - END FUNCTION outerProduct - FUNCTION tensorTrace(a) RESULT(t) IMPLICIT NONE @@ -56,6 +46,7 @@ MODULE moduleOutput USE moduleConstParam USE moduleRefParam USE moduleSpecies + USE moduleMath IMPLICIT NONE TYPE(outputNode), INTENT(in):: rawValues