diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index aaf4b08..3ff9e5b 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -337,7 +337,7 @@ MODULE moduleInput !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi INTEGER:: nNewPart = 0.D0 - CLASS(meshVol), POINTER:: vol + CLASS(meshCell), POINTER:: vol TYPE(particle), POINTER:: partNew REAL(8):: vTh TYPE(lNode), POINTER:: partCurr, partNext @@ -356,13 +356,13 @@ MODULE moduleInput filename = path // spFile CALL mesh%readInitial(sp, filename, density, velocity, temperature) !For each volume in the node, create corresponding particles - DO e = 1, mesh%numVols + DO e = 1, mesh%numCells !Scale variables !Density at centroid of cell - nodes = mesh%vols(e)%obj%getNodes() - nNodes = SIZE(nodes) + nodes = mesh%cells(e)%obj%getNodes() + nNodes = mesh%cells(e)%obj%nNodes ALLOCATE(fPsi(1:nNodes)) - CALL mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), fPsi) + fPsi = mesh%cells(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) ALLOCATE(source(1:nNodes)) DO j = 1, nNodes source(j) = density(nodes(j)) @@ -371,16 +371,16 @@ MODULE moduleInput densityCen = DOT_PRODUCT(fPsi, source) !Calculate number of particles - nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) !Allocate new particles DO p = 1, nNewPart ALLOCATE(partNew) partNew%species => species(sp)%obj - partNew%r = mesh%vols(e)%obj%randPos() - partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) + partNew%r = mesh%cells(e)%obj%randPos() + partNew%xi = mesh%cells(e)%obj%phy2log(partNew%r) !Get mean velocity at particle position - CALL mesh%vols(e)%obj%fPsi(partNew%xi, fPsi) + fPsi = mesh%cells(e)%obj%fPsi(partNew%xi) DO j = 1, nNodes source(j) = velocity(nodes(j), 1) @@ -426,7 +426,7 @@ MODULE moduleInput CALL partInitial%add(partNew) !Assign particle to list in volume - vol => meshforMCC%vols(partNew%volColl)%obj + vol => meshforMCC%cells(partNew%volColl)%obj CALL OMP_SET_LOCK(vol%lock) CALL vol%listPart_in(sp)%add(partNew) vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight @@ -643,7 +643,7 @@ MODULE moduleInput REAL(8):: energyThreshold, energyBinding CHARACTER(:), ALLOCATABLE:: electron INTEGER:: e - CLASS(meshVol), POINTER:: vol + CLASS(meshCell), POINTER:: vol !Firstly, checks if the object 'interactions' exists CALL config%info('interactions', found) @@ -739,8 +739,8 @@ MODULE moduleInput END DO !Init the required arrays in each volume to account for MCC. - DO e = 1, meshForMCC%numVols - vol => meshForMCC%vols(e)%obj + DO e = 1, meshForMCC%numCells + vol => meshForMCC%cells(e)%obj !Allocate Maximum cross section per collision pair and assign the initial collision rate ALLOCATE(vol%sigmaVrelMax(1:nCollPairs)) @@ -930,8 +930,8 @@ MODULE moduleInput CALL config%get(object // '.volume', volume, found) !Rescale the volumne IF (found) THEN - mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref - mesh%nodes(1)%obj%v = mesh%vols(1)%obj%volume + mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref + mesh%nodes(1)%obj%v = mesh%cells(1)%obj%volume END IF diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 0a14520..a57f244 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -11,22 +11,25 @@ MODULE moduleMesh0D END TYPE meshNode0D - TYPE, PUBLIC, EXTENDS(meshVol):: meshVol0D + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell0D CLASS(meshNode), POINTER:: n1 CONTAINS - PROCEDURE, PASS:: init => initVol0D + PROCEDURE, PASS:: init => initCell0D PROCEDURE, PASS:: getNodes => getNodes0D PROCEDURE, PASS:: randPos => randPos0D - PROCEDURE, NOPASS:: fPsi => fPsi0D - PROCEDURE, PASS:: gatherEF => gatherEF0D - PROCEDURE, PASS:: gatherMF => gatherMF0D + PROCEDURE, PASS:: fPsi => fPsi0D + PROCEDURE, PASS:: dPsi => dPsi0D + PROCEDURE, PASS:: detJac => detJ0D + PROCEDURE, PASS:: invJac => invJ0D PROCEDURE, PASS:: elemK => elemK0D PROCEDURE, PASS:: elemF => elemF0D + PROCEDURE, PASS:: gatherElectricField => gatherEF0D + PROCEDURE, PASS:: gatherMagneticField => gatherMF0D PROCEDURE, PASS:: phy2log => phy2log0D PROCEDURE, NOPASS:: inside => inside0D PROCEDURE, PASS:: nextElement => nextElement0D - END TYPE meshVol0D + END TYPE meshCell0D CONTAINS !NODE FUNCTIONS @@ -61,18 +64,20 @@ MODULE moduleMesh0D !VOLUME FUNCTIONS !Inits dummy 0D volume - SUBROUTINE initVol0D(self, n, p, nodes) + SUBROUTINE initCell0D(self, n, p, nodes) USE moduleRefParam USE moduleSpecies IMPLICIT NONE - CLASS(meshVol0D), INTENT(out):: self + CLASS(meshCell0D), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) self%n = n + self%nNodes = SIZE(p) + self%n1 => nodes(p(1))%obj self%volume = 1.D0 self%n1%v = 1.D0 @@ -82,15 +87,13 @@ MODULE moduleMesh0D ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVol0D + END SUBROUTINE initCell0D PURE FUNCTION getNodes0D(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:1)) + CLASS(meshCell0D), INTENT(in):: self + INTEGER:: n(1:self%nNodes) n = self%n1%n @@ -99,50 +102,65 @@ MODULE moduleMesh0D FUNCTION randPos0D(self) RESULT(r) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self + CLASS(meshCell0D), INTENT(in):: self REAL(8):: r(1:3) r = 0.D0 END FUNCTION randPos0D - PURE SUBROUTINE fPsi0D(xi, fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + PURE FUNCTION fPsi0D(self, Xi) RESULT(fPsi) + IMPLICIT NONE + + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: fPsi(1:self%nNodes) fPsi = 1.D0 - END SUBROUTINE fPsi0D + END FUNCTION fPsi0D - PURE FUNCTION gatherEF0D(self, xi) RESULT(EF) + PURE FUNCTION dPsi0D(self, Xi) RESULT(dPsi) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: EF(1:3) + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: dPsi(1:3,1:self%nNodes) - EF = 0.D0 + dPsi = 0.D0 - END FUNCTION gatherEF0D + END FUNCTION dPsi0D - PURE FUNCTION gatherMF0D(self, xi) RESULT(MF) + PURE FUNCTION detJ0D(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: MF(1:3) + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dJ - MF = 0.D0 + dJ = 0.D0 - END FUNCTION gatherMF0D + END FUNCTION detJ0D + + PURE FUNCTION invJ0D(self, Xi, dPsi_in) RESULT(invJ) + IMPLICIT NONE + + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: invJ(1:3,1:3) + + invJ = 0.D0 + + END FUNCTION invJ0D PURE FUNCTION elemK0D(self) RESULT(localK) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell0D), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) - ALLOCATE(localK(1:1, 1:1)) localK = 0.D0 END FUNCTION elemK0D @@ -150,19 +168,48 @@ MODULE moduleMesh0D PURE FUNCTION elemF0D(self, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) - ALLOCATE(localF(1:1)) localF = 0.D0 END FUNCTION elemF0D + PURE FUNCTION gatherEF0D(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:1) + + phi = (/ self%n1%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEF0D + + PURE FUNCTION gatherMF0D(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:1,1:3) + + B(:,1) = (/ self%n1%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMF0D + PURE FUNCTION phy2log0D(self,r) RESULT(xN) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self + CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: xN(1:3) @@ -170,21 +217,21 @@ MODULE moduleMesh0D END FUNCTION phy2log0D - PURE FUNCTION inside0D(xi) RESULT(ins) + PURE FUNCTION inside0D(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins ins = .TRUE. END FUNCTION inside0D - SUBROUTINE nextElement0D(self, xi, nextElement) + SUBROUTINE nextElement0D(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement nextElement => NULL() diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 3380cf0..0f46ba5 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -8,7 +8,7 @@ MODULE moduleMesh1DCart 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 /) + REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /) TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart !Element coordinates @@ -32,33 +32,26 @@ MODULE moduleMesh1DCart END TYPE meshEdge1DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DCart CONTAINS PROCEDURE, PASS:: detJac => detJ1DCart PROCEDURE, PASS:: invJac => invJ1DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - END TYPE meshVol1DCart + END TYPE meshCell1DCart ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - PURE SUBROUTINE partialDer_interface(self, dPsi, dx) - IMPORT meshVol1DCart - CLASS(meshVol1DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + IMPORT meshCell1DCart + CLASS(meshCell1DCart), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx END SUBROUTINE partialDer_interface END INTERFACE - TYPE, PUBLIC, EXTENDS(meshVol1DCart):: meshVol1DCartSegm + TYPE, PUBLIC, EXTENDS(meshCell1DCart):: meshCell1DCartSegm !Element coordinates REAL(8):: x(1:2) !Connectivity to nodes @@ -67,22 +60,22 @@ MODULE moduleMesh1DCart CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initVol1DCartSegm - PROCEDURE, PASS:: randPos => randPos1DCartSeg + PROCEDURE, PASS:: init => initCell1DCartSegm + PROCEDURE, PASS:: randPos => randPos1DCartSegm PROCEDURE, PASS:: area => areaSegm - PROCEDURE, NOPASS:: fPsi => fPsiSegm - PROCEDURE, NOPASS:: dPsi => dPsiSegm + PROCEDURE, PASS:: fPsi => fPsiSegm + PROCEDURE, PASS:: dPsi => dPsiSegm PROCEDURE, PASS:: partialDer => partialDerSegm PROCEDURE, PASS:: elemK => elemKSegm PROCEDURE, PASS:: elemF => elemFSegm + PROCEDURE, PASS:: gatherElectricField => gatherEFSegm + PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm PROCEDURE, NOPASS:: inside => insideSegm - PROCEDURE, PASS:: gatherEF => gatherEFSegm - PROCEDURE, PASS:: gatherMF => gatherMFSegm PROCEDURE, PASS:: getNodes => getNodesSegm PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: nextElement => nextElementSegm - END TYPE meshVol1DCartSegm + END TYPE meshCell1DCartSegm CONTAINS !NODE FUNCTIONS @@ -193,17 +186,18 @@ MODULE moduleMesh1DCart !VOLUME FUNCTIONS !SEGMENT FUNCTIONS !Init segment element - SUBROUTINE initVol1DCartSegm(self, n, p, nodes) + SUBROUTINE initCell1DCartSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(out):: self + CLASS(meshCell1DCartSegm), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2 self%n = n + self%nNodes = SIZE(p) self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj !Get element coordinates @@ -221,14 +215,14 @@ MODULE moduleMesh1DCart ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVol1DCartSegm + END SUBROUTINE initCell1DCartSegm !Calculates a random position in 1D volume - FUNCTION randPos1DCartSeg(self) RESULT(r) + FUNCTION randPos1DCartSegm(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self + CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) REAL(8), ALLOCATABLE:: fPsi(:) @@ -236,16 +230,16 @@ MODULE moduleMesh1DCart Xi(1) = random(-1.D0, 1.D0) Xi(2:3) = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%x) - END FUNCTION randPos1DCartSeg + END FUNCTION randPos1DCartSegm !Computes element area PURE SUBROUTINE areaSegm(self) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(inout):: self + CLASS(meshCell1DCartSegm), INTENT(inout):: self REAL(8):: l !element length REAL(8):: fPsi(1:2) REAL(8):: detJ @@ -255,7 +249,7 @@ MODULE moduleMesh1DCart self%arNodes = 0.D0 !1 point Gauss integral Xi = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) detJ = self%detJac(Xi) l = 2.D0*detJ self%volume = l @@ -264,26 +258,29 @@ MODULE moduleMesh1DCart END SUBROUTINE areaSegm !Computes element functions at point Xi - PURE SUBROUTINE fPsiSegm(xi, fPsi) + PURE FUNCTION fPsiSegm(self, xi) RESULT(fPsi) IMPLICIT NONE + CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + REAL(8):: fPsi(1:self%nNodes) fPsi(1) = 1.D0 - xi(1) fPsi(2) = 1.D0 + xi(1) + fPsi = fPsi * 5.D-1 - END SUBROUTINE fPsiSegm + END FUNCTION fPsiSegm !Computes element derivative shape function at Xi - PURE FUNCTION dPsiSegm(xi) RESULT(dPsi) + PURE FUNCTION dPsiSegm(self, xi) RESULT(dPsi) IMPLICIT NONE + CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:1, 1:2)) + dPsi = 0.D0 dPsi(1, 1) = -5.D-1 dPsi(1, 2) = 5.D-1 @@ -294,8 +291,8 @@ MODULE moduleMesh1DCart PURE SUBROUTINE partialDerSegm(self, dPsi, dx) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) @@ -306,14 +303,13 @@ MODULE moduleMesh1DCart PURE FUNCTION elemKSegm(self) RESULT(localK) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) REAL(8):: Xi(1:3) - REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ(1), detJ + REAL(8):: dPsi(1:3, 1:2) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l - ALLOCATE(localK(1:2,1:2)) localK = 0.D0 Xi = 0.D0 DO l = 1, 3 @@ -332,22 +328,21 @@ MODULE moduleMesh1DCart PURE FUNCTION elemFSegm(self, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: fPsi(1:2) REAL(8):: detJ, f REAL(8):: Xi(1:3) INTEGER:: l - ALLOCATE(localF(1:2)) localF = 0.D0 Xi = 0.D0 DO l = 1, 3 Xi(1) = corSeg(l) detJ = self%detJac(Xi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wSeg(l)*detJ @@ -355,6 +350,40 @@ MODULE moduleMesh1DCart END FUNCTION elemFSegm + PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:2) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEFSegm + + PURE FUNCTION gatherMFSegm(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:2,1:3) + + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFSegm + PURE FUNCTION insideSegm(xi) RESULT(ins) IMPLICIT NONE @@ -366,58 +395,13 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm - !Gathers EF at position Xi - PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1, 1:2) - REAL(8):: phi(1:2) - REAL(8):: EF(1:3) - REAL(8):: invJ - - phi = (/ self%n1%emData%phi, & - self%n2%emData%phi /) - - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ - EF(2) = 0.D0 - EF(3) = 0.D0 - - END FUNCTION gatherEFSegm - - PURE FUNCTION gatherMFSegm(self, xi) RESULT(MF) - IMPLICIT NONE - - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:2) - REAL(8):: MF_Nodes(1:2, 1:3) - REAL(8):: MF(1:3) - REAL(8):: invJ - - MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), & - self%n2%emData%B(1) /) - MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), & - self%n2%emData%B(2) /) - MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & - self%n2%emData%B(3) /) - - CALL self%fPsi(xi, fPsi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFSegm - !Get nodes from 1D volume PURE FUNCTION getNodesSegm(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:2)) n = (/ self%n1%n, self%n2%n /) END FUNCTION getNodesSegm @@ -425,7 +409,7 @@ MODULE moduleMesh1DCart PURE FUNCTION phy2logSegm(self, r) RESULT(xN) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self + CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: xN(1:3) @@ -438,7 +422,7 @@ MODULE moduleMesh1DCart SUBROUTINE nextElementSegm(self, xi, nextElement) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self + CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement @@ -459,10 +443,10 @@ MODULE moduleMesh1DCart PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol1DCart), INTENT(in):: self + CLASS(meshCell1DCart), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dJ REAL(8):: dx(1) @@ -483,12 +467,12 @@ MODULE moduleMesh1DCart PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol1DCart), INTENT(in):: self + CLASS(meshCell1DCart), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: invJ(1:3,1:3) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dx(1) - REAL(8):: invJ IF (PRESENT(dPsi_in)) THEN dPsi = dPsi_in @@ -498,8 +482,11 @@ MODULE moduleMesh1DCart END IF + invJ = 0.D0 + CALL self%partialDer(dPsi, dx) - invJ = 1.D0/dx(1) + + invJ(1,1) = 1.D0/dx(1) END FUNCTION invJ1DCart @@ -509,11 +496,11 @@ MODULE moduleMesh1DCart CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et - DO e = 1, self%numVols + DO e = 1, self%numCells !Connect Vol-Vol - DO et = 1, self%numVols + DO et = 1, self%numCells IF (e /= et) THEN - CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj) END IF @@ -523,7 +510,7 @@ MODULE moduleMesh1DCart TYPE IS(meshParticles) !Connect Vol-Edge DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj) END DO @@ -536,13 +523,13 @@ MODULE moduleMesh1DCart SUBROUTINE connectVolVol(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVol1DCartSegm) + TYPE IS(meshCell1DCartSegm) SELECT TYPE(elemB) - TYPE IS(meshVol1DCartSegm) + TYPE IS(meshCell1DCartSegm) CALL connectSegmSegm(elemA, elemB) END SELECT @@ -554,8 +541,8 @@ MODULE moduleMesh1DCart SUBROUTINE connectSegmSegm(elemA, elemB) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB + CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemA + CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemB IF (.NOT. ASSOCIATED(elemA%e1) .AND. & elemA%n2%n == elemB%n1%n) THEN @@ -577,11 +564,11 @@ MODULE moduleMesh1DCart SUBROUTINE connectVolEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemA CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS (meshVol1DCartSegm) + TYPE IS (meshCell1DCartSegm) SELECT TYPE(elemB) CLASS IS(meshEdge1DCart) CALL connectSegmEdge(elemA, elemB) @@ -595,7 +582,7 @@ MODULE moduleMesh1DCart SUBROUTINE connectSegmEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA + CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemA CLASS(meshEdge1DCart), INTENT(inout), TARGET:: elemB IF (.NOT. ASSOCIATED(elemA%e1) .AND. & diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 7b09e5b..8ebac17 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -32,34 +32,27 @@ MODULE moduleMesh1DRad END TYPE meshEdge1DRad - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DRad CONTAINS PROCEDURE, PASS:: detJac => detJ1DRad PROCEDURE, PASS:: invJac => invJ1DRad - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - END TYPE meshVol1DRad + END TYPE meshCell1DRad ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - PURE SUBROUTINE partialDer_interface(self, dPsi, dx) - IMPORT meshVol1DRad + IMPORT meshCell1DRad - CLASS(meshVol1DRad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell1DRad), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx END SUBROUTINE partialDer_interface END INTERFACE - TYPE, PUBLIC, EXTENDS(meshVol1DRad):: meshVol1DRadSegm + TYPE, PUBLIC, EXTENDS(meshCell1DRad):: meshCell1DRadSegm !Element coordinates REAL(8):: r(1:2) !Connectivity to nodes @@ -68,22 +61,22 @@ MODULE moduleMesh1DRad CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initVol1DRadSegm + PROCEDURE, PASS:: init => initCell1DRadSegm PROCEDURE, PASS:: randPos => randPos1DRadSeg PROCEDURE, PASS:: area => areaRad - PROCEDURE, NOPASS:: fPsi => fPsiRad - PROCEDURE, NOPASS:: dPsi => dPsiRad + PROCEDURE, PASS:: fPsi => fPsiRad + PROCEDURE, PASS:: dPsi => dPsiRad PROCEDURE, PASS:: partialDer => partialDerRad PROCEDURE, PASS:: elemK => elemKRad PROCEDURE, PASS:: elemF => elemFRad + PROCEDURE, PASS:: gatherElectricField => gatherEFRad + PROCEDURE, PASS:: gatherMagneticField => gatherMFRad PROCEDURE, NOPASS:: inside => insideRad - PROCEDURE, PASS:: gatherEF => gatherEFRad - PROCEDURE, PASS:: gatherMF => gatherMFRad PROCEDURE, PASS:: getNodes => getNodesRad PROCEDURE, PASS:: phy2log => phy2logRad PROCEDURE, PASS:: nextElement => nextElementRad - END TYPE meshVol1DRadSegm + END TYPE meshCell1DRadSegm CONTAINS !NODE FUNCTIONS @@ -195,17 +188,18 @@ MODULE moduleMesh1DRad !VOLUME FUNCTIONS !SEGMENT FUNCTIONS !Init segment element - SUBROUTINE initVol1DRadSegm(self, n, p, nodes) + SUBROUTINE initCell1DRadSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(out):: self + CLASS(meshCell1DRadSegm), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2 self%n = n + self%nNodes = SIZE(p) self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj !Get element coordinates @@ -223,22 +217,22 @@ MODULE moduleMesh1DRad ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVol1DRadSegm + END SUBROUTINE initCell1DRadSegm !Calculates a random position in 1D volume FUNCTION randPos1DRadSeg(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: fPsi(1:2) Xi(1) = random(-1.D0, 1.D0) Xi(2:3) = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%r) END FUNCTION randPos1DRadSeg @@ -247,7 +241,7 @@ MODULE moduleMesh1DRad PURE SUBROUTINE areaRad(self) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(inout):: self + CLASS(meshCell1DRadSegm), INTENT(inout):: self REAL(8):: l !element length REAL(8):: fPsi(1:2), fPsi_node(1:2) REAL(8):: r @@ -258,7 +252,7 @@ MODULE moduleMesh1DRad self%arNodes = 0.D0 !1 point Gauss integral Xi = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) detJ = self%detJac(Xi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) @@ -266,37 +260,40 @@ MODULE moduleMesh1DRad self%volume = r*l !Computes volume per node Xi = (/-5.D-1, 0.D0, 0.D0/) - CALL self%fPsi(Xi, fPsi_node) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(1) = fPsi(1)*r*l Xi = (/ 5.D-1, 0.D0, 0.D0/) - CALL self%fPsi(Xi, fPsi_node) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(2) = fPsi(2)*r*l END SUBROUTINE areaRad !Computes element functions at point Xi - PURE SUBROUTINE fPsiRad(xi, fPsi) + PURE FUNCTION fPsiRad(self, xi) RESULT(fPsi) IMPLICIT NONE + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + REAL(8):: fPsi(1:self%nNodes) fPsi(1) = 1.D0 - xi(1) fPsi(2) = 1.D0 + xi(1) + fPsi = fPsi * 5.D-1 - END SUBROUTINE fPsiRad + END FUNCTION fPsiRad !Computes element derivative shape function at Xi - PURE FUNCTION dPsiRad(xi) RESULT(dPsi) + PURE FUNCTION dPsiRad(self, xi) RESULT(dPsi) IMPLICIT NONE + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:1, 1:2)) + dPsi = 0.D0 dPsi(1, 1) = -5.D-1 dPsi(1, 2) = 5.D-1 @@ -307,8 +304,8 @@ MODULE moduleMesh1DRad PURE SUBROUTINE partialDerRad(self, dPsi, dx) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx dx(1) = DOT_PRODUCT(dPsi(1,:), self%r) @@ -320,15 +317,14 @@ MODULE moduleMesh1DRad USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) REAL(8):: Xi(1:3) - REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ(1), detJ + REAL(8):: dPsi(1:3, 1:2) + REAL(8):: invJ(1:3,1:3), detJ REAL(8):: r, fPsi(1:2) INTEGER:: l - ALLOCATE(localK(1:2, 1:2)) localK = 0.D0 Xi = 0.D0 DO l = 1, 3 @@ -336,7 +332,7 @@ MODULE moduleMesh1DRad dPsi = self%dPsi(Xi) detJ = self%detJac(Xi, dPsi) invJ = self%invJac(Xi, dPsi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r = DOT_PRODUCT(fPsi, self%r) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & @@ -352,22 +348,21 @@ MODULE moduleMesh1DRad USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: fPsi(1:2) REAL(8):: detJ, f, r REAL(8):: Xi(1:3) INTEGER:: l - ALLOCATE(localF(1:2)) localF = 0.D0 Xi = 0.D0 DO l = 1, 3 Xi(1) = corSeg(l) detJ = self%detJac(Xi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r = DOT_PRODUCT(fPsi, self%r) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*r*wSeg(l)*detJ @@ -376,6 +371,40 @@ MODULE moduleMesh1DRad END FUNCTION elemFRad + PURE FUNCTION gatherEFRad(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:2) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEFRad + + PURE FUNCTION gatherMFRad(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:2,1:3) + + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFRad + PURE FUNCTION insideRad(xi) RESULT(ins) IMPLICIT NONE @@ -387,58 +416,13 @@ MODULE moduleMesh1DRad END FUNCTION insideRad - !Gathers EF at position Xi - PURE FUNCTION gatherEFRad(self, xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1, 1:2) - REAL(8):: phi(1:2) - REAL(8):: EF(1:3) - REAL(8):: invJ - - phi = (/ self%n1%emData%phi, & - self%n2%emData%phi /) - - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ - EF(2) = 0.D0 - EF(3) = 0.D0 - - END FUNCTION gatherEFRad - - PURE FUNCTION gatherMFRad(self, xi) RESULT(MF) - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:2) - REAL(8):: MF_Nodes(1:2, 1:3) - REAL(8):: MF(1:3) - REAL(8):: invJ - - MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), & - self%n2%emData%B(1) /) - MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), & - self%n2%emData%B(2) /) - MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & - self%n2%emData%B(3) /) - - CALL self%fPsi(xi, fPsi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFRad - !Get nodes from 1D volume PURE FUNCTION getNodesRad(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:2)) n = (/ self%n1%n, self%n2%n /) END FUNCTION getNodesRad @@ -446,7 +430,7 @@ MODULE moduleMesh1DRad PURE FUNCTION phy2logRad(self, r) RESULT(xN) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: xN(1:3) @@ -459,7 +443,7 @@ MODULE moduleMesh1DRad SUBROUTINE nextElementRad(self, xi, nextElement) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement @@ -479,10 +463,10 @@ MODULE moduleMesh1DRad PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol1DRad), INTENT(in):: self + CLASS(meshCell1DRad), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dJ REAL(8):: dx(1) @@ -503,12 +487,12 @@ MODULE moduleMesh1DRad PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol1DRad), INTENT(in):: self + CLASS(meshCell1DRad), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dx(1) - REAL(8):: invJ + REAL(8):: invJ(1:3,1:3) IF (PRESENT(dPsi_in)) THEN dPsi = dPsi_in @@ -518,8 +502,11 @@ MODULE moduleMesh1DRad END IF + invJ = 0.D0 + CALL self%partialDer(dPsi, dx) - invJ = 1.D0/dx(1) + + invJ(1,1) = 1.D0/dx(1) END FUNCTION invJ1DRad @@ -529,11 +516,11 @@ MODULE moduleMesh1DRad CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et - DO e = 1, self%numVols + DO e = 1, self%numCells !Connect Vol-Vol - DO et = 1, self%numVols + DO et = 1, self%numCells IF (e /= et) THEN - CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj) END IF @@ -543,7 +530,7 @@ MODULE moduleMesh1DRad TYPE IS(meshParticles) !Connect Vol-Edge DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj) END DO @@ -556,13 +543,13 @@ MODULE moduleMesh1DRad SUBROUTINE connectVolVol(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVol1DRadSegm) + TYPE IS(meshCell1DRadSegm) SELECT TYPE(elemB) - TYPE IS(meshVol1DRadSegm) + TYPE IS(meshCell1DRadSegm) CALL connectSegmSegm(elemA, elemB) END SELECT @@ -574,8 +561,8 @@ MODULE moduleMesh1DRad SUBROUTINE connectSegmSegm(elemA, elemB) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB + CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemA + CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemB IF (.NOT. ASSOCIATED(elemA%e1) .AND. & elemA%n2%n == elemB%n1%n) THEN @@ -597,11 +584,11 @@ MODULE moduleMesh1DRad SUBROUTINE connectVolEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemA CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS (meshVol1DRadSegm) + TYPE IS (meshCell1DRadSegm) SELECT TYPE(elemB) CLASS IS(meshEdge1DRad) CALL connectSegmEdge(elemA, elemB) @@ -615,7 +602,7 @@ MODULE moduleMesh1DRad SUBROUTINE connectSegmEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA + CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemA CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB IF (.NOT. ASSOCIATED(elemA%e1) .AND. & diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index c57cecc..13c5901 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -37,26 +37,19 @@ MODULE moduleMesh2DCart END TYPE meshEdge2DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCart + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCart CONTAINS PROCEDURE, PASS:: detJac => detJ2DCart PROCEDURE, PASS:: invJac => invJ2DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer - END TYPE meshVol2DCart + END TYPE meshCell2DCart ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi) - REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - PURE SUBROUTINE partialDer_interface(self, dPsi, dx, dy) - IMPORT meshVol2DCart - CLASS(meshVol2DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + IMPORT meshCell2DCart + CLASS(meshCell2DCart), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy END SUBROUTINE partialDer_interface @@ -64,7 +57,7 @@ MODULE moduleMesh2DCart END INTERFACE !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartQuad + TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartQuad !Element coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0 !Connectivity to nodes @@ -73,27 +66,27 @@ MODULE moduleMesh2DCart CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolQuad2DCart - PROCEDURE, PASS:: randPos => randPosVolQuad + PROCEDURE, PASS:: init => initCellQuad2DCart + PROCEDURE, PASS:: randPos => randPosCellQuad PROCEDURE, PASS:: area => areaQuad - PROCEDURE, NOPASS:: fPsi => fPsiQuad - PROCEDURE, NOPASS:: dPsi => dPsiQuad - PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2 - PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, PASS:: fPsi => fPsiQuad + PROCEDURE, PASS:: dPsi => dPsiQuad + PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiQuadXi1 + PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiQuadXi2 + PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad PROCEDURE, PASS:: elemK => elemKQuad PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, PASS:: gatherElectricField => gatherEFQuad + PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: gatherEF => gatherEFQuad - PROCEDURE, PASS:: gatherMF => gatherMFQuad PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: nextElement => nextElementQuad - END TYPE meshVol2DCartQuad + END TYPE meshCell2DCartQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartTria + TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartTria !Element coordinates REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0 !Connectivity to nodes @@ -103,24 +96,24 @@ MODULE moduleMesh2DCart REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolTria2DCart - PROCEDURE, PASS:: randPos => randPosVolTria + PROCEDURE, PASS:: init => initCellTria2DCart + PROCEDURE, PASS:: randPos => randPosCellTria PROCEDURE, PASS:: area => areaTria - PROCEDURE, NOPASS:: fPsi => fPsiTria - PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: fPsi => fPsiTria + PROCEDURE, PASS:: dPsi => dPsiTria PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 PROCEDURE, PASS:: partialDer => partialDerTria PROCEDURE, PASS:: elemK => elemKTria PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, PASS:: gatherElectricField => gatherEFTria + PROCEDURE, PASS:: gatherMagneticField => gatherMFTria PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: gatherEF => gatherEFTria - PROCEDURE, PASS:: gatherMF => gatherMFTria PROCEDURE, PASS:: getNodes => getNodesTria PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: nextElement => nextElementTria - END TYPE meshVol2DCartTria + END TYPE meshCell2DCartTria CONTAINS !NODE FUNCTIONS @@ -204,26 +197,26 @@ MODULE moduleMesh2DCart END SUBROUTINE initEdge2DCart !Random position in quadrilateral volume - FUNCTION randPosVolQuad(self) RESULT(r) + FUNCTION randPosCellQuad(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: fPsi(1:4) Xi(1) = random(-1.D0, 1.D0) Xi(2) = random(-1.D0, 1.D0) Xi(3) = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) r(3) = 0.D0 - END FUNCTION randposVolQuad + END FUNCTION randposCellQuad !Get nodes from edge PURE FUNCTION getNodes2DCart(self) RESULT(n) @@ -232,7 +225,6 @@ MODULE moduleMesh2DCart CLASS(meshEdge2DCart), INTENT(in):: self INTEGER, ALLOCATABLE:: n(:) - ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) END FUNCTION getNodes2DCart @@ -277,17 +269,18 @@ MODULE moduleMesh2DCart !VOLUME FUNCTIONS !QUAD FUNCTIONS !Inits quadrilateral element - SUBROUTINE initVolQuad2DCart(self, n, p, nodes) + SUBROUTINE initCellQuad2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(out):: self + CLASS(meshCell2DCartQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 self%n = n + self%nNodes = SIZE(p) self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -312,13 +305,13 @@ MODULE moduleMesh2DCart ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolQuad2DCart + END SUBROUTINE initCellQuad2DCart !Computes element area PURE SUBROUTINE areaQuad(self) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(inout):: self + CLASS(meshCell2DCartQuad), INTENT(inout):: self REAL(8):: Xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:4) @@ -328,18 +321,19 @@ MODULE moduleMesh2DCart !2D 1 point Gauss Quad Integral Xi = 0.D0 detJ = self%detJac(Xi)*4.D0 !4 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) self%volume = detJ self%arNodes = fPsi*detJ END SUBROUTINE areaQuad !Computes element functions in point Xi - PURE SUBROUTINE fPsiQuad(Xi, fPsi) + PURE FUNCTION fPsiQuad(self, Xi) RESULT(fPsi) IMPLICIT NONE + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + REAL(8):: fPsi(1:self%nNodes) fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2)) fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2)) @@ -348,16 +342,17 @@ MODULE moduleMesh2DCart fPsi = fPsi*0.25D0 - END SUBROUTINE fPsiQuad + END FUNCTION fPsiQuad !Derivative element function at coordinates Xi - PURE FUNCTION dPsiQuad(Xi) RESULT(dPsi) + PURE FUNCTION dPsiQuad(self, Xi) RESULT(dPsi) IMPLICIT NONE + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:2,1:4)) + dPsi = 0.D0 dPsi(1,:) = dPsiQuadXi1(Xi(2)) dPsi(2,:) = dPsiQuadXi2(Xi(1)) @@ -397,8 +392,8 @@ MODULE moduleMesh2DCart PURE SUBROUTINE partialDerQuad(self, dPsi, dx, dy) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) @@ -412,14 +407,13 @@ MODULE moduleMesh2DCart PURE FUNCTION elemKQuad(self) RESULT(localK) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) REAL(8):: Xi(1:3) - REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: invJ(1:2,1:2), detJ + REAL(8):: fPsi(1:4), dPsi(1:3,1:4) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m - ALLOCATE(localK(1:4, 1:4)) localK=0.D0 Xi=0.D0 !Start 2D Gauss Quad Integral @@ -429,7 +423,7 @@ MODULE moduleMesh2DCart DO m = 1, 3 Xi(1) = corQuad(m) dPsi(2,:) = self%dPsiXi2(Xi(1)) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) detJ = self%detJac(Xi,dPsi) invJ = self%invJac(Xi,dPsi) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ @@ -443,15 +437,14 @@ MODULE moduleMesh2DCart PURE FUNCTION elemFQuad(self, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: detJ, f INTEGER:: l, m - ALLOCATE(localF(1:4)) localF = 0.D0 Xi = 0.D0 DO l=1, 3 @@ -459,7 +452,7 @@ MODULE moduleMesh2DCart DO m = 1, 3 Xi(2) = corQuad(m) detJ = self%detJac(Xi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -468,6 +461,48 @@ MODULE moduleMesh2DCart END FUNCTION elemFQuad + PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:4) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEFQuad + + PURE FUNCTION gatherMFQuad(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:4,1:3) + + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFQuad + !Checks if a particle is inside a quad element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -480,97 +515,42 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Gathers the electric field at position Xi - PURE FUNCTION gatherEFQuad(self,Xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: dPsi(1:2,1:4) - REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates - REAL(8):: invJ(1:2,1:2), detJ - REAL(8):: phi(1:4) - REAL(8):: EF(1:3) - - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi, & - self%n4%emData%phi /) - - dPsi = self%dPsi(Xi) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = 0.D0 - - END FUNCTION gatherEFQuad - - PURE FUNCTION gatherMFQuad(self,Xi) RESULT(MF) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: fPsi(1:4) - REAL(8):: MF_Nodes(1:4,1:3) - REAL(8):: MF(1:3) - - MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & - self%n2%emData%B(1), & - self%n3%emData%B(1), & - self%n4%emData%B(1) /) - MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & - self%n2%emData%B(2), & - self%n3%emData%B(2), & - self%n4%emData%B(2) /) - MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & - self%n2%emData%B(3), & - self%n3%emData%B(3), & - self%n4%emData%B(3) /) - - CALL self%fPsi(Xi, fPsi) - MF = MATMUL(fPsi(:), MF_Nodes) - - END FUNCTION gatherMFQuad - !Gets nodes from quadrilateral element PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:4)) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) END FUNCTION getNodesQuad !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(XiN) + PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: XiN(1:3) + REAL(8):: Xi(1:3) REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2) - REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: conv !Iterative newton method to transform coordinates conv=1.D0 XiO=0.D0 - DO WHILE(conv>1.D-4) + DO WHILE(conv>1.D-3) dPsi = self%dPsi(XiO) invJ = self%invJac(XiO, dPsi) - CALL self%fPsi(XiO, fPsi) + fPsi = self%fPsi(XiO) f(1) = DOT_PRODUCT(fPsi,self%x)-r(1) f(2) = DOT_PRODUCT(fPsi,self%y)-r(2) detJ = self%detJac(XiO,dPsi) - XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(XiN-XiO),1) - XiO=XiN + Xi(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ + conv=MAXVAL(DABS(Xi-XiO),1) + XiO=Xi END DO @@ -580,7 +560,7 @@ MODULE moduleMesh2DCart SUBROUTINE nextElementQuad(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement REAL(8):: XiArray(1:4) @@ -605,11 +585,11 @@ MODULE moduleMesh2DCart !TRIA ELEMENT !Init tria element - SUBROUTINE initVolTria2DCart(self, n, p, nodes) + SUBROUTINE initCellTria2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(out):: self + CLASS(meshCell2DCartTria), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) @@ -618,6 +598,9 @@ MODULE moduleMesh2DCart !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 @@ -639,14 +622,14 @@ MODULE moduleMesh2DCart ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolTria2DCart + END SUBROUTINE initCellTria2DCart !Random position in quadrilateral volume - FUNCTION randPosVolTria(self) RESULT(r) + FUNCTION randPosCellTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self + CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:3) @@ -655,19 +638,19 @@ MODULE moduleMesh2DCart Xi(2) = random( 0.D0, 1.D0 - Xi(1)) Xi(3) = 0.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) r(3) = 0.D0 - END FUNCTION randposVolTria + END FUNCTION randposCellTria !Calculates area for triangular element PURE SUBROUTINE areaTria(self) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(inout):: self + CLASS(meshCell2DCartTria), INTENT(inout):: self REAL(8):: Xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:3) @@ -677,33 +660,35 @@ MODULE moduleMesh2DCart !2D 1 point Gauss Quad Integral Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) detJ = self%detJac(Xi)/2.D0 - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) self%volume = detJ self%arNodes = fPsi*detJ END SUBROUTINE areaTria !Shape functions for triangular element - PURE SUBROUTINE fPsiTria(Xi, fPsi) + PURE FUNCTION fPsiTria(self, Xi) RESULT(fPsi) IMPLICIT NONE + CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + REAL(8):: fPsi(1:self%nNodes) fPsi(1) = 1.D0 - Xi(1) - Xi(2) fPsi(2) = Xi(1) fPsi(3) = Xi(2) - END SUBROUTINE fPsiTria + END FUNCTION fPsiTria !Derivative element function at coordinates Xi - PURE FUNCTION dPsiTria(Xi) RESULT(dPsi) + PURE FUNCTION dPsiTria(self, Xi) RESULT(dPsi) IMPLICIT NONE + CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:2,1:3)) + dPsi = 0.D0 dPsi(1,:) = dPsiTriaXi1(Xi(2)) dPsi(2,:) = dPsiTriaXi2(Xi(1)) @@ -739,8 +724,8 @@ MODULE moduleMesh2DCart PURE SUBROUTINE partialDerTria(self, dPsi, dx, dy) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) @@ -754,14 +739,13 @@ MODULE moduleMesh2DCart PURE FUNCTION elemKTria(self) RESULT(localK) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) REAL(8):: Xi(1:3) - REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: invJ(1:2,1:2), detJ + REAL(8):: fPsi(1:3), dPsi(1:3,1:3) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l - ALLOCATE(localK(1:4, 1:4)) localK=0.D0 Xi=0.D0 !Start 2D Gauss Quad Integral @@ -771,7 +755,7 @@ MODULE moduleMesh2DCart dPsi = self%dPsi(Xi) detJ = self%detJac(Xi,dPsi) invJ = self%invJac(Xi,dPsi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -782,15 +766,14 @@ MODULE moduleMesh2DCart PURE FUNCTION elemFTria(self, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: fPsi(1:3) REAL(8):: Xi(1:3) REAL(8):: detJ, f INTEGER:: l - ALLOCATE(localF(1:3)) localF = 0.D0 Xi = 0.D0 !Start 2D Gauss Quad Integral @@ -798,7 +781,7 @@ MODULE moduleMesh2DCart Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) detJ = self%detJac(Xi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wTria(l)*detJ @@ -806,6 +789,44 @@ MODULE moduleMesh2DCart END FUNCTION elemFTria + PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:3) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEFTria + + PURE FUNCTION gatherMFTria(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:3,1:3) + + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFTria + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE @@ -818,64 +839,13 @@ MODULE moduleMesh2DCart END FUNCTION insideTria - !Gathers the electric field at position Xi - PURE FUNCTION gatherEFTria(self,Xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: dPsi(1:2,1:3) - REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates - REAL(8):: invJ(1:2,1:2), detJ - REAL(8):: phi(1:3) - REAL(8):: EF(1:3) - - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi /) - - dPsi = self%dPsi(Xi) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = 0.D0 - - END FUNCTION gatherEFTria - - PURE FUNCTION gatherMFTria(self,Xi) RESULT(MF) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: fPsi(1:3) - REAL(8):: MF_Nodes(1:3,1:3) - REAL(8):: MF(1:3) - - MF_Nodes(1:3,1) = (/self%n1%emData%B(1), & - self%n2%emData%B(1), & - self%n3%emData%B(1) /) - MF_Nodes(1:3,2) = (/self%n1%emData%B(2), & - self%n2%emData%B(2), & - self%n3%emData%B(2) /) - MF_Nodes(1:3,3) = (/self%n1%emData%B(3), & - self%n2%emData%B(3), & - self%n3%emData%B(3) /) - - CALL self%fPsi(Xi, fPsi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFTria - !Gets node indexes from triangular element PURE FUNCTION getNodesTria(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:3)) n = (/self%n1%n, self%n2%n, self%n3%n /) END FUNCTION getNodesTria @@ -884,27 +854,27 @@ MODULE moduleMesh2DCart PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self + CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: invJ(1:2,1:2), detJ - REAL(8):: deltaR(1:2) - REAL(8):: dPsi(1:2,1:3) + REAL(8):: invJ(1:3,1:3), detJ + REAL(8):: deltaR(1:3) + REAL(8):: dPsi(1:3,1:3) !Direct method to convert coordinates - Xi = 0.D0 !Irrelevant, required for input - deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) - dPsi = self%dPsi(Xi) - invJ = self%invJac(Xi, dPsi) - detJ = self%detJac(Xi, dPsi) - Xi(1:2) = MATMUL(invJ,deltaR)/detJ + Xi = 0.D0 + deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) + dPsi = self%dPsi(Xi) + invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, dPsi) + Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria SUBROUTINE nextElementTria(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self + CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement REAL(8):: XiArray(1:3) @@ -929,10 +899,10 @@ MODULE moduleMesh2DCart PURE FUNCTION detJ2DCart(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol2DCart), INTENT(in):: self + CLASS(meshCell2DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dJ REAL(8):: dx(1:2), dy(1:2) @@ -953,12 +923,12 @@ MODULE moduleMesh2DCart PURE FUNCTION invJ2DCart(self,Xi,dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol2DCart), INTENT(in):: self + CLASS(meshCell2DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dx(1:2), dy(1:2) - REAL(8):: invJ(1:2,1:2) + REAL(8):: invJ(1:3,1:3) IF(PRESENT(dPsi_in)) THEN dPsi=dPsi_in @@ -968,9 +938,12 @@ MODULE moduleMesh2DCart END IF + invJ = 0.D0 + CALL self%partialDer(dPsi, dx, dy) - invJ(1,:) = (/ dy(2), -dx(2) /) - invJ(2,:) = (/ -dy(1), dx(1) /) + + invJ(1,1:2) = (/ dy(2), -dx(2) /) + invJ(2,1:2) = (/ -dy(1), dx(1) /) END FUNCTION invJ2DCart @@ -980,11 +953,11 @@ MODULE moduleMesh2DCart CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et - DO e = 1, self%numVols - !Connect Vol-Vol - DO et = 1, self%numVols + DO e = 1, self%numCells + !Connect Cell-Cell + DO et = 1, self%numCells IF (e /= et) THEN - CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj) END IF @@ -992,9 +965,9 @@ MODULE moduleMesh2DCart SELECT TYPE(self) TYPE IS(meshParticles) - !Connect Vol-Edge + !Connect Cell-Edge DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj) END DO @@ -1005,34 +978,34 @@ MODULE moduleMesh2DCart END SUBROUTINE connectMesh2DCart !Selects type of elements to build connection - SUBROUTINE connectVolVol(elemA, elemB) + SUBROUTINE connectCellCell(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVol2DCartQuad) + TYPE IS(meshCell2DCartQuad) !Element A is a quadrilateral SELECT TYPE(elemB) - TYPE IS(meshVol2DCartQuad) + TYPE IS(meshCell2DCartQuad) !Element B is a quadrilateral CALL connectQuadQuad(elemA, elemB) - TYPE IS(meshVol2DCartTria) + TYPE IS(meshCell2DCartTria) !Element B is a triangle CALL connectQuadTria(elemA, elemB) END SELECT - TYPE IS(meshVol2DCartTria) + TYPE IS(meshCell2DCartTria) !Element A is a Triangle SELECT TYPE(elemB) - TYPE IS(meshVol2DCartQuad) + TYPE IS(meshCell2DCartQuad) !Element B is a quadrilateral CALL connectQuadTria(elemB, elemA) - TYPE IS(meshVol2DCartTria) + TYPE IS(meshCell2DCartTria) !Element B is a triangle CALL connectTriaTria(elemA, elemB) @@ -1040,22 +1013,22 @@ MODULE moduleMesh2DCart END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell - SUBROUTINE connectVolEdge(elemA, elemB) + SUBROUTINE connectCellEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemA CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemB) CLASS IS(meshEdge2DCart) SELECT TYPE(elemA) - TYPE IS(meshVol2DCartQuad) + TYPE IS(meshCell2DCartQuad) !Element A is a quadrilateral CALL connectQuadEdge(elemA, elemB) - TYPE IS(meshVol2DCartTria) + TYPE IS(meshCell2DCartTria) !Element A is a triangle CALL connectTriaEdge(elemA, elemB) @@ -1063,13 +1036,13 @@ MODULE moduleMesh2DCart END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge SUBROUTINE connectQuadQuad(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1) .AND. & @@ -1112,8 +1085,8 @@ MODULE moduleMesh2DCart SUBROUTINE connectQuadTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -1204,8 +1177,8 @@ MODULE moduleMesh2DCart SUBROUTINE connectTriaTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -1277,7 +1250,7 @@ MODULE moduleMesh2DCart SUBROUTINE connectQuadEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB !Check direction 1 @@ -1361,7 +1334,7 @@ MODULE moduleMesh2DCart SUBROUTINE connectTriaEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemA CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB !Check direction 1 diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 2869cb3..9032d61 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -1,4 +1,4 @@ -!moduleMesh2DCyl: 2D axial symmetric extension of generic mesh from GMSH format. +!moduleMesh2DCyl: 2D aXial symmetric extension of generic mesh from GMSH format. ! x == z ! y == r ! z == theta (unused) @@ -11,8 +11,8 @@ MODULE moduleMesh2DCyl REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) - REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) - REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) + REAL(8), PARAMETER:: Xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) + REAL(8), PARAMETER:: Xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /) TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl @@ -37,26 +37,19 @@ MODULE moduleMesh2DCyl END TYPE meshEdge2DCyl - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCyl CONTAINS PROCEDURE, PASS:: detJac => detJ2DCyl PROCEDURE, PASS:: invJac => invJ2DCyl - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer - END TYPE meshVol2DCyl + END TYPE meshCell2DCyl ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr) - IMPORT meshVol2DCyl - CLASS(meshVol2DCyl), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + IMPORT meshCell2DCyl + CLASS(meshCell2DCyl), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr END SUBROUTINE partialDer_interface @@ -64,7 +57,7 @@ MODULE moduleMesh2DCyl END INTERFACE !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylQuad + TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylQuad !Element coordinates REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -74,27 +67,27 @@ MODULE moduleMesh2DCyl REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolQuad2DCyl + PROCEDURE, PASS:: init => initCellQuad2DCyl PROCEDURE, PASS:: randPos => randPosVolQuad PROCEDURE, PASS:: area => areaQuad - PROCEDURE, NOPASS:: fPsi => fPsiQuad - PROCEDURE, NOPASS:: dPsi => dPsiQuad - PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2 - PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, PASS:: fPsi => fPsiQuad + PROCEDURE, PASS:: dPsi => dPsiQuad + PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiQuadXi1 + PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiQuadXi2 + PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad PROCEDURE, PASS:: elemK => elemKQuad PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, PASS:: gatherElectricField => gatherEFQuad + PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: gatherEF => gatherEFQuad - PROCEDURE, PASS:: gatherMF => gatherMFQuad PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: nextElement => nextElementQuad - END TYPE meshVol2DCylQuad + END TYPE meshCell2DCylQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylTria + TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylTria !Element coordinates REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 !Connectivity to nodes @@ -104,24 +97,24 @@ MODULE moduleMesh2DCyl REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolTria2DCyl + PROCEDURE, PASS:: init => initCellTria2DCyl PROCEDURE, PASS:: randPos => randPosVolTria PROCEDURE, PASS:: area => areaTria - PROCEDURE, NOPASS:: fPsi => fPsiTria - PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: fPsi => fPsiTria + PROCEDURE, PASS:: dPsi => dPsiTria PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 - PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, PASS, PRIVATE:: partialDer => partialDerTria PROCEDURE, PASS:: elemK => elemKTria PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, PASS:: gatherElectricField => gatherEFTria + PROCEDURE, PASS:: gatherMagneticField => gatherMFTria PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: gatherEF => gatherEFTria - PROCEDURE, PASS:: gatherMF => gatherMFTria PROCEDURE, PASS:: getNodes => getNodesTria PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: nextElement => nextElementTria - END TYPE meshVol2DCylTria + END TYPE meshCell2DCylTria CONTAINS !NODE FUNCTIONS @@ -265,17 +258,18 @@ MODULE moduleMesh2DCyl !VOLUME FUNCTIONS !QUAD FUNCTIONS !Inits quadrilateral element - SUBROUTINE initVolQuad2DCyl(self, n, p, nodes) + SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(out):: self + CLASS(meshCell2DCylQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 self%n = n + self%nNodes = SIZE(p) self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -300,104 +294,106 @@ MODULE moduleMesh2DCyl ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolQuad2DCyl + END SUBROUTINE initCellQuad2DCyl !Computes element area PURE SUBROUTINE areaQuad(self) USE moduleConstParam, ONLY: PI8 IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(inout):: self - REAL(8):: r, xi(1:3) + CLASS(meshCell2DCylQuad), INTENT(inout):: self + REAL(8):: r, Xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:4), fPsi_node(1:4) self%volume = 0.D0 self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral - xi = 0.D0 - detJ = self%detJac(xi)*PI8 !4*2*pi - CALL self%fPsi(xi, fPsi) + Xi = 0.D0 + detJ = self%detJac(Xi)*PI8 !4*2*pi + fPsi = self%fPsi(Xi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ !Computes volume per node - xi = (/-5.D-1, -5.D-1, 0.D0/) - CALL self%fPsi(xi, fPsi_node) + Xi = (/-5.D-1, -5.D-1, 0.D0/) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(1) = fPsi(1)*r*detJ - xi = (/ 5.D-1, -5.D-1, 0.D0/) - CALL self%fPsi(xi, fPsi_node) + Xi = (/ 5.D-1, -5.D-1, 0.D0/) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(2) = fPsi(2)*r*detJ - xi = (/ 5.D-1, 5.D-1, 0.D0/) - CALL self%fPsi(xi, fPsi_node) + Xi = (/ 5.D-1, 5.D-1, 0.D0/) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(3) = fPsi(3)*r*detJ - xi = (/-5.D-1, 5.D-1, 0.D0/) - CALL self%fPsi(xi, fPsi_node) + Xi = (/-5.D-1, 5.D-1, 0.D0/) + fPsi_node = self%fPsi(Xi) r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(4) = fPsi(4)*r*detJ END SUBROUTINE areaQuad - !Computes element functions in point xi - PURE SUBROUTINE fPsiQuad(xi, fPsi) + !Computes element functions in point Xi + PURE FUNCTION fPsiQuad(self, Xi) RESULT(fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: fPsi(1:self%nNodes) - fPsi(1) = (1.D0-xi(1)) * (1.D0-xi(2)) - fPsi(2) = (1.D0+xi(1)) * (1.D0-xi(2)) - fPsi(3) = (1.D0+xi(1)) * (1.D0+xi(2)) - fPsi(4) = (1.D0-xi(1)) * (1.D0+xi(2)) + fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2)) + fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2)) + fPsi(3) = (1.D0+Xi(1)) * (1.D0+Xi(2)) + fPsi(4) = (1.D0-Xi(1)) * (1.D0+Xi(2)) fPsi = fPsi*0.25D0 - END SUBROUTINE fPsiQuad + END FUNCTION fPsiQuad - !Derivative element function at coordinates xi - PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiQuad(self, Xi) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:2,1:4)) + dPsi = 0.D0 - dPsi(1,:) = dPsiQuadXi1(xi(2)) - dPsi(2,:) = dPsiQuadXi2(xi(1)) + dPsi(1,:) = dPsiQuadXi1(Xi(2)) + dPsi(2,:) = dPsiQuadXi2(Xi(1)) END FUNCTION dPsiQuad - !Derivative element function (xi1) - PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) + !Derivative element function (Xi1) + PURE FUNCTION dPsiQuadXi1(Xi2) RESULT(dPsiXi1) IMPLICIT NONE - REAL(8),INTENT(in):: xi2 + REAL(8),INTENT(in):: Xi2 REAL(8):: dPsiXi1(1:4) - dPsiXi1(1) = -(1.D0 - xi2) - dPsiXi1(2) = (1.D0 - xi2) - dPsiXi1(3) = (1.D0 + xi2) - dPsiXi1(4) = -(1.D0 + xi2) + dPsiXi1(1) = -(1.D0 - Xi2) + dPsiXi1(2) = (1.D0 - Xi2) + dPsiXi1(3) = (1.D0 + Xi2) + dPsiXi1(4) = -(1.D0 + Xi2) dPsiXi1 = dPsiXi1*0.25D0 END FUNCTION dPsiQuadXi1 - !Derivative element function (xi2) - PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) + !Derivative element function (Xi2) + PURE FUNCTION dPsiQuadXi2(Xi1) RESULT(dPsiXi2) IMPLICIT NONE - REAL(8),INTENT(in):: xi1 + REAL(8),INTENT(in):: Xi1 REAL(8):: dPsiXi2(1:4) - dPsiXi2(1) = -(1.D0 - xi1) - dPsiXi2(2) = -(1.D0 + xi1) - dPsiXi2(3) = (1.D0 + xi1) - dPsiXi2(4) = (1.D0 - xi1) + dPsiXi2(1) = -(1.D0 - Xi1) + dPsiXi2(2) = -(1.D0 + Xi1) + dPsiXi2(3) = (1.D0 + Xi1) + dPsiXi2(4) = (1.D0 - Xi1) dPsiXi2 = dPsiXi2 * 0.25D0 @@ -407,8 +403,8 @@ MODULE moduleMesh2DCyl PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) @@ -423,16 +419,16 @@ MODULE moduleMesh2DCyl USE moduleRandom IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self + CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:4) - xii(1) = random(-1.D0, 1.D0) - xii(2) = random(-1.D0, 1.D0) - xii(3) = 0.D0 + Xi(1) = random(-1.D0, 1.D0) + Xi(2) = random(-1.D0, 1.D0) + Xi(3) = 0.D0 - CALL self%fPsi(xii, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) @@ -445,26 +441,25 @@ MODULE moduleMesh2DCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: r, xi(1:3) - REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) + REAL(8):: r, Xi(1:3) + REAL(8):: fPsi(1:4), dPsi(1:3,1:4) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m - ALLOCATE(localK(1:4, 1:4)) localK=0.D0 - xi=0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 - xi(2) = corQuad(l) - dPsi(1,:) = self%dPsiXi1(xi(2)) + Xi(2) = corQuad(l) + dPsi(1,:) = self%dPsiXi1(Xi(2)) DO m = 1, 3 - xi(1) = corQuad(m) - dPsi(2,:) = self%dPsiXi2(xi(1)) - CALL self%fPsi(xi, fPsi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) + Xi(1) = corQuad(m) + dPsi(2,:) = self%dPsiXi2(Xi(1)) + fPsi = self%fPsi(Xi) + detJ = self%detJac(Xi,dPsi) + invJ = self%invJac(Xi,dPsi) r = DOT_PRODUCT(fPsi,self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & MATMUL(invJ,dPsi))* & @@ -481,23 +476,22 @@ MODULE moduleMesh2DCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: r, xi(1:3) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) + REAL(8):: r, Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: detJ, f INTEGER:: l, m - ALLOCATE(localF(1:4)) localF = 0.D0 - xi = 0.D0 + Xi = 0.D0 DO l=1, 3 - xi(1) = corQuad(l) + Xi(1) = corQuad(l) DO m = 1, 3 - xi(2) = corQuad(m) - detJ = self%detJac(xi) - CALL self%fPsi(xi, fPsi) + Xi(2) = corQuad(m) + detJ = self%detJac(Xi) + fPsi = self%fPsi(Xi) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -508,93 +502,80 @@ MODULE moduleMesh2DCyl END FUNCTION elemFQuad - !Checks if a particle is inside a quad element - PURE FUNCTION insideQuad(xi) RESULT(ins) + PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array) IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = (xi(1) >= -1.D0 .AND. xi(1) <= 1.D0) .AND. & - (xi(2) >= -1.D0 .AND. xi(2) <= 1.D0) - - END FUNCTION insideQuad - - !Gathers the electric field at position xi - PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1:2,1:4) - REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) REAL(8):: phi(1:4) - REAL(8):: EF(1:3) - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi, & - self%n4%emData%phi /) + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = 0.D0 + array = -self%gatherDF(Xi, phi) END FUNCTION gatherEFQuad - PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF) + PURE FUNCTION gatherMFQuad(self, Xi) RESULT(array) IMPLICIT NONE + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:4,1:3) - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:4) - REAL(8):: MF_Nodes(1:4,1:3) - REAL(8):: MF(1:3) + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) - MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & - self%n2%emData%B(1), & - self%n3%emData%B(1), & - self%n4%emData%B(1) /) - MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & - self%n2%emData%B(2), & - self%n3%emData%B(2), & - self%n4%emData%B(2) /) - MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & - self%n2%emData%B(3), & - self%n3%emData%B(3), & - self%n4%emData%B(3) /) + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) - CALL self%fPsi(xi, fPsi) - MF = MATMUL(fPsi(:), MF_Nodes) + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) END FUNCTION gatherMFQuad + !Checks if a particle is inside a quad element + PURE FUNCTION insideQuad(Xi) RESULT(ins) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + LOGICAL:: ins + + ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. & + (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0) + + END FUNCTION insideQuad + !Gets nodes from quadrilateral element PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell2DCylQuad), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:4)) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) END FUNCTION getNodesQuad !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(XiN) + PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self + CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: XiN(1:3) + REAL(8):: Xi(1:3) REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2) - REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: conv !Iterative newton method to transform coordinates @@ -602,32 +583,33 @@ MODULE moduleMesh2DCyl XiO=0.D0 DO WHILE(conv>1.D-3) - CALL self%fPsi(XiO, fPsi) - f = (/ DOT_PRODUCT(fPsi,self%z)-r(1), & - DOT_PRODUCT(fPsi,self%r)-r(2) /) dPsi = self%dPsi(XiO) invJ = self%invJac(XiO, dPsi) - detJ = self%detJac(XiO,dPsi) - XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(XiN-XiO),1) - XiO=XiN + detJ = self%detJac(XiO, dPsi) + fPsi = self%fPsi(XiO) + f = (/ DOT_PRODUCT(fPsi,self%z), & + DOT_PRODUCT(fPsi,self%r) /) + f = f - r(1:2) + Xi(1:2) = XiO(1:2) - MATMUL(invJ, f)/detJ + conv = MAXVAL(DABS(Xi-XiO),1) + XiO = Xi END DO END FUNCTION phy2logQuad - !Gets the next element for a logical position xi - SUBROUTINE nextElementQuad(self, xi, nextElement) + !Gets the next element for a logical position Xi + SUBROUTINE nextElementQuad(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + REAL(8):: XiArray(1:4) INTEGER:: nextInt - xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /) - nextInt = MAXLOC(xiArray,1) + XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) + nextInt = MAXLOC(XiArray,1) !Selects the higher value of directions and searches in that direction NULLIFY(nextElement) SELECT CASE (nextInt) @@ -645,11 +627,11 @@ MODULE moduleMesh2DCyl !TRIA ELEMENT !Init tria element - SUBROUTINE initVolTria2DCyl(self, n, p, nodes) + SUBROUTINE initCellTria2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(out):: self + CLASS(meshCell2DCylTria), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) @@ -658,6 +640,9 @@ MODULE moduleMesh2DCyl !Assign node index self%n = n + !Assign nomber of nodes to cell + self%nNodes = SIZE(p) + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj @@ -679,23 +664,23 @@ MODULE moduleMesh2DCyl ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolTria2DCyl + END SUBROUTINE initCellTria2DCyl !Random position in quadrilateral volume FUNCTION randPosVolTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self + CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:3) - xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0 - xii(1)) - xii(3) = 0.D0 + Xi(1) = random( 0.D0, 1.D0) + Xi(2) = random( 0.D0, 1.D0 - Xi(1)) + Xi(3) = 0.D0 - CALL self%fPsi(xii, fPsi) + fPsi = self%fPsi(Xi) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) @@ -708,17 +693,17 @@ MODULE moduleMesh2DCyl USE moduleConstParam, ONLY: PI IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(inout):: self - REAL(8):: r, xi(1:3) + CLASS(meshCell2DCylTria), INTENT(inout):: self + REAL(8):: r, Xi(1:3) REAL(8):: detJ 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 /) - detJ = self%detJac(xi)*PI !2PI*1/2 - CALL self%fPsi(xi, fPsi) + Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + detJ = self%detJac(Xi)*PI !2PI*1/2 + fPsi = self%fPsi(Xi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ @@ -728,37 +713,39 @@ MODULE moduleMesh2DCyl END SUBROUTINE areaTria !Shape functions for triangular element - PURE SUBROUTINE fPsiTria(xi, fPsi) + PURE FUNCTION fPsiTria(self, Xi) RESULT(fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: fPsi(1:self%nNodes) - fPsi(1) = 1.D0 - xi(1) - xi(2) - fPsi(2) = xi(1) - fPsi(3) = xi(2) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) - END SUBROUTINE fPsiTria + END FUNCTION fPsiTria - !Derivative element function at coordinates xi - PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiTria(self, Xi) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: dPsi(1:3,1:self%nNodes) - ALLOCATE(dPsi(1:2,1:3)) + dPsi = 0.D0 - dPsi(1,:) = dPsiTriaXi1(xi(2)) - dPsi(2,:) = dPsiTriaXi2(xi(1)) + dPsi(1,:) = dPsiTriaXi1(Xi(2)) + dPsi(2,:) = dPsiTriaXi2(Xi(1)) END FUNCTION dPsiTria - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + !Derivative element function (Xi1) + PURE FUNCTION dPsiTriaXi1(Xi2) RESULT(dPsiXi1) IMPLICIT NONE - REAL(8), INTENT(in):: xi2 + REAL(8), INTENT(in):: Xi2 REAL(8):: dPsiXi1(1:3) dPsiXi1(1) = -1.D0 @@ -767,11 +754,11 @@ MODULE moduleMesh2DCyl END FUNCTION dPsiTriaXi1 - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + !Derivative element function (Xi1) + PURE FUNCTION dPsiTriaXi2(Xi1) RESULT(dPsiXi2) IMPLICIT NONE - REAL(8), INTENT(in):: xi1 + REAL(8), INTENT(in):: Xi1 REAL(8):: dPsiXi2(1:3) dPsiXi2(1) = -1.D0 @@ -783,8 +770,8 @@ MODULE moduleMesh2DCyl PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) @@ -799,24 +786,23 @@ MODULE moduleMesh2DCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: r, xi(1:3) - REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) + REAL(8):: r, Xi(1:3) + REAL(8):: fPsi(1:3), dPsi(1:3,1:3) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l - ALLOCATE(localK(1:4, 1:4)) localK=0.D0 - xi=0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - CALL self%fPsi(xi, fPsi) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi,dPsi) + invJ = self%invJac(Xi,dPsi) + fPsi = self%fPsi(Xi) r = DOT_PRODUCT(fPsi,self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ @@ -830,23 +816,22 @@ MODULE moduleMesh2DCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: fPsi(1:3) - REAL(8):: r, xi(1:3) + REAL(8):: r, Xi(1:3) REAL(8):: detJ, f INTEGER:: l - ALLOCATE(localF(1:3)) localF = 0.D0 - xi = 0.D0 + Xi = 0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - detJ = self%detJac(xi) - CALL self%fPsi(xi, fPsi) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + detJ = self%detJac(Xi) + fPsi = self%fPsi(Xi) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wTria(l)*detJ @@ -856,112 +841,99 @@ MODULE moduleMesh2DCyl END FUNCTION elemFTria - PURE FUNCTION insideTria(xi) RESULT(ins) + PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = xi(1) >= 0.D0 .AND. & - xi(2) >= 0.D0 .AND. & - 1.D0 - xi(1) - xi(2) >= 0.D0 - - END FUNCTION insideTria - - !Gathers the electric field at position xi - PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1:2,1:3) - REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) REAL(8):: phi(1:3) - REAL(8):: EF(1:3) - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi /) + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi /) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = 0.D0 + array = -self%gatherDF(Xi, phi) END FUNCTION gatherEFTria - PURE FUNCTION gatherMFTria(self,xi) RESULT(MF) + PURE FUNCTION gatherMFTria(self, Xi) RESULT(array) IMPLICIT NONE + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:3,1:3) - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:3) - REAL(8):: MF_Nodes(1:3,1:3) - REAL(8):: MF(1:3) + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1) /) - MF_Nodes(1:3,1) = (/self%n1%emData%B(1), & - self%n2%emData%B(1), & - self%n3%emData%B(1) /) - MF_Nodes(1:3,2) = (/self%n1%emData%B(2), & - self%n2%emData%B(2), & - self%n3%emData%B(2) /) - MF_Nodes(1:3,3) = (/self%n1%emData%B(3), & - self%n2%emData%B(3), & - self%n3%emData%B(3) /) + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2) /) - CALL self%fPsi(xi, fPsi) - MF = MATMUL(fPsi, MF_Nodes) + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) END FUNCTION gatherMFTria + PURE FUNCTION insideTria(Xi) RESULT(ins) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + LOGICAL:: ins + + ins = Xi(1) >= 0.D0 .AND. & + Xi(2) >= 0.D0 .AND. & + 1.D0 - Xi(1) - Xi(2) >= 0.D0 + + END FUNCTION insideTria + !Gets node indexes from triangular element PURE FUNCTION getNodesTria(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell2DCylTria), INTENT(in):: self + INTEGER:: n(1:self%nNodes) - ALLOCATE(n(1:3)) n = (/self%n1%n, self%n2%n, self%n3%n /) END FUNCTION getNodesTria !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logTria(self,r) RESULT(xi) + PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self + CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: invJ(1:2,1:2), detJ REAL(8):: deltaR(1:2) - REAL(8):: dPsi(1:2,1:3) + REAL(8):: dPsi(1:3,1:3) !Direct method to convert coordinates - xi = 0.D0 !Irrelevant, required for input + Xi = 0.D0 !Irrelevant, required for input deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - detJ = self%detJac(xi, dPsi) - xi(1:2) = MATMUL(invJ,deltaR)/detJ + dPsi = self%dPsi(Xi) + invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, dPsi) + Xi(1:2) = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria - SUBROUTINE nextElementTria(self, xi, nextElement) + SUBROUTINE nextElementTria(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:3) + REAL(8):: XiArray(1:3) INTEGER:: nextInt - xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) - nextInt = MINLOC(xiArray,1) + XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray,1) NULLIFY(nextElement) SELECT CASE (nextInt) CASE (1) @@ -976,21 +948,21 @@ MODULE moduleMesh2DCyl !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS !Computes element Jacobian determinant - PURE FUNCTION detJ2DCyl(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCyl(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol2DCyl), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell2DCyl), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) REAL(8):: dJ + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dz(1:2), dr(1:2) IF(PRESENT(dPsi_in)) THEN dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF @@ -1000,27 +972,30 @@ MODULE moduleMesh2DCyl END FUNCTION detJ2DCyl !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCyl(self,xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCyl(self,Xi,dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol2DCyl), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell2DCyl), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: invJ(1:3,1:3) + REAL(8):: dPsi(1:3,1:self%nNodes) REAL(8):: dz(1:2), dr(1:2) - REAL(8):: invJ(1:2,1:2) IF(PRESENT(dPsi_in)) THEN dPsi=dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF + invJ = 0.D0 + CALL self%partialDer(dPsi, dz, dr) - invJ(1,:) = (/ dr(2), -dz(2) /) - invJ(2,:) = (/ -dr(1), dz(1) /) + + invJ(1,1:2) = (/ dr(2), -dz(2) /) + invJ(2,1:2) = (/ -dr(1), dz(1) /) END FUNCTION invJ2DCyl @@ -1030,11 +1005,11 @@ MODULE moduleMesh2DCyl CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et - DO e = 1, self%numVols + DO e = 1, self%numCells !Connect Vol-Vol - DO et = 1, self%numVols + DO et = 1, self%numCells IF (e /= et) THEN - CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj) END IF @@ -1044,7 +1019,7 @@ MODULE moduleMesh2DCyl TYPE IS(meshParticles) !Connect Vol-Edge DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj) END DO @@ -1058,31 +1033,31 @@ MODULE moduleMesh2DCyl SUBROUTINE connectVolVol(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVol2DCylQuad) + TYPE IS(meshCell2DCylQuad) !Element A is a quadrilateral SELECT TYPE(elemB) - TYPE IS(meshVol2DCylQuad) + TYPE IS(meshCell2DCylQuad) !Element B is a quadrilateral CALL connectQuadQuad(elemA, elemB) - TYPE IS(meshVol2DCylTria) + TYPE IS(meshCell2DCylTria) !Element B is a triangle CALL connectQuadTria(elemA, elemB) END SELECT - TYPE IS(meshVol2DCylTria) + TYPE IS(meshCell2DCylTria) !Element A is a Triangle SELECT TYPE(elemB) - TYPE IS(meshVol2DCylQuad) + TYPE IS(meshCell2DCylQuad) !Element B is a quadrilateral CALL connectQuadTria(elemB, elemA) - TYPE IS(meshVol2DCylTria) + TYPE IS(meshCell2DCylTria) !Element B is a triangle CALL connectTriaTria(elemA, elemB) @@ -1095,17 +1070,17 @@ MODULE moduleMesh2DCyl SUBROUTINE connectVolEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemA CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemB) CLASS IS(meshEdge2DCyl) SELECT TYPE(elemA) - TYPE IS(meshVol2DCylQuad) + TYPE IS(meshCell2DCylQuad) !Element A is a quadrilateral CALL connectQuadEdge(elemA, elemB) - TYPE IS(meshVol2DCylTria) + TYPE IS(meshCell2DCylTria) !Element A is a triangle CALL connectTriaEdge(elemA, elemB) @@ -1118,8 +1093,8 @@ MODULE moduleMesh2DCyl SUBROUTINE connectQuadQuad(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1) .AND. & @@ -1162,8 +1137,8 @@ MODULE moduleMesh2DCyl SUBROUTINE connectQuadTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -1254,8 +1229,8 @@ MODULE moduleMesh2DCyl SUBROUTINE connectTriaTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB + CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -1327,7 +1302,7 @@ MODULE moduleMesh2DCyl SUBROUTINE connectQuadEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB !Check direction 1 @@ -1411,7 +1386,7 @@ MODULE moduleMesh2DCyl SUBROUTINE connectTriaEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemA CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB !Check direction 1 diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 4b29c16..9bd6468 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -31,26 +31,19 @@ MODULE moduleMesh3DCart END TYPE meshEdge3DCartTria - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol3DCart + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell3DCart CONTAINS PROCEDURE, PASS:: detJac => detJ3DCart PROCEDURE, PASS:: invJac => invJ3DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - END TYPE meshVol3DCart + END TYPE meshCell3DCart ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi) - REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - PURE SUBROUTINE partialDer_interface(self, dPsi, dx, dy, dz) - IMPORT meshVol3DCart - CLASS(meshVol3DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) + IMPORT meshCell3DCart + CLASS(meshCell3DCart), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz END SUBROUTINE partialDer_interface @@ -58,7 +51,7 @@ MODULE moduleMesh3DCart END INTERFACE !Tetrahedron volume element - TYPE, PUBLIC, EXTENDS(meshVol3DCart):: meshVol3DCartTetra + TYPE, PUBLIC, EXTENDS(meshCell3DCart):: meshCell3DCartTetra !Element Coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -66,24 +59,24 @@ MODULE moduleMesh3DCart !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() CONTAINS - PROCEDURE, PASS:: init => initVolTetra3DCart - PROCEDURE, PASS:: randPos => randPosVolTetra - PROCEDURE, PASS:: calcVol => volumeTetra - PROCEDURE, NOPASS:: fPsi => fPsiTetra - PROCEDURE, NOPASS:: dPsi => dPsiTetra - PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXi2 + PROCEDURE, PASS:: init => initCellTetra3DCart + PROCEDURE, PASS:: randPos => randPosCellTetra + PROCEDURE, PASS:: calcCell => volumeTetra + PROCEDURE, PASS:: fPsi => fPsiTetra + PROCEDURE, PASS:: dPsi => dPsiTetra + PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiTetraXi1 + PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiTetraXi2 PROCEDURE, PASS:: partialDer => partialDerTetra PROCEDURE, PASS:: elemK => elemKTetra PROCEDURE, PASS:: elemF => elemFTetra + PROCEDURE, PASS:: gatherElectricField => gatherEFTetra + PROCEDURE, PASS:: gatherMagneticField => gatherMFTetra PROCEDURE, NOPASS:: inside => insideTetra - PROCEDURE, PASS:: gatherEF => gatherEFTetra - PROCEDURE, PASS:: gatherMF => gatherMFTetra PROCEDURE, PASS:: getNodes => getNodesTetra PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: nextElement => nextElementTetra - END TYPE meshVol3DCartTetra + END TYPE meshCell3DCartTetra CONTAINS !NODE FUNCTIONS @@ -245,19 +238,20 @@ MODULE moduleMesh3DCart !VOLUME FUNCTIONS !TETRA FUNCTIONS !Inits tetrahedron element - SUBROUTINE initVolTetra3DCart(self, n, p, nodes) + SUBROUTINE initCellTetra3DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(out):: self + CLASS(meshCell3DCartTetra), INTENT(out):: self INTEGER, INTENT(in):: n 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) !Volume of each node + REAL(8):: volNodes(1:4) !Cellume of each node self%n = n + self%nNodes = SIZE(p) self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -272,11 +266,11 @@ MODULE moduleMesh3DCart self%z = (/r1(3), r2(3), r3(3), r4(3)/) !Computes the element volume - CALL self%calcVol() + CALL self%calcCell() !Assign proportional volume to each node Xi = (/0.25D0, 0.25D0, 0.25D0/) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) volNodes = fPsi*self%volume self%n1%v = self%n1%v + volNodes(1) self%n2%v = self%n2%v + volNodes(2) @@ -288,14 +282,14 @@ MODULE moduleMesh3DCart ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolTetra3DCart + END SUBROUTINE initCellTetra3DCart !Random position in volume tetrahedron - FUNCTION randPosVolTetra(self) RESULT(r) + FUNCTION randPosCellTetra(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self + CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) @@ -304,19 +298,19 @@ MODULE moduleMesh3DCart Xi(2) = random( 0.D0, 1.D0 - Xi(1)) Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2)) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) r = (/ DOT_PRODUCT(fPsi, self%x), & DOT_PRODUCT(fPsi, self%y), & DOT_PRODUCT(fPsi, self%z) /) - END FUNCTION randPosVolTetra + END FUNCTION randPosCellTetra !Computes the element volume PURE SUBROUTINE volumeTetra(self) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(inout):: self + CLASS(meshCell3DCartTetra), INTENT(inout):: self REAL(8):: Xi(1:3) self%volume = 0.D0 @@ -326,27 +320,29 @@ MODULE moduleMesh3DCart END SUBROUTINE volumeTetra !Computes element functions in point Xi - PURE SUBROUTINE fPsiTetra(Xi, fPsi) + PURE FUNCTION fPsiTetra(self, Xi) RESULT(fPsi) IMPLICIT NONE + CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + REAL(8):: fPsi(1:self%nNodes) fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3) fPsi(2) = Xi(1) fPsi(3) = Xi(2) fPsi(4) = Xi(3) - END SUBROUTINE fPsiTetra + END FUNCTION fPsiTetra !Derivative element function at coordinates Xi - PURE FUNCTION dPsiTetra(Xi) RESULT(dPsi) + PURE FUNCTION dPsiTetra(self, Xi) RESULT(dPsi) IMPLICIT NONE + CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3, 1:self%nNodes) - ALLOCATE(dPsi(1:3,1:4)) + dPsi = 0.D0 dPsi(1,:) = dPsiTetraXi1(Xi(2), Xi(3)) dPsi(2,:) = dPsiTetraXi2(Xi(1), Xi(3)) @@ -397,8 +393,8 @@ MODULE moduleMesh3DCart PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:, 1:) + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:3, 1:self%nNodes) REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) @@ -418,13 +414,12 @@ MODULE moduleMesh3DCart PURE FUNCTION elemKTetra(self) RESULT(localK) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) REAL(8):: invJ(1:3,1:3), detJ - ALLOCATE(localK(1:4,1:4)) localK = 0.D0 Xi = 0.D0 !TODO: One point Gauss integral. Upgrade when possible @@ -432,7 +427,7 @@ MODULE moduleMesh3DCart dPsi = self%dPsi(Xi) detJ = self%detJac(Xi, dPsi) invJ = self%invJac(Xi, dPsi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra @@ -440,26 +435,66 @@ MODULE moduleMesh3DCart PURE FUNCTION elemFTetra(self, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) REAL(8):: Xi(1:3) REAL(8):: detJ, f - ALLOCATE(localF(1:4)) - localF = 0.D0 Xi = 0.D0 Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) dPsi = self%dPsi(Xi) detJ = self%detJac(Xi, dPsi) - CALL self%fPsi(Xi, fPsi) + fPsi = self%fPsi(Xi) f = DOT_PRODUCT(fPsi, source) localF = f*fPsi*1.D0*detJ END FUNCTION elemFTetra + PURE FUNCTION gatherEFTetra(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: phi(1:4) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) + + array = -self%gatherDF(Xi, phi) + + END FUNCTION gatherEFTetra + + PURE FUNCTION gatherMFTetra(self, Xi) RESULT(array) + IMPLICIT NONE + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + REAL(8):: B(1:4,1:3) + + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) + + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) + + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFTetra + PURE FUNCTION insideTetra(Xi) RESULT(ins) IMPLICIT NONE @@ -473,66 +508,12 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra - PURE FUNCTION gatherEFTetra(self, Xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: dPsi(1:3, 1:4) - REAL(8):: dPsiR(1:3, 1:4) - REAL(8):: invJ(1:3, 1:3), detJ - REAL(8):: phi(1:4) - REAL(8):: EF(1:3) - - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi, & - self%n4%emData%phi /) - - dPsi = self%dPsi(Xi) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = -DOT_PRODUCT(dPsiR(3,:), phi) - - END FUNCTION gatherEFTetra - - PURE FUNCTION gatherMFTetra(self, Xi) RESULT(MF) - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - REAL(8):: fPsi(1:4) - REAL(8):: MF_Nodes(1:4,1:3) - REAL(8):: MF(1:3) - - MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & - self%n2%emData%B(1), & - self%n3%emData%B(1), & - self%n4%emData%B(1) /) - MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & - self%n2%emData%B(2), & - self%n3%emData%B(2), & - self%n4%emData%B(2) /) - MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & - self%n2%emData%B(3), & - self%n3%emData%B(3), & - self%n4%emData%B(3) /) - - CALL self%fPsi(Xi, fPsi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFTetra - PURE FUNCTION getNodesTetra(self) RESULT(n) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + CLASS(meshCell3DCartTetra), INTENT(in):: self + INTEGER:: n(1:self%nnodes) - ALLOCATE(n(1:4)) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) END FUNCTION getNodesTetra @@ -540,7 +521,7 @@ MODULE moduleMesh3DCart PURE FUNCTION phy2logTetra(self,r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self + CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) REAL(8):: invJ(1:3, 1:3), detJ @@ -559,7 +540,7 @@ MODULE moduleMesh3DCart SUBROUTINE nextElementTetra(self, Xi, nextElement) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self + CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement REAL(8):: XiArray(1:4) @@ -587,11 +568,11 @@ MODULE moduleMesh3DCart PURE FUNCTION detJ3DCart(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol3DCart), INTENT(in)::self + CLASS(meshCell3DCart), INTENT(in)::self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:, 1:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:self%nNodes) REAL(8):: dJ - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dPsi(1:3, 1:self%nNodes) REAL(8):: dx(1:3), dy(1:3), dz(1:3) IF (PRESENT(dPsi_in)) THEN @@ -612,10 +593,10 @@ MODULE moduleMesh3DCart PURE FUNCTION invJ3DCart(self,Xi,dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol3DCart), INTENT(in):: self + CLASS(meshCell3DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:self%nNodes) + REAL(8):: dPsi(1:3, 1:self%nNodes) REAL(8), DIMENSION(1:3):: dx, dy, dz REAL(8):: invJ(1:3,1:3) @@ -645,17 +626,17 @@ MODULE moduleMesh3DCart END FUNCTION invJ3DCart !Selects type of elements to build connection - SUBROUTINE connectVolVol(elemA, elemB) + SUBROUTINE connectCellCell(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) + TYPE IS(meshCell3DCartTetra) !Element A is a tetrahedron SELECT TYPE(elemB) - TYPE IS(meshVol3DCartTetra) + TYPE IS(meshCell3DCartTetra) !Element B is a tetrahedron CALL connectTetraTetra(elemA, elemB) @@ -663,18 +644,18 @@ MODULE moduleMesh3DCart END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell - SUBROUTINE connectVolEdge(elemA, elemB) + SUBROUTINE connectCellEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemA CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemB) CLASS IS(meshEdge3DCartTria) SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) + TYPE IS(meshCell3DCartTetra) !Element A is a tetrahedron CALL connectTetraEdge(elemA, elemB) @@ -682,7 +663,7 @@ MODULE moduleMesh3DCart END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge SUBROUTINE connectMesh3DCart(self) IMPLICIT NONE @@ -690,11 +671,11 @@ MODULE moduleMesh3DCart CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et - DO e = 1, self%numVols - !Connect Vol-Vol - DO et = 1, self%numVols + DO e = 1, self%numCells + !Connect Cell-Cell + DO et = 1, self%numCells IF (e /= et) THEN - CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj) + CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj) END IF @@ -702,9 +683,9 @@ MODULE moduleMesh3DCart SELECT TYPE(self) TYPE IS(meshParticles) - !Connect Vol-Edge + !Connect Cell-Edge DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj) END DO @@ -740,8 +721,8 @@ MODULE moduleMesh3DCart SUBROUTINE connectTetraTetra(elemA, elemB) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB + CLASS(meshCell3DCartTetra), INTENT(inout), TARGET:: elemA + CLASS(meshCell3DCartTetra), INTENT(inout), TARGET:: elemB !Check surface 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -869,11 +850,11 @@ MODULE moduleMesh3DCart USE moduleMath IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA + CLASS(meshCell3DCartTetra), 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) + REAL(8):: normCell(1:3) nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /) @@ -888,10 +869,10 @@ MODULE moduleMesh3DCart 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) + normCell = crossProduct(vec1, vec2) + normCell = normalize(normCell) - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN elemA%e1 => elemB elemB%e1 => elemA @@ -921,10 +902,10 @@ MODULE moduleMesh3DCart 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) + normCell = crossProduct(vec1, vec2) + normCell = normalize(normCell) - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN elemA%e2 => elemB elemB%e1 => elemA @@ -954,10 +935,10 @@ MODULE moduleMesh3DCart 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) + normCell = crossProduct(vec1, vec2) + normCell = normalize(normCell) - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN elemA%e3 => elemB elemB%e1 => elemA @@ -987,10 +968,10 @@ MODULE moduleMesh3DCart 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) + normCell = crossProduct(vec1, vec2) + normCell = normalize(normCell) - IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN elemA%e4 => elemB elemB%e1 => elemA diff --git a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 index 5ac0682..37dbf82 100644 --- a/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshInput0D.f90 @@ -41,8 +41,8 @@ MODULE moduleMeshInput0D self%numNodes = 1 ALLOCATE(self%nodes(1:1)) !Allocates one volume - self%numVols = 1 - ALLOCATE(self%vols(1:1)) + self%numCells = 1 + ALLOCATE(self%cells(1:1)) !Allocates matrix K SELECT TYPE(self) TYPE IS(meshParticles) @@ -59,8 +59,8 @@ MODULE moduleMeshInput0D CALL self%nodes(1)%obj%init(1, r) !Creates the volume - ALLOCATE(meshVol0D:: self%vols(1)%obj) - CALL self%vols(1)%obj%init(1, (/ 1/), self%nodes) + ALLOCATE(meshCell0D:: self%cells(1)%obj) + CALL self%cells(1)%obj%init(1, (/ 1/), self%nodes) END SUBROUTINE read0D diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index 15045b1..dfe9605 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -57,7 +57,7 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%vols(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) + WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 7832843..aae2216 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -136,7 +136,7 @@ MODULE moduleMeshInputGmsh2 !Substract the number of edges to the total number of elements !to obtain the number of volume elements - self%numVols = TotalnumElem - self%numEdges + self%numCells = TotalnumElem - self%numEdges ALLOCATE(self%edges(1:self%numEdges)) numEdges = self%numEdges @@ -146,13 +146,13 @@ MODULE moduleMeshInputGmsh2 END DO TYPE IS(meshCollisions) - self%numVols = TotalnumElem + self%numCells = TotalnumElem numEdges = 0 END SELECT !Allocates arrays - ALLOCATE(self%vols(1:self%numVols)) + ALLOCATE(self%cells(1:self%numCells)) SELECT TYPE(self) TYPE IS(meshParticles) @@ -232,7 +232,7 @@ MODULE moduleMeshInputGmsh2 END SELECT !Read and initialize volumes - DO e = 1, self%numVols + DO e = 1, self%numCells !Reads the volume according to the geometry SELECT CASE(self%dimen) CASE(3) @@ -244,7 +244,7 @@ MODULE moduleMeshInputGmsh2 !Tetrahedron element ALLOCATE(p(1:4)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol3DCartTetra:: self%vols(e)%obj) + ALLOCATE(meshCell3DCartTetra:: self%cells(e)%obj) END SELECT @@ -259,13 +259,13 @@ MODULE moduleMeshInputGmsh2 !Triangular element ALLOCATE(p(1:3)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj) + ALLOCATE(meshCell2DCylTria:: self%cells(e)%obj) CASE (3) !Quadrilateral element ALLOCATE(p(1:4)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj) + ALLOCATE(meshCell2DCylQuad:: self%cells(e)%obj) END SELECT @@ -278,13 +278,13 @@ MODULE moduleMeshInputGmsh2 !Triangular element ALLOCATE(p(1:3)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj) + ALLOCATE(meshCell2DCartTria:: self%cells(e)%obj) CASE (3) !Quadrilateral element ALLOCATE(p(1:4)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj) + ALLOCATE(meshCell2DCartQuad:: self%cells(e)%obj) END SELECT @@ -296,19 +296,19 @@ MODULE moduleMeshInputGmsh2 ALLOCATE(p(1:2)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) - ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj) + ALLOCATE(meshCell1DRadSegm:: self%cells(e)%obj) CASE("Cart") ALLOCATE(p(1:2)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) - ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj) + ALLOCATE(meshCell1DCartSegm:: self%cells(e)%obj) END SELECT END SELECT - CALL self%vols(e)%obj%init(n - numEdges, p, self%nodes) + CALL self%cells(e)%obj%init(n - numEdges, p, self%nodes) DEALLOCATE(p) END DO diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 7eade3e..53485f4 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -181,9 +181,9 @@ MODULE moduleMeshOutputGmsh2 DO c = 1, interactionMatrix(k)%amount WRITE(cString, "(I2)") c title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString - CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numVols) - DO n=1, self%numVols - WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c) + CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells) + DO n=1, self%numCells + WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c) END DO CALL writeGmsh2FooterElementData(60) @@ -211,9 +211,9 @@ MODULE moduleMeshOutputGmsh2 REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName CHARACTER (LEN=iterationDigits):: tstring - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) - xi = (/ 0.D0, 0.D0, 0.D0 /) + Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN time = DBLE(t)*tauMin*ti_ref @@ -231,9 +231,9 @@ MODULE moduleMeshOutputGmsh2 END DO CALL writeGmsh2FooterNodeData(20) - CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numVols) - DO e=1, self%numVols - WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells) + DO e=1, self%numCells + WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref END DO CALL writeGmsh2FooterElementData(20) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 11e8bc7..d0a03f4 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -66,10 +66,10 @@ MODULE moduleMesh !Parent of Edge element TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge - !Connectivity to vols - CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL() - !Connectivity to vols in meshColl - CLASS(meshVol), POINTER:: eColl => NULL() + !Connectivity to cells + CLASS(meshCell), POINTER:: e1 => NULL(), e2 => NULL() + !Connectivity to cells in meshColl + CLASS(meshCell), POINTER:: eColl => NULL() !Normal vector REAL(8):: normal(1:3) !Weight for random injection of particles @@ -146,8 +146,10 @@ MODULE moduleMesh END TYPE meshEdgeCont - !Parent of Volume element - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol + !Parent of cell element + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshCell + !Number of nodes in the cell + INTEGER:: nNodes !Maximum collision rate REAL(8), ALLOCATABLE:: sigmaVrelMax(:) !Arrays for counting number of collisions @@ -161,114 +163,152 @@ MODULE moduleMesh !Total weight of particles inside cell REAL(8), ALLOCATABLE:: totalWeight(:) CONTAINS - PROCEDURE(initVol_interface), DEFERRED, PASS:: init - PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes - PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos - PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi - PROCEDURE, PASS:: scatter - PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF - PROCEDURE(gatherMF_interface), DEFERRED, PASS:: gatherMF - PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK - PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF + !Init the cell + PROCEDURE(initCell_interface), DEFERRED, PASS:: init + !Get the index of the nodes + PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes + !Calculate random position on the cell + PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos + !Obtain functions and values of cell natural functions + PROCEDURE(fPsi_interface), DEFERRED, PASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, PASS:: dPsi + PROCEDURE(detJac_interface), DEFERRED, PASS:: detJac + PROCEDURE(invJac_interface), DEFERRED, PASS:: invJac + !Scatter properties of particles on cell nodes + PROCEDURE, PASS:: scatter + !Gather value and spatial derivative on the nodes at position Xi + PROCEDURE, PASS, PRIVATE:: gatherF_scalar + PROCEDURE, PASS, PRIVATE:: gatherF_array + PROCEDURE, PASS, PRIVATE:: gatherDF_scalar + GENERIC:: gatherF => gatherF_scalar, gatherF_array + GENERIC:: gatherDF => gatherDF_scalar + !Procedures to get specific values in the node + PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherElectricField + PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherMagneticField + !Compute K and F to solve PDE on the mesh + PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK + PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF + !Subroutines to find in which cell a particle is located PROCEDURE, PASS:: findCell - PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement + !Convert physical coordinates (r) into logical coordinates (Xi) + PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log - END TYPE meshVol + END TYPE meshCell ABSTRACT INTERFACE - SUBROUTINE initVol_interface(self, n, p, nodes) - IMPORT:: meshVol + SUBROUTINE initCell_interface(self, n, p, nodes) + IMPORT:: meshCell IMPORT meshNodeCont - CLASS(meshVol), INTENT(out):: self + CLASS(meshCell), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) - END SUBROUTINE initVol_interface - - PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: EF(1:3) - - END FUNCTION gatherEF_interface - - PURE FUNCTION gatherMF_interface(self, xi) RESULT(MF) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: MF(1:3) - - END FUNCTION gatherMF_interface + END SUBROUTINE initCell_interface PURE FUNCTION getNodesVol_interface(self) RESULT(n) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + INTEGER:: n(1:self%nNodes) END FUNCTION getNodesVol_interface - PURE SUBROUTINE fPsi_interface(xi, fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(out):: fPsi(:) + PURE FUNCTION fPsi_interface(self, Xi) RESULT(fPsi) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: fPsi(1:self%nNodes) - END SUBROUTINE fPsi_interface + END FUNCTION fPsi_interface + + PURE FUNCTION dPsi_interface(self, Xi) RESULT(dPsi) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:self%nNodes) + + END FUNCTION dPsi_interface + + PURE FUNCTION detJac_interface(self, Xi, dPsi_in) RESULT(dJ) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: dJ + + END FUNCTION detJac_interface + + PURE FUNCTION invJac_interface(self, Xi, dPsi_in) RESULT(invJ) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + REAL(8):: invJ(1:3,1:3) + + END FUNCTION invJac_interface + + PURE FUNCTION gatherArray_interface(self, Xi) RESULT(array) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) + + END FUNCTION gatherArray_interface PURE FUNCTION elemK_interface(self) RESULT(localK) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8):: localK(1:self%nNodes,1:self%nNodes) END FUNCTION elemK_interface PURE FUNCTION elemF_interface(self, source) RESULT(localF) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: source(1:self%nNodes) + REAL(8):: localF(1:self%nNodes) END FUNCTION elemF_interface - SUBROUTINE nextElement_interface(self, xi, nextElement) - IMPORT:: meshVol, meshElement - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + SUBROUTINE nextElement_interface(self, Xi, nextElement) + IMPORT:: meshCell, meshElement + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement END SUBROUTINE nextElement_interface - PURE FUNCTION phy2log_interface(self,r) RESULT(xN) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self + PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) + REAL(8):: Xi(1:3) END FUNCTION phy2log_interface - PURE FUNCTION inside_interface(xi) RESULT(ins) - IMPORT:: meshVol - REAL(8), INTENT(in):: xi(1:3) + PURE FUNCTION inside_interface(Xi) RESULT(ins) + IMPORT:: meshCell + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins END FUNCTION inside_interface FUNCTION randPosVol_interface(self) RESULT(r) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self REAL(8):: r(1:3) END FUNCTION randPosVol_interface END INTERFACE - !Containers for volumes in the mesh - TYPE:: meshVolCont - CLASS(meshVol), ALLOCATABLE:: obj + !Containers for cells in the mesh + TYPE:: meshCellCont + CLASS(meshCell), ALLOCATABLE:: obj - END TYPE meshVolCont + END TYPE meshCellCont !Generic mesh type TYPE, ABSTRACT:: meshGeneric @@ -277,11 +317,11 @@ MODULE moduleMesh !Geometry of the mesh CHARACTER(:), ALLOCATABLE:: geometry !Number of elements - INTEGER:: numNodes, numVols + INTEGER:: numNodes, numCells !Array of nodes TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) - !Array of volume elements - TYPE(meshVolCont), ALLOCATABLE:: vols(:) + !Array of cell elements + TYPE(meshCellCont), ALLOCATABLE:: cells(:) PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() @@ -310,7 +350,7 @@ MODULE moduleMesh END SUBROUTINE readInitial_interface - !Connects volume and edges to the mesh + !Connects cell and edges to the mesh SUBROUTINE connectMesh_interface(self) IMPORT meshGeneric @@ -318,7 +358,7 @@ MODULE moduleMesh END SUBROUTINE connectMesh_interface - !Prints number of collisions in each volume + !Prints number of collisions in each cell SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric @@ -416,7 +456,7 @@ MODULE moduleMesh !Pointer to mesh used for MC collisions CLASS(meshGeneric), POINTER:: meshForMCC => NULL() - !Procedure to find a volume for a particle in meshColl + !Procedure to find a cell for a particle in meshColl PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL() ABSTRACT INTERFACE @@ -445,9 +485,9 @@ MODULE moduleMesh REAL(8), ALLOCATABLE:: localK(:,:) INTEGER:: nNodes, i, j - DO e = 1, self%numVols - n = self%vols(e)%obj%getNodes() - localK = self%vols(e)%obj%elemK() + DO e = 1, self%numCells + n = self%cells(e)%obj%getNodes() + localK = self%cells(e)%obj%elemK() nNodes = SIZE(n) DO i = 1, nNodes @@ -480,33 +520,84 @@ MODULE moduleMesh END SUBROUTINE resetOutput - !Scatters particle properties into vol nodes + !Gather the value of valNodes (scalar) at position Xi + PURE FUNCTION gatherF_scalar(self, Xi, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in):: valNodes(1:self%nNodes) + REAL(8):: f + REAL(8):: fPsi(1:self%nNodes) + + fPsi = self%fPsi(Xi) + f = DOT_PRODUCT(fPsi, valNodes) + + END FUNCTION gatherF_scalar + + !Gather the value of valNodes (array) at position Xi + PURE FUNCTION gatherF_array(self, Xi, n, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: valNodes(1:self%nNodes, 1:n) + REAL(8):: f(1:n) + REAL(8):: fPsi(1:self%nNodes) + + fPsi = self%fPsi(Xi) + f = MATMUL(fPsi, valNodes) + + END FUNCTION gatherF_array + + !Gather the spatial derivative of valNodes (scalar) at position Xi + PURE FUNCTION gatherDF_scalar(self, Xi, valNodes) RESULT(df) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(in):: valNodes(1:self%nNodes) + REAL(8):: df(1:3) + REAL(8):: dPsi(1:3, 1:self%nNodes) + REAL(8):: dPsiR(1:3, 1:self%nNodes) + REAL(8):: invJ(1:3, 1:3), detJ + + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) + dPsiR = MATMUL(invJ, dPsi)/detJ + df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), & + DOT_PRODUCT(dPsiR(2,:), valNodes), & + DOT_PRODUCT(dPsiR(3,:), valNodes) /) + + END FUNCTION gatherDF_scalar + + !Scatters particle properties into cell nodes SUBROUTINE scatter(self, part) USE moduleMath USE moduleSpecies USE OMP_LIB IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: self + CLASS(meshCell), INTENT(inout):: self CLASS(particle), INTENT(in):: part - REAL(8), ALLOCATABLE:: fPsi(:) - INTEGER, ALLOCATABLE:: volNodes(:) + REAL(8):: fPsi(1:self%nNodes) + INTEGER:: cellNodes(1:self%nNodes) REAL(8):: tensorS(1:3, 1:3) INTEGER:: sp - INTEGER:: i, nNodes + INTEGER:: i CLASS(meshNode), POINTER:: node - volNodes = self%getNodes() - nNodes = SIZE(volNodes) - ALLOCATE(fPsi(1:nNodes)) - CALL self%fPsi(part%xi, fPsi) + cellNodes = self%getNodes() + fPsi = self%fPsi(part%Xi) tensorS = outerProduct(part%v, part%v) sp = part%species%n - DO i = 1, nNodes - node => mesh%nodes(volNodes(i))%obj + DO i = 1, self%nNodes + node => mesh%nodes(cellNodes(i))%obj CALL OMP_SET_LOCK(node%lock) node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:) @@ -524,18 +615,18 @@ MODULE moduleMesh USE OMP_LIB IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: self + CLASS(meshCell), INTENT(inout):: self CLASS(particle), INTENT(inout), TARGET:: part - CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell - REAL(8):: xi(1:3) + CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell + REAL(8):: Xi(1:3) CLASS(meshElement), POINTER:: nextElement INTEGER:: sp - xi = self%phy2log(part%r) + Xi = self%phy2log(part%r) !Checks if particle is inside 'self' cell - IF (self%inside(xi)) THEN + IF (self%inside(Xi)) THEN part%vol = self%n - part%xi = xi + part%Xi = Xi part%n_in = .TRUE. !Assign particle to listPart_in CALL OMP_SET_LOCK(self%lock) @@ -546,10 +637,10 @@ MODULE moduleMesh ELSE !If not, searches for a neighbour and repeats the process. - CALL self%nextElement(xi, nextElement) + CALL self%nextElement(Xi, nextElement) !Defines the next step SELECT TYPE(nextElement) - CLASS IS(meshVol) + CLASS IS(meshCell) !Particle moved to new cell, repeat find procedure CALL nextElement%findCell(part, self) @@ -598,31 +689,31 @@ MODULE moduleMesh TYPE(particle), INTENT(inout):: part LOGICAL:: found - CLASS(meshVol), POINTER:: vol - REAL(8), DIMENSION(1:3):: xii + CLASS(meshCell), POINTER:: cell + REAL(8), DIMENSION(1:3):: Xi CLASS(meshElement), POINTER:: nextElement INTEGER:: sp found = .FALSE. - vol => meshColl%vols(part%volColl)%obj + cell => meshColl%cells(part%volColl)%obj DO WHILE(.NOT. found) - xii = vol%phy2log(part%r) - IF (vol%inside(xii)) THEN - part%volColl = vol%n - CALL OMP_SET_LOCK(vol%lock) + Xi = cell%phy2log(part%r) + IF (cell%inside(Xi)) THEN + part%volColl = cell%n + CALL OMP_SET_LOCK(cell%lock) sp = part%species%n - CALL vol%listPart_in(sp)%add(part) - vol%totalWeight(sp) = vol%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(vol%lock) + CALL cell%listPart_in(sp)%add(part) + cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(cell%lock) found = .TRUE. ELSE - CALL vol%nextElement(xii, nextElement) + CALL cell%nextElement(Xi, nextElement) SELECT TYPE(nextElement) - CLASS IS(meshVol) + CLASS IS(meshCell) !Try next element - vol => nextElement + cell => nextElement CLASS DEFAULT !Should never happend, but just in case, stops loops @@ -647,15 +738,15 @@ MODULE moduleMesh REAL(8), DIMENSION(1:3), INTENT(in):: r INTEGER:: nVol INTEGER:: e - REAL(8), DIMENSION(1:3):: xii + REAL(8), DIMENSION(1:3):: Xi !Inits RESULT nVol = 0 - DO e = 1, self%numVols - xii = self%vols(e)%obj%phy2log(r) - IF(self%vols(e)%obj%inside(xii)) THEN - nVol = self%vols(e)%obj%n + DO e = 1, self%numCells + Xi = self%cells(e)%obj%phy2log(r) + IF(self%cells(e)%obj%inside(Xi)) THEN + nVol = self%cells(e)%obj%n EXIT END IF @@ -678,7 +769,7 @@ MODULE moduleMesh CLASS(meshGeneric), INTENT(inout), TARGET:: self INTEGER, INTENT(in):: t INTEGER:: e - CLASS(meshVol), POINTER:: vol + CLASS(meshCell), POINTER:: cell INTEGER:: k, i, j INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision @@ -689,21 +780,22 @@ MODULE moduleMesh REAL(8):: vRel, rMass, eRel REAL(8):: sigmaVrelTotal REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) - REAL(8):: rnd !Random number for collision + REAL(8):: rnd_real !Random number for collision + INTEGER:: rnd_int !Random number for collision IF (MOD(t, everyColl) == 0) THEN !Collisions need to be performed in this iteration !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) - DO e=1, self%numVols + DO e=1, self%numCells - vol => self%vols(e)%obj + cell => self%cells(e)%obj !TODO: Simplify this, to many sublevels !Iterate over the number of pairs DO k = 1, nCollPairs !Reset tally of collisions IF (collOutput) THEN - vol%tallyColl(k)%tally = 0 + cell%tallyColl(k)%tally = 0 END IF @@ -713,8 +805,8 @@ MODULE moduleMesh j = interactionMatrix(k)%sp_j%n !Number of particles per species in the collision pair - nPart_i = vol%listPart_in(i)%amount - nPart_j = vol%listPart_in(j)%amount + nPart_i = cell%listPart_in(i)%amount + nPart_j = cell%listPart_in(j)%amount IF (nPart_i > 0 .AND. nPart_j > 0) THEN !Total number of particles for the collision pair @@ -724,15 +816,15 @@ MODULE moduleMesh nColl = 0 !Probability of collision for pair i-j - pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax(k)*tauColl/vol%volume + pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume !Number of collisions in the cell nColl = NINT(REAL(nPart)*pMax*0.5D0) !Converts the list of particles to an array for easy access IF (nColl > 0) THEN - partTemp_i = vol%listPart_in(i)%convert2Array() - partTemp_j = vol%listPart_in(j)%convert2Array() + partTemp_i = cell%listPart_in(i)%convert2Array() + partTemp_j = cell%listPart_in(j)%convert2Array() END IF @@ -740,10 +832,10 @@ MODULE moduleMesh !Select random particles part_i => NULL() part_j => NULL() - rnd = random(1, nPart_i) - part_i => partTemp_i(rnd)%part - rnd = random(1, nPart_j) - part_j => partTemp_j(rnd)%part + rnd_int = random(1, nPart_i) + part_i => partTemp_i(rnd_int)%part + rnd_int = random(1, nPart_j) + part_j => partTemp_j(rnd_int)%part !If they are the same particle, skip !TODO: Maybe try to improve this IF (ASSOCIATED(part_i, part_j)) THEN @@ -767,32 +859,32 @@ MODULE moduleMesh CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) !Update maximum sigma*v_rel - IF (sigmaVrelTotal > vol%sigmaVrelMax(k)) THEN - vol%sigmaVrelMax(k) = sigmaVrelTotal + IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN + cell%sigmaVrelMax(k) = sigmaVrelTotal END IF ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) probabilityColl = 0.0 DO c = 1, interactionMatrix(k)%amount - probabilityColl(c) = sigmaVrel(c)/vol%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1)) + probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1)) END DO !Selects random number between 0 and 1 - rnd = random() + rnd_real = random() !If the random number is below the total probability of collision, collide particles - IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax(k)) THEN + IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN !Loop over collisions DO c = 1, interactionMatrix(k)%amount - IF (rnd <= probabilityColl(c)) THEN + IF (rnd_real <= probabilityColl(c)) THEN CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) !If collisions are gonna be output, count the collision IF (collOutput) THEN - vol%tallyColl(k)%tally(c) = vol%tallyColl(k)%tally(c) + 1 + cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1 END IF diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 2e15c5c..cc5d78c 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -1,4 +1,4 @@ -!moduleMeshBoundary: Boundary functions +!moduleMeshBoundary: Boundary functions for the mesh edges MODULE moduleMeshBoundary USE moduleMesh @@ -159,7 +159,7 @@ MODULE moduleMeshBoundary newElectron%vol = part%vol newIon%vol = part%vol - newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) + newElectron%xi = mesh%cells(part%vol)%obj%phy2log(newElectron%r) newIon%xi = newElectron%xi newElectron%weight = part%weight diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index ba1edfd..224bfbb 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -439,7 +439,6 @@ MODULE moduleCollisions REAL(8), INTENT(in):: vRel TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), POINTER:: electron => NULL(), ion => NULL() - REAL(8):: sigmaVrel REAL(8), DIMENSION(1:3):: vp_i TYPE(particle), POINTER:: remainingIon => NULL() diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 6b8c7a0..050ad52 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -132,7 +132,7 @@ MODULE moduleInject IF (doubleMesh) THEN nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) IF (nVolColl > 0) THEN - mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj + mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj ELSE CALL criticalError("No connection between edge and meshColl", "initInject") @@ -305,7 +305,7 @@ MODULE moduleInject self%v(3)%obj%randomVel() /) !Obtain natural coordinates of particle in cell - partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r) + partInj(n)%Xi = mesh%cells(partInj(n)%vol)%obj%phy2log(partInj(n)%r) !Push new particle with the minimum time step CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) !Assign cell to new particle diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index ae73ebc..55eb618 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -46,40 +46,6 @@ MODULE moduleEM END SUBROUTINE - PURE FUNCTION gatherElecField(part) RESULT(elField) - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - TYPE(particle), INTENT(in):: part - REAl(8):: xi(1:3) !Logical coordinates of particle in element - REAL(8):: elField(1:3) - - elField = 0.D0 - - xi = part%xi - - elField = mesh%vols(part%vol)%obj%gatherEF(xi) - - END FUNCTION gatherElecField - - PURE FUNCTION gatherMagnField(part) RESULT(BField) - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - TYPE(particle), INTENT(in):: part - REAl(8):: xi(1:3) !Logical coordinates of particle in element - REAL(8):: BField(1:3) - - BField = 0.D0 - - xi = part%xi - - BField = mesh%vols(part%vol)%obj%gatherMF(xi) - - END FUNCTION gatherMagnField - !Assemble the source vector based on the charge density to solve Poisson's equation SUBROUTINE assembleSourceVector(vectorF) USE moduleMesh @@ -99,8 +65,8 @@ MODULE moduleEM !$OMP END SINGLE !$OMP DO REDUCTION(+:vectorF) - DO e = 1, mesh%numVols - nodes = mesh%vols(e)%obj%getNodes() + DO e = 1, mesh%numCells + nodes = mesh%cells(e)%obj%getNodes() nNodes = SIZE(nodes) !Calculates charge density (rho) in element nodes ALLOCATE(rho(1:nNodes)) @@ -113,7 +79,7 @@ MODULE moduleEM END DO !Calculates local F vector - localF = mesh%vols(e)%obj%elemF(rho) + localF = mesh%cells(e)%obj%elemF(rho) !Assign local F to global F DO i = 1, nNodes diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 135b08f..6962075 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -49,8 +49,8 @@ MODULE moduleSolver IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: volOld + CLASS(meshCell), POINTER, INTENT(inout):: volNew END SUBROUTINE weightingScheme_interface @@ -314,10 +314,10 @@ MODULE moduleSolver !$OMP SECTION !Erase the list of particles inside the cell if particles have been pushed DO s = 1, nSpecies - DO e = 1, mesh%numVols + DO e = 1, mesh%numCells IF (solver%pusher(s)%pushSpecies) THEN - CALL mesh%vols(e)%obj%listPart_in(s)%erase() - mesh%vols(e)%obj%totalWeight(s) = 0.D0 + CALL mesh%cells(e)%obj%listPart_in(s)%erase() + mesh%cells(e)%obj%totalWeight(s) = 0.D0 END IF @@ -328,10 +328,10 @@ MODULE moduleSolver !$OMP SECTION !Erase the list of particles inside the cell in coll mesh DO s = 1, nSpecies - DO e = 1, meshColl%numVols + DO e = 1, meshColl%numCells IF (solver%pusher(s)%pushSpecies) THEN - CALL meshColl%vols(e)%obj%listPart_in(s)%erase() - meshColl%vols(e)%obj%totalWeight(s) = 0.D0 + CALL meshColl%cells(e)%obj%listPart_in(s)%erase() + meshColl%cells(e)%obj%totalWeight(s) = 0.D0 END IF @@ -358,7 +358,7 @@ MODULE moduleSolver !Loops over the particles to scatter them !$OMP DO DO n = 1, nPartOld - CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) + CALL mesh%cells(partOld(n)%vol)%obj%scatter(partOld(n)) END DO !$OMP END DO @@ -383,8 +383,8 @@ MODULE moduleSolver IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: volOld + CLASS(meshCell), POINTER, INTENT(inout):: volNew REAL(8):: fractionVolume, pSplit !If particle changes volume to smaller cell @@ -416,7 +416,7 @@ MODULE moduleSolver TYPE(particle), INTENT(inout):: part INTEGER, INTENT(in):: nSplit - CLASS(meshVol), INTENT(inout):: vol + CLASS(meshCell), INTENT(inout):: vol REAL(8):: newWeight TYPE(particle), POINTER:: newPart INTEGER:: p @@ -454,15 +454,15 @@ MODULE moduleSolver CLASS(solverGeneric), INTENT(in):: self TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER:: volOld, volNew + CLASS(meshCell), POINTER:: volOld, volNew !Assume that particle is outside the domain part%n_in = .FALSE. - volOld => mesh%vols(part%vol)%obj + volOld => mesh%cells(part%vol)%obj CALL volOld%findCell(part) CALL findCellColl(part) - volNew => mesh%vols(part%vol)%obj + volNew => mesh%cells(part%vol)%obj !Call the NA shcme IF (ASSOCIATED(self%weightingScheme)) THEN CALL self%weightingScheme(part, volOld, volNew) diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90 index 69fcaab..c2aa46a 100644 --- a/src/modules/solver/pusher/modulePusher.f90 +++ b/src/modules/solver/pusher/modulePusher.f90 @@ -15,7 +15,7 @@ MODULE modulePusher PURE SUBROUTINE pushCartElectrostatic(part, tauIn) USE moduleSPecies - USE moduleEM + USE moduleMesh IMPLICIT NONE TYPE(particle), INTENT(inout):: part @@ -23,7 +23,8 @@ MODULE modulePusher REAL(8):: qmEFt(1:3) !Get the electric field at particle position - qmEFt = part%species%qm*gatherElecField(part)*tauIn + qmEFT = mesh%cells(part%vol)%obj%gatherElectricField(part%Xi) + qmEFt = qmEFt*part%species%qm*tauMin !Update velocity part%v = part%v + qmEFt @@ -34,8 +35,8 @@ MODULE modulePusher END SUBROUTINE pushCartElectrostatic PURE SUBROUTINE pushCartElectromagnetic(part, tauIn) - USE moduleSPecies - USE moduleEM + USE moduleSpecies + USE moduleMesh USE moduleMath IMPLICIT NONE @@ -49,13 +50,14 @@ MODULE modulePusher tauInHalf = tauIn *0.5D0 !Half of the force o f the electric field - qmEFt = part%species%qm*gatherElecField(part)*tauInHalf + qmEFT = mesh%cells(part%vol)%obj%gatherElectricField(part%Xi) + qmEFt = qmEFt*part%species%qm*tauInHalf !Half step for electrostatic v_minus = part%v + qmEFt !Full step rotation - B = gatherMagnField(part) + B = mesh%cells(part%vol)%obj%gatherMagneticField(part%Xi) BNorm = NORM2(B) IF (BNorm > 0.D0) THEN fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm @@ -112,7 +114,7 @@ MODULE modulePusher !Push one particle. Boris pusher for 2D Cyl Electrostatic particle PURE SUBROUTINE push2DCylElectrostatic(part, tauIn) USE moduleSpecies - USE moduleEM + USE moduleMesh IMPLICIT NONE TYPE(particle), INTENT(inout):: part @@ -124,7 +126,8 @@ MODULE modulePusher part_temp = part !Get electric field at particle position - qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn + qmEFT = mesh%cells(part_temp%vol)%obj%gatherElectricField(part_temp%Xi) + qmEFt = qmEFt*part_temp%species%qm*tauMin !z part_temp%v(1) = part%v(1) + qmEFt(1) part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn @@ -153,7 +156,6 @@ MODULE modulePusher !Push one particle. Boris pusher for 1D Radial Neutral particle PURE SUBROUTINE push1DRadNeutral(part, tauIn) USE moduleSpecies - USE moduleEM IMPLICIT NONE TYPE(particle), INTENT(inout):: part @@ -188,7 +190,7 @@ MODULE modulePusher !Push one particle. Boris pusher for 1D Radial Electrostatic particle PURE SUBROUTINE push1DRadElectrostatic(part, tauIn) USE moduleSpecies - USE moduleEM + USE moduleMesh IMPLICIT NONE TYPE(particle), INTENT(inout):: part @@ -200,7 +202,8 @@ MODULE modulePusher part_temp = part !Get electric field at particle position - qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin + qmEFT = mesh%cells(part_temp%vol)%obj%gatherElectricField(part_temp%Xi) + qmEFt = qmEFt*part_temp%species%qm*tauMin !r,theta v_p_oh_star(1) = part%v(1) + qmEFt(1) x_new = part%r(1) + v_p_oh_star(1)*tauIn @@ -226,7 +229,6 @@ MODULE modulePusher !Dummy pusher for 0D geometry PURE SUBROUTINE push0D(part, tauIn) USE moduleSpecies - USE moduleEM IMPLICIT NONE TYPE(particle), INTENT(inout):: part