diff --git a/src/modules/common/moduleConstParam.f90 b/src/modules/common/moduleConstParam.f90 index cb12a23..58a3cc8 100644 --- a/src/modules/common/moduleConstParam.f90 +++ b/src/modules/common/moduleConstParam.f90 @@ -6,7 +6,8 @@ MODULE moduleConstParam REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi REAL(8), PARAMETER:: PI2 = 2.D0*PI !2*pi - REAL(8), PARAMETER:: PI8 = 8.D0*PI !2*pi + REAL(8), PARAMETER:: PI4 = 4.D0*PI !4*pi + REAL(8), PARAMETER:: PI8 = 8.D0*PI !8*pi REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1 REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 0f409fb..ff2cc4b 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -354,7 +354,7 @@ MODULE moduleInput CALL config%get(object // '.file', spFile, found) !Reads node values at the nodes filename = path // spFile - CALL mesh%readInitial(sp, filename, density, velocity, temperature) + CALL mesh%readInitial(filename, density, velocity, temperature) !For each volume in the node, create corresponding particles DO e = 1, mesh%numCells !Scale variables @@ -378,9 +378,9 @@ MODULE moduleInput ALLOCATE(partNew) partNew%species => species(sp)%obj partNew%r = mesh%cells(e)%obj%randPos() - partNew%xi = mesh%cells(e)%obj%phy2log(partNew%r) + partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r) !Get mean velocity at particle position - fPsi = mesh%cells(e)%obj%fPsi(partNew%xi, nNodes) + fPsi = mesh%cells(e)%obj%fPsi(partNew%Xi, nNodes) DO j = 1, nNodes source(j) = velocity(nodes(j), 1) @@ -645,7 +645,7 @@ MODULE moduleInput INTEGER:: e CLASS(meshCell), POINTER:: vol - !Firstly, checks if the object 'interactions' exists + !Firstly, check if the object 'interactions' exists CALL config%info('interactions', found) IF (found) THEN !Checks if MC collisions have been defined diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 6a43ca1..82f43c2 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -41,7 +41,6 @@ MODULE moduleMesh1DCart CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() - REAL(8):: arNodes(1:2) CONTAINS !meshCell DEFERRED PROCEDURES PROCEDURE, PASS:: init => initCell1DCartSegm @@ -60,7 +59,7 @@ MODULE moduleMesh1DCart PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaSegm + PROCEDURE, PASS, PRIVATE:: vol => volumeSegm END TYPE meshCell1DCartSegm @@ -100,7 +99,7 @@ MODULE moduleMesh1DCart END FUNCTION getCoord1DCart !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -133,7 +132,7 @@ MODULE moduleMesh1DCart CALL pointBoundaryFunction(self, s) END DO - + !Physical Surface self%physicalSurface = physicalSurface @@ -162,7 +161,7 @@ MODULE moduleMesh1DCart END FUNCTION intersection1DCart - !Calculates a 'random' position in edge + !Calculate a 'random' position in edge FUNCTION randPosEdge(self) RESULT(r) CLASS(meshEdge1DCart), INTENT(in):: self REAL(8):: r(1:3) @@ -173,7 +172,7 @@ MODULE moduleMesh1DCart !VOLUME FUNCTIONS !SEGMENT FUNCTIONS - !Init segment element + !Init element SUBROUTINE initCell1DCartSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -194,9 +193,7 @@ MODULE moduleMesh1DCart self%x = (/ r1(1), r2(1) /) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -237,7 +234,7 @@ MODULE moduleMesh1DCart END FUNCTION randPos1DCartSegm - !Computes element functions at point Xi + !Compute element functions at point Xi PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -317,7 +314,7 @@ MODULE moduleMesh1DCart END FUNCTION gatherMFSegm - !Computes element local stiffness matrix + !Compute element local stiffness matrix PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) IMPLICIT NONE @@ -348,7 +345,7 @@ MODULE moduleMesh1DCart END FUNCTION elemKSegm - !Computes the local source vector for a force f + !Compute the local source vector for a force f PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) IMPLICIT NONE @@ -422,8 +419,8 @@ MODULE moduleMesh1DCart END SUBROUTINE neighbourElementSegm - !Computes element area - PURE SUBROUTINE areaSegm(self) + !Compute element vol + PURE SUBROUTINE volumeSegm(self) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(inout):: self @@ -433,22 +430,22 @@ MODULE moduleMesh1DCart REAL(8):: fPsi(1:2) self%volume = 0.D0 - self%arNodes = 0.D0 !1D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 2) pDer = self%partialDer(2, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 2) - !Computes total volume of the cell + !Compute total volume of the cell self%volume = detJ*2.D0 - !Computes volume per node - self%arNodes = fPsi*self%volume + !Compute volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume - END SUBROUTINE areaSegm + END SUBROUTINE volumeSegm !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Computes element Jacobian determinant + !Compute element Jacobian determinant PURE FUNCTION detJ1DCart(pDer) RESULT(dJ) IMPLICIT NONE @@ -459,7 +456,7 @@ MODULE moduleMesh1DCart END FUNCTION detJ1DCart - !Computes element Jacobian inverse matrix (without determinant) + !Compute element Jacobian inverse matrix (without determinant) PURE FUNCTION invJ1DCart(pDer) RESULT(invJ) IMPLICIT NONE @@ -575,7 +572,7 @@ MODULE moduleMesh1DCart elemA%e1 => elemB elemB%e2 => elemA - !Revers the normal to point inside the domain + !Rever the normal to point inside the domain elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 8230901..c8ff414 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -41,7 +41,6 @@ MODULE moduleMesh1DRad CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() - REAL(8):: arNodes(1:2) CONTAINS !meshCell DEFERRED PROCEDURES PROCEDURE, PASS:: init => initCell1DRadSegm @@ -60,7 +59,7 @@ MODULE moduleMesh1DRad PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaSegm + PROCEDURE, PASS, PRIVATE:: vol => volumeSegm END TYPE meshCell1DRadSegm @@ -82,7 +81,7 @@ MODULE moduleMesh1DRad !Node volume, to be determined in mesh self%v = 0.D0 - !Allocates output + !Allocate output ALLOCATE(self%output(1:nSpecies)) CALL OMP_INIT_LOCK(self%lock) @@ -100,7 +99,7 @@ MODULE moduleMesh1DRad END FUNCTION getCoord1DRad !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -162,7 +161,7 @@ MODULE moduleMesh1DRad END FUNCTION intersection1DRad - !Calculates a 'random' position in edge + !Calculate a 'random' position in edge FUNCTION randPos1DRad(self) RESULT(r) CLASS(meshEdge1DRad), INTENT(in):: self REAL(8):: r(1:3) @@ -173,7 +172,7 @@ MODULE moduleMesh1DRad !VOLUME FUNCTIONS !SEGMENT FUNCTIONS - !Init segment element + !Init element SUBROUTINE initCell1DRadSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -194,9 +193,7 @@ MODULE moduleMesh1DRad self%r = (/ r1(1), r2(1) /) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -237,7 +234,7 @@ MODULE moduleMesh1DRad END FUNCTION randPos1DRadSegm - !Computes element functions at point Xi + !Compute element functions at point Xi PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -317,7 +314,7 @@ MODULE moduleMesh1DRad END FUNCTION gatherMFSegm - !Computes element local stiffness matrix + !Compute element local stiffness matrix PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE @@ -352,7 +349,7 @@ MODULE moduleMesh1DRad END FUNCTION elemKSegm - !Computes the local source vector for a force f + !Compute the local source vector for a force f PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE @@ -430,9 +427,9 @@ MODULE moduleMesh1DRad END SUBROUTINE neighbourElementSegm - !Computes element area - PURE SUBROUTINE areaSegm(self) - USE moduleConstParam, ONLY: PI + !Compute element vol + PURE SUBROUTINE volumeSegm(self) + USE moduleConstParam, ONLY: PI4 IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(inout):: self @@ -443,28 +440,27 @@ MODULE moduleMesh1DRad REAL(8):: r self%volume = 0.D0 - self%arNodes = 0.D0 !1D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 2) pDer = self%partialDer(2, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 2) - !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) - self%volume = r*detJ*2.D0*PI !2PI - !Computes volume per node + !Compute total volume of the cell + self%volume = r*detJ*PI4 !2*2PI + !Compute volume per node Xi = (/-5.D-1, 0.D0, 0.D0/) r = self%gatherF(Xi, 2, self%r) - self%arNodes(1) = fPsi(1)*self%volume + self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI4 Xi = (/ 5.D-1, 0.D0, 0.D0/) r = self%gatherF(Xi, 2, self%r) - self%arNodes(2) = fPsi(2)*self%volume + self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI4 - END SUBROUTINE areaSegm + END SUBROUTINE volumeSegm !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Computes element Jacobian determinant + !Compute element Jacobian determinant PURE FUNCTION detJ1DRad(pDer) RESULT(dJ) IMPLICIT NONE @@ -475,7 +471,7 @@ MODULE moduleMesh1DRad END FUNCTION detJ1DRad - !Computes element Jacobian inverse matrix (without determinant) + !Compute element Jacobian inverse matrix (without determinant) PURE FUNCTION invJ1DRad(pDer) RESULT(invJ) IMPLICIT NONE @@ -591,7 +587,7 @@ MODULE moduleMesh1DRad elemA%e1 => elemB elemB%e2 => elemA - !Revers the normal to point inside the domain + !Rever the normal to point inside the domain elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index eab1266..cba0cdf 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -47,7 +47,6 @@ MODULE moduleMesh2DCart CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() - REAL(8):: arNodes(1:4) = 0.D0 CONTAINS !meshCell DEFERRED PROCEDURES @@ -67,7 +66,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaQuad + PROCEDURE, PASS, PRIVATE:: vol => volumeQuad END TYPE meshCell2DCartQuad @@ -79,7 +78,6 @@ MODULE moduleMesh2DCart CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() - REAL(8):: arNodes(1:3) = 0.D0 CONTAINS !meshCell DEFERRED PROCEDURES @@ -99,7 +97,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaTria + PROCEDURE, PASS, PRIVATE:: vol => volumeTria END TYPE meshCell2DCartTria @@ -141,7 +139,7 @@ MODULE moduleMesh2DCart END FUNCTION getCoord2DCart !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -198,6 +196,7 @@ MODULE moduleMesh2DCart END FUNCTION getNodes2DCart + !Calculate intersection between position and edge PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r) IMPLICIT NONE @@ -216,7 +215,7 @@ MODULE moduleMesh2DCart END FUNCTION intersection2DCartEdge - !Calculates a random position in edge + !Calculate a random position in edge FUNCTION randPosEdge(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -237,7 +236,7 @@ MODULE moduleMesh2DCart !VOLUME FUNCTIONS !QUAD FUNCTIONS - !Inits quadrilateral element + !Init element SUBROUTINE initCellQuad2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -268,11 +267,7 @@ MODULE moduleMesh2DCart self%y = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - self%n4%v = self%n4%v + self%arNodes(4) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -303,19 +298,19 @@ MODULE moduleMesh2DCart REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) - Xi = 0.D0 Xi(1) = random(-1.D0, 1.D0) Xi(2) = random(-1.D0, 1.D0) + Xi(3) = 0.D0 fPsi = self%fPsi(Xi, 4) - r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) + r(3) = 0.D0 END FUNCTION randPosCellQuad - !Computes element functions in point Xi + !Compute element functions in point Xi PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -323,10 +318,10 @@ MODULE moduleMesh2DCart INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi = (/ (1.D0-Xi(1)) * (1.D0-Xi(2)), & - (1.D0+Xi(1)) * (1.D0-Xi(2)), & - (1.D0+Xi(1)) * (1.D0+Xi(2)), & - (1.D0-Xi(1)) * (1.D0+Xi(2)) /) + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & + (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /) fPsi = fPsi * 0.25D0 @@ -417,7 +412,7 @@ MODULE moduleMesh2DCart END FUNCTION gatherMFQuad - !Computes element local stiffness matrix + !Compute element local stiffness matrix PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) IMPLICIT NONE @@ -427,7 +422,6 @@ MODULE moduleMesh2DCart REAL(8):: Xi(1:3) REAL(8):: dPsi(1:3, 1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: r REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m @@ -478,7 +472,7 @@ MODULE moduleMesh2DCart pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) - f = DOT_PRODUCT(fPsi,source) + f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ END DO @@ -486,7 +480,7 @@ MODULE moduleMesh2DCart END FUNCTION elemFQuad - !Checks if a particle is inside a quad element + !Check if Xi is inside the element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -498,7 +492,7 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transforms physical coordinates to element coordinates + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE @@ -532,7 +526,7 @@ MODULE moduleMesh2DCart END FUNCTION phy2logQuad - !Gets the next element for a logical position Xi + !Get the neighbour element for a logical position Xi SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE @@ -544,7 +538,7 @@ MODULE moduleMesh2DCart XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) nextInt = MAXLOC(XiArray,1) - !Selects the higher value of directions and searches in that direction + !Select the higher value of directions and searches in that direction NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) @@ -559,8 +553,8 @@ MODULE moduleMesh2DCart END SUBROUTINE neighbourElementQuad - !Computes element area - PURE SUBROUTINE areaQuad(self) + !Compute element volume + PURE SUBROUTINE volumeQuad(self) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(inout):: self @@ -570,22 +564,24 @@ MODULE moduleMesh2DCart REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 - self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) - !Computes total volume of the cell - self%volume = detJ - !Computes volume per node - self%arNodes = fPsi*detJ + !Compute total volume of the cell + self%volume = detJ*4.D0 + !Compute volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + self%n4%v = self%n4%v + fPsi(4)*self%volume - END SUBROUTINE areaQuad + END SUBROUTINE volumeQuad - !TRIA ELEMENT - !Init tria element + !TRIA FUNCTIONS + !Init element SUBROUTINE initCellTria2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -613,10 +609,7 @@ MODULE moduleMesh2DCart self%x = (/r1(1), r2(1), r3(1)/) self%y = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -625,7 +618,7 @@ MODULE moduleMesh2DCart END SUBROUTINE initCellTria2DCart - !Gets node indexes from triangular element + !Random position in cell PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) IMPLICIT NONE @@ -637,7 +630,7 @@ MODULE moduleMesh2DCart END FUNCTION getNodesTria - !Random position in quadrilateral volume + !Random position in cell FUNCTION randPosCellTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -659,7 +652,7 @@ MODULE moduleMesh2DCart END FUNCTION randPosCellTria - !Shape functions for triangular element + !Compute element functions in point Xi PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -673,7 +666,7 @@ MODULE moduleMesh2DCart END FUNCTION fPsiTria - !Derivative element function at coordinates Xi + !Compute element derivative functions in point Xi PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE @@ -688,6 +681,7 @@ MODULE moduleMesh2DCart END FUNCTION dPsiTria + !Compute the derivatives in global coordinates PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE @@ -705,6 +699,7 @@ MODULE moduleMesh2DCart END FUNCTION partialDerTria + !Gather electric field at position Xi PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self @@ -720,6 +715,7 @@ MODULE moduleMesh2DCart END FUNCTION gatherEFTria + !Gather magnetic field at position Xi PURE FUNCTION gatherMFTria(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self @@ -743,7 +739,7 @@ MODULE moduleMesh2DCart END FUNCTION gatherMFTria - !Computes element local stiffness matrix + !Compute cell local stiffness matrix PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) IMPLICIT NONE @@ -756,7 +752,8 @@ MODULE moduleMesh2DCart REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l - localK=0.D0 + localK = 0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 @@ -772,7 +769,7 @@ MODULE moduleMesh2DCart END FUNCTION elemKTria - !Computes element local source vector + !Compute element local source vector PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) IMPLICIT NONE @@ -787,22 +784,24 @@ MODULE moduleMesh2DCart INTEGER:: l localF = 0.D0 - Xi = 0.D0 + + Xi = 0.D0 !Start 2D Gauss Quad Integral - DO l=1, 4 + DO l = 1, 4 Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 3) - f = DOT_PRODUCT(fPsi,source) + f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wTria(l)*detJ END DO END FUNCTION elemFTria + !Check if Xi is inside the element PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE @@ -815,7 +814,7 @@ MODULE moduleMesh2DCart END FUNCTION insideTria - !Transforms physical coordinates to element coordinates + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE @@ -838,6 +837,7 @@ MODULE moduleMesh2DCart END FUNCTION phy2logTria + !Get the neighbour cell for a logical position Xi SUBROUTINE neighbourElementTria(self, Xi, neighbourElement) IMPLICIT NONE @@ -861,8 +861,8 @@ MODULE moduleMesh2DCart END SUBROUTINE neighbourElementTria - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) + !Calculate volume for triangular element + PURE SUBROUTINE volumeTria(self) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(inout):: self @@ -872,22 +872,23 @@ MODULE moduleMesh2DCart REAL(8):: fPsi(1:3) self%volume = 0.D0 - self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) - fPsi = self%fPsi(Xi, 4) + fPsi = self%fPsi(Xi, 3) !Computes total volume of the cell - self%volume = detJ + self%volume = detJ !Computes volume per node - self%arNodes = fPsi*detJ + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume - END SUBROUTINE areaTria + END SUBROUTINE volumeTria !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D - !Computes element Jacobian determinant + !Compute element Jacobian determinant PURE FUNCTION detJ2DCart(pDer) RESULT(dJ) IMPLICIT NONE @@ -898,7 +899,7 @@ MODULE moduleMesh2DCart END FUNCTION detJ2DCart - !Computes element Jacobian inverse matrix (without determinant) + !Compute element Jacobian inverse matrix (without determinant) PURE FUNCTION invJ2DCart(pDer) RESULT(invJ) IMPLICIT NONE diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d9a7f32..c2b0674 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -47,7 +47,6 @@ MODULE moduleMesh2DCyl CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() - REAL(8):: arNodes(1:4) = 0.D0 CONTAINS !meshCell DEFERRED PROCEDURES @@ -67,7 +66,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaQuad + PROCEDURE, PASS, PRIVATE:: vol => volumeQuad END TYPE meshCell2DCylQuad @@ -79,7 +78,6 @@ MODULE moduleMesh2DCyl CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() - REAL(8):: arNodes(1:3) = 0.D0 CONTAINS !meshCell DEFERRED PROCEDURES @@ -99,13 +97,13 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: area => areaTria + PROCEDURE, PASS, PRIVATE:: vol => volumeTria END TYPE meshCell2DCylTria CONTAINS !NODE FUNCTIONS - !Inits node element + !Init node element SUBROUTINE initNode2DCyl(self, n, r) USE moduleSpecies USE moduleRefParam @@ -141,7 +139,7 @@ MODULE moduleMesh2DCyl END FUNCTION getCoord2DCyl !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -198,6 +196,7 @@ MODULE moduleMesh2DCyl END FUNCTION getNodes2DCyl + !Calculate intersection between position and edge PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r) IMPLICIT NONE @@ -216,15 +215,15 @@ MODULE moduleMesh2DCyl END FUNCTION intersection2DCylEdge - !Calculates a random position in edge + !Calculate a random position in edge FUNCTION randPosEdge(self) RESULT(r) USE moduleRandom IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd - REAL(8):: dr, dz REAL(8):: r(1:3) + REAL(8):: dr, dz rnd = random() dr = self%r(2) - self%r(1) @@ -245,7 +244,7 @@ MODULE moduleMesh2DCyl !VOLUME FUNCTIONS !QUAD FUNCTIONS - !Inits quadrilateral element + !Init element SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -276,11 +275,7 @@ MODULE moduleMesh2DCyl self%r = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - self%n4%v = self%n4%v + self%arNodes(4) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -289,7 +284,7 @@ MODULE moduleMesh2DCyl END SUBROUTINE initCellQuad2DCyl - !Gets nodes from quadrilateral element + !Get nodes from quadrilateral element PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE @@ -297,7 +292,7 @@ MODULE moduleMesh2DCyl INTEGER, INTENT(in):: nNodes INTEGER:: n(1:nNodes) - n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) + n = (/ self%n1%n, self%n2%n, self%n3%n, self%n4%n /) END FUNCTION getNodesQuad @@ -331,12 +326,12 @@ MODULE moduleMesh2DCyl INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi = (/ (1.D0-Xi(1)) * (1.D0-Xi(2)), & - (1.D0+Xi(1)) * (1.D0-Xi(2)), & - (1.D0+Xi(1)) * (1.D0+Xi(2)), & - (1.D0-Xi(1)) * (1.D0+Xi(2)) /) + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & + (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /) - fPsi = fPsi*0.25D0 + fPsi = fPsi * 0.25D0 END FUNCTION fPsiQuad @@ -350,15 +345,15 @@ MODULE moduleMesh2DCyl dPsi = 0.D0 - dPsi(1,:) = (/ -(1.D0 - Xi(2)), & - (1.D0 - Xi(2)), & - (1.D0 + Xi(2)), & - -(1.D0 + Xi(2)) /) + dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), & + (1.D0 - Xi(2)), & + (1.D0 + Xi(2)), & + -(1.D0 + Xi(2)) /) - dPsi(2,:) = (/ -(1.D0 - Xi(1)), & - -(1.D0 + Xi(1)), & - (1.D0 + Xi(1)), & - (1.D0 - Xi(1)) /) + dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), & + -(1.D0 + Xi(1)), & + (1.D0 + Xi(1)), & + (1.D0 - Xi(1)) /) dPsi = dPsi * 0.25D0 @@ -379,6 +374,7 @@ MODULE moduleMesh2DCyl DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /) pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), & DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /) + pDer(3,3) = 1.D0 END FUNCTION partialDerQuad @@ -439,10 +435,11 @@ MODULE moduleMesh2DCyl REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m - localK=0.D0 - Xi=0.D0 + localK = 0.D0 + + Xi = 0.D0 !Start 2D Gauss Quad Integral - DO l=1, 3 + DO l = 1, 3 Xi(2) = corQuad(l) DO m = 1, 3 Xi(1) = corQuad(m) @@ -479,8 +476,9 @@ MODULE moduleMesh2DCyl INTEGER:: l, m localF = 0.D0 - Xi = 0.D0 - DO l=1, 3 + + Xi = 0.D0 + DO l = 1, 3 Xi(1) = corQuad(l) DO m = 1, 3 Xi(2) = corQuad(m) @@ -489,7 +487,7 @@ MODULE moduleMesh2DCyl detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi, self%r) - f = DOT_PRODUCT(fPsi,source) + f = DOT_PRODUCT(fPsi, source) localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ END DO @@ -498,7 +496,7 @@ MODULE moduleMesh2DCyl END FUNCTION elemFQuad - !Checks if a particle is inside a quad element + !Checks if Xi is inside the element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -510,7 +508,7 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transforms physical coordinates to element coordinates + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE @@ -544,7 +542,7 @@ MODULE moduleMesh2DCyl END FUNCTION phy2logQuad - !Get the next element for a logical position Xi + !Get the neighbour element for a logical position Xi SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE @@ -571,8 +569,8 @@ MODULE moduleMesh2DCyl END SUBROUTINE neighbourElementQuad - !Computes element area - PURE SUBROUTINE areaQuad(self) + !Compute element volume + PURE SUBROUTINE volumeQuad(self) USE moduleConstParam, ONLY: PI8 IMPLICIT NONE @@ -584,34 +582,33 @@ MODULE moduleMesh2DCyl REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 - self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) - !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) + !Computes total volume of the cell self%volume = r*detJ*PI8 !4*2*pi !Computes volume per node Xi = (/-5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(1) = fPsi(1)*self%volume + self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 Xi = (/ 5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(2) = fPsi(2)*self%volume + self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 Xi = (/ 5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(3) = fPsi(3)*self%volume + self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 Xi = (/-5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(4) = fPsi(4)*self%volume + self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 - END SUBROUTINE areaQuad + END SUBROUTINE volumeQuad - !TRIA ELEMENT - !Init tria element + !TRIA FUNCTIONS + !Init element SUBROUTINE initCellTria2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -639,10 +636,7 @@ MODULE moduleMesh2DCyl self%z = (/r1(1), r2(1), r3(1)/) self%r = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -651,7 +645,7 @@ MODULE moduleMesh2DCyl END SUBROUTINE initCellTria2DCyl - !Gets node indexes from triangular element + !Random position in cell PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) IMPLICIT NONE @@ -663,7 +657,7 @@ MODULE moduleMesh2DCyl END FUNCTION getNodesTria - !Random position in quadrilateral volume + !Random position in cell FUNCTION randPosCellTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -685,7 +679,7 @@ MODULE moduleMesh2DCyl END FUNCTION randPosCellTria - !Shape functions for triangular element + !Compute element functions in point Xi PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -699,7 +693,7 @@ MODULE moduleMesh2DCyl END FUNCTION fPsiTria - !Derivative element function at coordinates Xi + !Compute element derivative functions in point Xi PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE @@ -714,6 +708,7 @@ MODULE moduleMesh2DCyl END FUNCTION dPsiTria + !Compute the derivatives in global coordinates PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE @@ -731,6 +726,7 @@ MODULE moduleMesh2DCyl END FUNCTION partialDerTria + !Gather electric field at position Xi PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self @@ -746,6 +742,7 @@ MODULE moduleMesh2DCyl END FUNCTION gatherEFTria + !Gather magnetic field at position Xi PURE FUNCTION gatherMFTria(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self @@ -769,7 +766,7 @@ MODULE moduleMesh2DCyl END FUNCTION gatherMFTria - !Computes element local stiffness matrix + !Compute cell local stiffness matrix PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE @@ -784,7 +781,8 @@ MODULE moduleMesh2DCyl REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l - localK=0.D0 + localK = 0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 @@ -802,7 +800,7 @@ MODULE moduleMesh2DCyl END FUNCTION elemKTria - !Computes element local source vector + !Compute element local source vector PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE @@ -838,6 +836,7 @@ MODULE moduleMesh2DCyl END FUNCTION elemFTria + !Check if Xi is inside the element PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE @@ -850,7 +849,7 @@ MODULE moduleMesh2DCyl END FUNCTION insideTria - !Transforms physical coordinates to element coordinates + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE @@ -873,6 +872,7 @@ MODULE moduleMesh2DCyl END FUNCTION phy2logTria + !Get the neighbour cell for a logical position Xi SUBROUTINE neighbourElementTria(self, Xi, neighbourElement) IMPLICIT NONE @@ -896,8 +896,8 @@ MODULE moduleMesh2DCyl END SUBROUTINE neighbourElementTria - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) + !Calculate volume for triangular element + PURE SUBROUTINE volumeTria(self) USE moduleConstParam, ONLY: PI IMPLICIT NONE @@ -909,23 +909,24 @@ MODULE moduleMesh2DCyl REAL(8):: r self%volume = 0.D0 - self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 3) - !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) + !Computes total volume of the cell self%volume = r*detJ*PI !2PI*1/2 !Computes volume per node - self%arNodes = fPsi*self%volume + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume - END SUBROUTINE areaTria + END SUBROUTINE volumeTria !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS - !Computes element Jacobian determinant + !Compute element Jacobian determinant PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ) IMPLICIT NONE @@ -936,7 +937,7 @@ MODULE moduleMesh2DCyl END FUNCTION detJ2DCyl - !Computes element Jacobian inverse matrix (without determinant) + !Compute element Jacobian inverse matrix (without determinant) PURE FUNCTION invJ2DCyl(pDer) RESULT(invJ) IMPLICIT NONE diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 34474cd..fcd4647 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -60,13 +60,13 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: neighbourElement => neighbourElementTetra !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: calcVol => volumeTetra + PROCEDURE, PASS, PRIVATE:: vol => volumeTetra END TYPE meshCell3DCartTetra CONTAINS !NODE FUNCTIONS - !Inits node element + !Init node element SUBROUTINE initNode3DCart(self, n, r) USE moduleSpecies USE moduleRefParam @@ -102,8 +102,8 @@ MODULE moduleMesh3DCart END FUNCTION getCoord3DCart - !SURFACE FUNCTIONS - !Inits surface element + !EDGE FUNCTIONS + !Init surface element SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -168,6 +168,7 @@ MODULE moduleMesh3DCart END FUNCTION getNodes3DCartTria + !Calculate intersection between position and edge PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r) IMPLICIT NONE @@ -186,7 +187,7 @@ MODULE moduleMesh3DCart END FUNCTION intersection3DCartTria - !Calculates a random position in the surface + !Calculate a random position in the surface FUNCTION randPosEdgeTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -222,7 +223,7 @@ MODULE moduleMesh3DCart !VOLUME FUNCTIONS !TETRA FUNCTIONS - !Inits tetrahedron element + !Init element SUBROUTINE initCellTetra(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -232,11 +233,14 @@ MODULE moduleMesh3DCart INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node - REAL(8):: Xi(1:3), fPsi(1:4) - REAL(8):: volNodes(1:4) !Cellume of each node + !Assign node index self%n = n + + !Assign number of nodes of cell self%nNodes = SIZE(p) + + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -251,16 +255,7 @@ MODULE moduleMesh3DCart self%z = (/r1(3), r2(3), r3(3), r4(3)/) !Computes the element volume - CALL self%calcVol() - - !Assign proportional volume to each node - Xi = (/0.25D0, 0.25D0, 0.25D0/) - fPsi = self%fPsi(Xi, 4) - volNodes = fPsi*self%volume - self%n1%v = self%n1%v + volNodes(1) - self%n2%v = self%n2%v + volNodes(2) - self%n3%v = self%n3%v + volNodes(3) - self%n4%v = self%n4%v + volNodes(4) + CALL self%vol() CALL OMP_INIT_LOCK(self%lock) @@ -269,6 +264,7 @@ MODULE moduleMesh3DCart END SUBROUTINE initCellTetra + !Gets node indexes from cell PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n) IMPLICIT NONE @@ -280,7 +276,7 @@ MODULE moduleMesh3DCart END FUNCTION getNodesTetra - !Random position in volume tetrahedron + !Random position in cell FUNCTION randPosCellTetra(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -302,7 +298,7 @@ MODULE moduleMesh3DCart END FUNCTION randPosCellTetra - !Computes element functions in point Xi + !Compute element functions in point Xi PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE @@ -317,7 +313,7 @@ MODULE moduleMesh3DCart END FUNCTION fPsiTetra - !Derivative element function at coordinates Xi + !Compute element derivative functions in point Xi PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE @@ -333,7 +329,7 @@ MODULE moduleMesh3DCart END FUNCTION dPsiTetra - !Computes the derivatives in global coordinates + !Compute the derivatives in global coordinates PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE @@ -358,6 +354,7 @@ MODULE moduleMesh3DCart END FUNCTION partialDerTetra + !Gather electric field at position Xi PURE FUNCTION gatherEFTetra(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self @@ -374,6 +371,7 @@ MODULE moduleMesh3DCart END FUNCTION gatherEFTetra + !Gather magnetic field at position Xi PURE FUNCTION gatherMFTetra(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self @@ -400,6 +398,7 @@ MODULE moduleMesh3DCart END FUNCTION gatherMFTetra + !Compute cell local stiffness matrix PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK) IMPLICIT NONE @@ -424,6 +423,7 @@ MODULE moduleMesh3DCart END FUNCTION elemKTetra + !Compute element local source vector PURE FUNCTION elemFTetra(self, nNodes, source) RESULT(localF) IMPLICIT NONE @@ -448,6 +448,7 @@ MODULE moduleMesh3DCart END FUNCTION elemFTetra + !Check if Xi is inside the element PURE FUNCTION insideTetra(Xi) RESULT(ins) IMPLICIT NONE @@ -461,6 +462,7 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra + !Transform physical coordinates to element coordinates PURE FUNCTION phy2logTetra(self,r) RESULT(Xi) IMPLICIT NONE @@ -472,6 +474,7 @@ MODULE moduleMesh3DCart REAL(8):: invJ(1:3, 1:3), detJ REAL(8):: deltaR(1:3) + !Direct method to convert coordinates Xi = 0.D0 deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /) dPsi = self%dPsi(Xi, 4) @@ -482,6 +485,7 @@ MODULE moduleMesh3DCart END FUNCTION phy2logTetra + !Get the neighbour cell for a logical position Xi SUBROUTINE neighbourElementTetra(self, Xi, neighbourElement) IMPLICIT NONE @@ -508,25 +512,35 @@ MODULE moduleMesh3DCart END SUBROUTINE neighbourElementTetra - !Computes the element volume + !Calculate volume for triangular element PURE SUBROUTINE volumeTetra(self) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(inout):: self REAL(8):: Xi(1:3) - REAL(8):: dPsi(1:3, 1:4) - REAL(8):: pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 + !2D 1 point Gauss Quad Integral Xi = (/0.25D0, 0.25D0, 0.25D0/) dPsi = self%dPsi(Xi, 4) pDer = self%partialDer(4, dPsi) - self%volume = self%detJac(pDer) + detJ = self%detJac(pDer) + !Computes total volume of the cell + self%volume = detJ + !Computes volume per node + fPsi = self%fPsi(Xi, 4) + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + self%n4%v = self%n4%v + fPsi(4)*self%volume END SUBROUTINE volumeTetra !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D - !Computes element Jacobian determinant + !Compute element Jacobian determinant PURE FUNCTION detJ3DCart(pDer) RESULT(dJ) IMPLICIT NONE @@ -539,6 +553,7 @@ MODULE moduleMesh3DCart END FUNCTION detJ3DCart + !Compute element Jacobian inverse matrix (without determinant) PURE FUNCTION invJ3DCart(pDer) RESULT(invJ) IMPLICIT NONE @@ -561,7 +576,37 @@ MODULE moduleMesh3DCart END FUNCTION invJ3DCart - !Selects type of elements to build connection + SUBROUTINE connectMesh3DCart(self) + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(inout):: self + INTEGER:: e, et + + DO e = 1, self%numCells + !Connect Cell-Cell + DO et = 1, self%numCells + IF (e /= et) THEN + CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj) + + END IF + + END DO + + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Cell-Edge + DO et = 1, self%numEdges + CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj) + + END DO + + END SELECT + + END DO + + END SUBROUTINE connectMesh3DCart + + !Select type of elements to build connection SUBROUTINE connectCellCell(elemA, elemB) IMPLICIT NONE @@ -601,36 +646,6 @@ MODULE moduleMesh3DCart END SUBROUTINE connectCellEdge - SUBROUTINE connectMesh3DCart(self) - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(inout):: self - INTEGER:: e, et - - DO e = 1, self%numCells - !Connect Cell-Cell - DO et = 1, self%numCells - IF (e /= et) THEN - CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj) - - END IF - - END DO - - SELECT TYPE(self) - TYPE IS(meshParticles) - !Connect Cell-Edge - DO et = 1, self%numEdges - CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj) - - END DO - - END SELECT - - END DO - - END SUBROUTINE connectMesh3DCart - !Checks if two sets of nodes are coincidend in any order PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) IMPLICIT NONE diff --git a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 index 37dbf82..d968ba3 100644 --- a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 @@ -64,10 +64,9 @@ MODULE moduleMeshInput0D END SUBROUTINE read0D - SUBROUTINE readInitial0D(sp, filename, density, velocity, temperature) + SUBROUTINE readInitial0D(filename, density, velocity, temperature) IMPLICIT NONE - INTEGER, INTENT(in):: sp CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index aae2216..ae1fe05 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -321,10 +321,9 @@ MODULE moduleMeshInputGmsh2 END SUBROUTINE readGmsh2 !Reads the initial information from an output file for an species - SUBROUTINE readInitialGmsh2(sp, filename, density, velocity, temperature) + SUBROUTINE readInitialGmsh2(filename, density, velocity, temperature) IMPLICIT NONE - INTEGER, INTENT(in):: sp CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 97ce691..f15f880 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -356,8 +356,7 @@ MODULE moduleMesh END SUBROUTINE readMesh_interface - SUBROUTINE readInitial_interface(sp, filename, density, velocity, temperature) - INTEGER, INTENT(in):: sp + SUBROUTINE readInitial_interface(filename, density, velocity, temperature) CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 713c091..517835e 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -159,8 +159,8 @@ MODULE moduleMeshBoundary newElectron%vol = part%vol newIon%vol = part%vol - newElectron%xi = mesh%cells(part%vol)%obj%phy2log(newElectron%r) - newIon%xi = newElectron%xi + newElectron%Xi = mesh%cells(part%vol)%obj%phy2log(newElectron%r) + newIon%Xi = newElectron%Xi newElectron%weight = part%weight newIon%weight = newElectron%weight diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 224bfbb..ccca930 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -111,13 +111,13 @@ MODULE moduleCollisions IMPLICIT NONE REAL(8):: n(1:3) - REAL(8):: cosXii, sinXii, eps + REAL(8):: cosXi, sinXi, eps - cosXii = random(-1.D0, 1.D0) - sinXii = DSQRT(1.D0 - cosXii**2) + cosXi = random(-1.D0, 1.D0) + sinXi = DSQRT(1.D0 - cosXi**2) eps = random(0.D0, PI2) - n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /) + n = (/ cosXi, sinXi*DCOS(eps), sinXi*DSIN(eps) /) END FUNCTION randomDirectionVHS diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index d19ff28..ca7858c 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -40,7 +40,7 @@ MODULE moduleSpecies CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle INTEGER:: vol !Index of element in which the particle is located INTEGER:: volColl !Index of element in which the particle is located in the Collision Mesh - REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p. + REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle