diff --git a/runs/1D_Cathode/inputRad.json b/runs/1D_Cathode/inputRad.json index 34a09ce..824c0a8 100644 --- a/runs/1D_Cathode/inputRad.json +++ b/runs/1D_Cathode/inputRad.json @@ -31,10 +31,10 @@ ]} ], "boundaryEM": [ - {"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} + {"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} ], "inject": [ - {"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-2, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0], + {"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-5, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0], "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1} ], "solver": { diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 1224f53..3a8a033 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -74,7 +74,10 @@ PROGRAM fpakc tColl = omp_get_wtime() !$OMP END SINGLE - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions(t) + IF (doMCC) THEN + CALL meshForMCC%doCollisions(t) + + END IF !$OMP SINGLE tColl = omp_get_wtime() - tColl @@ -83,7 +86,10 @@ PROGRAM fpakc tCoul = omp_get_wTime() !$OMP END SINGLE - IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb() + IF (ASSOCIATED(mesh%doCoulomb)) THEN + CALL mesh%doCoulomb() + + END IF !$OMP SINGLE tCoul = omp_get_wTime() - tCoul diff --git a/src/modules/common/moduleConstParam.f90 b/src/modules/common/moduleConstParam.f90 index cb12a23..58a3cc8 100644 --- a/src/modules/common/moduleConstParam.f90 +++ b/src/modules/common/moduleConstParam.f90 @@ -6,7 +6,8 @@ MODULE moduleConstParam REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi REAL(8), PARAMETER:: PI2 = 2.D0*PI !2*pi - REAL(8), PARAMETER:: PI8 = 8.D0*PI !2*pi + REAL(8), PARAMETER:: PI4 = 4.D0*PI !4*pi + REAL(8), PARAMETER:: PI8 = 8.D0*PI !8*pi REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1 REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index ac40c3b..81ababd 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -331,13 +331,14 @@ MODULE moduleInput REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity INTEGER, ALLOCATABLE, DIMENSION(:):: nodes INTEGER:: nNodes - REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi + REAL(8), ALLOCATABLE, DIMENSION(:):: sourceScalar + REAL(8), ALLOCATABLE, DIMENSION(:,:):: sourceArray !Density at the volume centroid REAL(8):: densityCen !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi INTEGER:: nNewPart = 0.D0 - CLASS(meshVol), POINTER:: vol + CLASS(meshCell), POINTER:: cell TYPE(particle), POINTER:: partNew REAL(8):: vTh TYPE(lNode), POINTER:: partCurr, partNext @@ -354,87 +355,74 @@ MODULE moduleInput CALL config%get(object // '.file', spFile, found) !Reads node values at the nodes filename = path // spFile - CALL mesh%readInitial(sp, filename, density, velocity, temperature) + CALL mesh%readInitial(filename, density, velocity, temperature) !For each volume in the node, create corresponding particles - DO e = 1, mesh%numVols + DO e = 1, mesh%numCells !Scale variables !Density at centroid of cell - nodes = mesh%vols(e)%obj%getNodes() - nNodes = SIZE(nodes) - fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) - ALLOCATE(source(1:nNodes)) - DO j = 1, nNodes - source(j) = density(nodes(j)) - - END DO - densityCen = DOT_PRODUCT(fPsi, source) - DEALLOCATE(fPsi) + nNodes = mesh%cells(e)%obj%nNodes + nodes = mesh%cells(e)%obj%getNodes(nNodes) + ALLOCATE(sourceScalar(1:nNodes)) + ALLOCATE(sourceArray(1:nNodes, 1:3)) + sourceScalar = density(nodes) + densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar) !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 - fPsi = mesh%vols(e)%obj%fPsi(partNew%xi) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 1) - - END DO - velocityXi(1) = DOT_PRODUCT(fPsi, source) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 2) - - END DO - velocityXi(2) = DOT_PRODUCT(fPsi, source) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 3) - - END DO - velocityXi(3) = DOT_PRODUCT(fPsi, source) + sourceArray(:,1) = velocity((nodes), 1) + sourceArray(:,2) = velocity((nodes), 2) + sourceArray(:,3) = velocity((nodes), 3) + velocityXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceArray) velocityXi = velocityXi / v_ref !Get temperature at particle position - DO j = 1, nNodes - source(j) = temperature(nodes(j)) - - END DO - temperatureXi = DOT_PRODUCT(fPsi, source) + sourceScalar = temperature(nodes) + temperatureXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceScalar) temperatureXi = temperatureXi / T_ref vTh = DSQRT(temperatureXi / species(sp)%obj%m) partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian() partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() - partNew%vol = e + + partNew%cell = e IF (doubleMesh) THEN - partNew%volColl = findCellBrute(meshColl, partNew%r) + partNew%cellColl = findCellBrute(meshColl, partNew%r) ELSE - partNew%volColl = partNew%vol + partNew%cellColl = partNew%cell END IF + partNew%n_in = .TRUE. + partNew%weight = species(sp)%obj%weight !Assign particle to temporal list of particles CALL partInitial%add(partNew) !Assign particle to list in volume - vol => meshforMCC%vols(partNew%volColl)%obj - CALL OMP_SET_LOCK(vol%lock) - CALL vol%listPart_in(sp)%add(partNew) - vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight - CALL OMP_UNSET_LOCK(vol%lock) + IF (doMCC) THEN + cell => meshforMCC%cells(partNew%cellColl)%obj + CALL OMP_SET_LOCK(cell%lock) + CALL cell%listPart_in(sp)%add(partNew) + cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight + CALL OMP_UNSET_LOCK(cell%lock) + + END IF END DO - DEALLOCATE(source) + DEALLOCATE(sourceScalar, sourceArray) END DO @@ -643,9 +631,9 @@ MODULE moduleInput REAL(8):: energyThreshold, energyBinding CHARACTER(:), ALLOCATABLE:: electron INTEGER:: e - CLASS(meshVol), POINTER:: vol + CLASS(meshCell), POINTER:: cell - !Firstly, checks if the object 'interactions' exists + !Firstly, check if the object 'interactions' exists CALL config%info('interactions', found) IF (found) THEN !Checks if MC collisions have been defined @@ -739,18 +727,18 @@ 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 + cell => meshForMCC%cells(e)%obj !Allocate Maximum cross section per collision pair and assign the initial collision rate - ALLOCATE(vol%sigmaVrelMax(1:nCollPairs)) - vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) + ALLOCATE(cell%sigmaVrelMax(1:nCollPairs)) + cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) IF (collOutput) THEN - ALLOCATE(vol%tallyColl(1:nCollPairs)) + ALLOCATE(cell%tallyColl(1:nCollPairs)) DO k = 1, nCollPairs - ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) - vol%tallyColl(k)%tally = 0 + ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) + cell%tallyColl(k)%tally = 0 END DO @@ -907,6 +895,8 @@ MODULE moduleInput END IF + doMCC = ASSOCIATED(meshForMCC) + !Get the dimension of the geometry CALL config%get(object // '.dimension', mesh%dimen, found) IF (.NOT. found) THEN @@ -930,8 +920,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 8e8072d..133ae8e 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -6,27 +6,32 @@ MODULE moduleMesh0D TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D INTEGER:: n1 CONTAINS - PROCEDURE, PASS:: init => initNode0D + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode0D PROCEDURE, PASS:: getCoordinates => getCoord0D END TYPE meshNode0D - TYPE, PUBLIC, EXTENDS(meshVol):: meshVol0D + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell0D CLASS(meshNode), POINTER:: n1 CONTAINS - PROCEDURE, PASS:: init => initVol0D - PROCEDURE, PASS:: getNodes => getNodes0D - PROCEDURE, PASS:: randPos => randPos0D - PROCEDURE, NOPASS:: fPsi => fPsi0D - PROCEDURE, PASS:: gatherEF => gatherEF0D - PROCEDURE, PASS:: gatherMF => gatherMF0D - PROCEDURE, PASS:: elemK => elemK0D - PROCEDURE, PASS:: elemF => elemF0D - PROCEDURE, PASS:: phy2log => phy2log0D - PROCEDURE, NOPASS:: inside => inside0D - PROCEDURE, PASS:: nextElement => nextElement0D + PROCEDURE, PASS:: init => initCell0D + PROCEDURE, PASS:: getNodes => getNodes0D + PROCEDURE, PASS:: randPos => randPos0D + PROCEDURE, NOPASS:: fPsi => fPsi0D + PROCEDURE, NOPASS:: dPsi => dPsi0D + PROCEDURE, PASS:: partialDer => partialDer0D + PROCEDURE, NOPASS:: detJac => detJ0D + PROCEDURE, NOPASS:: invJac => invJ0D + PROCEDURE, PASS:: gatherElectricField => gatherEF0D + PROCEDURE, PASS:: gatherMagneticField => gatherMF0D + PROCEDURE, PASS:: elemK => elemK0D + PROCEDURE, PASS:: elemF => elemF0D + PROCEDURE, PASS:: phy2log => phy2log0D + PROCEDURE, NOPASS:: inside => inside0D + PROCEDURE, PASS:: neighbourElement => neighbourElement0D - END TYPE meshVol0D + END TYPE meshCell0D CONTAINS !NODE FUNCTIONS @@ -61,18 +66,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,88 +89,122 @@ 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) + !Get the nodes of the volume + PURE FUNCTION getNodes0D(self, nNodes) RESULT(n) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:1)) + CLASS(meshCell0D), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = self%n1%n END FUNCTION getNodes0D + !Calculate random position inside the volume 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 FUNCTION fPsi0D(xi) RESULT(fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + PURE FUNCTION fPsi0D(Xi, nNodes) RESULT(fPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - ALLOCATE(fPsi(1:1)) fPsi = 1.D0 END FUNCTION fPsi0D - PURE FUNCTION gatherEF0D(self, xi) RESULT(EF) + PURE FUNCTION dPsi0D(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: EF(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - EF = 0.D0 + dPsi = 0.D0 - END FUNCTION gatherEF0D + END FUNCTION dPsi0D - PURE FUNCTION gatherMF0D(self, xi) RESULT(MF) + PURE FUNCTION partialDer0D(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: MF(1:3) + CLASS(meshCell0D), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - MF = 0.D0 + pDer = 0.D0 - END FUNCTION gatherMF0D + END FUNCTION partialDer0D - PURE FUNCTION elemK0D(self) RESULT(localK) + PURE FUNCTION elemK0D(self, nNodes) RESULT(localK) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + CLASS(meshCell0D), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) - ALLOCATE(localK(1:1, 1:1)) localK = 0.D0 END FUNCTION elemK0D - PURE FUNCTION elemF0D(self, source) RESULT(localF) + PURE FUNCTION elemF0D(self, nNodes, 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1: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, 1, 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, 1, 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) @@ -171,25 +212,46 @@ 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 neighbourElement0D(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol0D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshCell0D), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - nextElement => NULL() + neighbourElement => NULL() - END SUBROUTINE nextElement0D + END SUBROUTINE neighbourElement0D + + !COMMON FUNCTIONS + PURE FUNCTION detJ0D(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: dJ + + dJ = 0.D0 + + END FUNCTION detJ0D + + PURE FUNCTION invJ0D(pDer) RESULT(invJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) + + invJ = 0.D0 + + END FUNCTION invJ0D END MODULE moduleMesh0D diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 0a846d9..269f157 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -8,13 +8,14 @@ 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 REAL(8):: x = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode1DCart + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode1DCart PROCEDURE, PASS:: getCoordinates => getCoord1DCart END TYPE meshNode1DCart @@ -25,6 +26,7 @@ MODULE moduleMesh1DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge1DCart PROCEDURE, PASS:: getNodes => getNodes1DCart PROCEDURE, PASS:: intersection => intersection1DCart @@ -32,57 +34,34 @@ MODULE moduleMesh1DCart END TYPE meshEdge1DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ1DCart - PROCEDURE, PASS:: invJac => invJ1DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol1DCart - - 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:) - REAL(8), INTENT(out), DIMENSION(1):: dx - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshVol1DCart):: meshVol1DCartSegm + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DCartSegm !Element coordinates REAL(8):: x(1:2) !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() - REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initVol1DCartSegm - PROCEDURE, PASS:: randPos => randPos1DCartSeg - PROCEDURE, PASS:: area => areaSegm - PROCEDURE, NOPASS:: fPsi => fPsiSegm - PROCEDURE, NOPASS:: dPsi => dPsiSegm - PROCEDURE, PASS:: partialDer => partialDerSegm - PROCEDURE, PASS:: elemK => elemKSegm - PROCEDURE, PASS:: elemF => elemFSegm - PROCEDURE, NOPASS:: inside => insideSegm - PROCEDURE, PASS:: gatherEF => gatherEFSegm - PROCEDURE, PASS:: gatherMF => gatherMFSegm - PROCEDURE, PASS:: getNodes => getNodesSegm - PROCEDURE, PASS:: phy2log => phy2logSegm - PROCEDURE, PASS:: nextElement => nextElementSegm + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCell1DCartSegm + PROCEDURE, PASS:: getNodes => getNodesSegm + PROCEDURE, PASS:: randPos => randPos1DCartSegm + PROCEDURE, NOPASS:: fPsi => fPsiSegm + PROCEDURE, NOPASS:: dPsi => dPsiSegm + PROCEDURE, PASS:: partialDer => partialDerSegm + PROCEDURE, NOPASS:: detJac => detJ1DCart + PROCEDURE, NOPASS:: invJac => invJ1DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFSegm + PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm + PROCEDURE, PASS:: elemK => elemKSegm + PROCEDURE, PASS:: elemF => elemFSegm + PROCEDURE, NOPASS:: inside => insideSegm + PROCEDURE, PASS:: phy2log => phy2logSegm + PROCEDURE, PASS:: neighbourElement => neighbourElementSegm + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm - END TYPE meshVol1DCartSegm + END TYPE meshCell1DCartSegm CONTAINS !NODE FUNCTIONS @@ -120,7 +99,7 @@ MODULE moduleMesh1DCart END FUNCTION getCoord1DCart !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -136,6 +115,7 @@ MODULE moduleMesh1DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj !Get element coordinates r1 = self%n1%getCoordinates() @@ -152,20 +132,20 @@ MODULE moduleMesh1DCart CALL pointBoundaryFunction(self, s) END DO - + !Physical Surface self%physicalSurface = physicalSurface END SUBROUTINE initEdge1DCart !Get nodes from edge - PURE FUNCTION getNodes1DCart(self) RESULT(n) + PURE FUNCTION getNodes1DCart(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1)) n = (/ self%n1%n /) END FUNCTION getNodes1DCart @@ -181,7 +161,7 @@ MODULE moduleMesh1DCart END FUNCTION intersection1DCart - !Calculates a 'random' position in edge + !Calculate a 'random' position in edge FUNCTION randPosEdge(self) RESULT(r) CLASS(meshEdge1DCart), INTENT(in):: self REAL(8):: r(1:3) @@ -192,18 +172,19 @@ MODULE moduleMesh1DCart !VOLUME FUNCTIONS !SEGMENT FUNCTIONS - !Init segment element - SUBROUTINE initVol1DCartSegm(self, n, p, nodes) + !Init element + 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 @@ -212,144 +193,182 @@ MODULE moduleMesh1DCart self%x = (/ r1(1), r2(1) /) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) 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) + !Get nodes from 1D volume + PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/ self%n1%n, self%n2%n /) + + END FUNCTION getNodesSegm + + !Random position in 1D volume + 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):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:2) - xii(1) = random(-1.D0, 1.D0) - xii(2:3) = 0.D0 + Xi = 0.D0 + Xi(1) = random(-1.D0, 1.D0) - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 2) + + r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%x) - END FUNCTION randPos1DCartSeg + END FUNCTION randPos1DCartSegm - !Computes element area - PURE SUBROUTINE areaSegm(self) + !Compute element functions at point Xi + PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(inout):: self - REAL(8):: l !element length - REAL(8):: fPsi(1:2) - REAL(8):: detJ - REAL(8):: Xii(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - self%volume = 0.D0 - self%arNodes = 0.D0 - !1 point Gauss integral - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - l = 2.D0*detJ - self%volume = l - self%arNodes = fPsi*l + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) - END SUBROUTINE areaSegm - - !Computes element functions at point xii - PURE FUNCTION fPsiSegm(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:2)) - - fPsi(1) = 1.D0 - xi(1) - fPsi(2) = 1.D0 + xi(1) - fPsi = fPsi * 5.D-1 + fPsi = fPsi * 0.50D0 END FUNCTION fPsiSegm - !Computes element derivative shape function at Xii - PURE FUNCTION dPsiSegm(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - ALLOCATE(dPsi(1:1, 1:2)) + dPsi = 0.D0 - dPsi(1, 1) = -5.D-1 - dPsi(1, 2) = 5.D-1 + dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /) END FUNCTION dPsiSegm - !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerSegm(self, dPsi, dx) + !Partial derivative in global coordinates + PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1):: dx + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) + pDer = 0.D0 - END SUBROUTINE partialDerSegm + pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%x(1:2)) + pDer(2,2) = 1.D0 + pDer(3,3) = 1.D0 - !Computes local stiffness matrix - PURE FUNCTION elemKSegm(self) RESULT(localK) + END FUNCTION partialDerSegm + + 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, 2, 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, 2, B) + + END FUNCTION gatherMFSegm + + !Compute element local stiffness matrix + PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: Xii(1:3) - REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ(1), detJ + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ INTEGER:: l - ALLOCATE(localK(1:2,1:2)) - localK = 0.D0 - Xii = 0.D0 + localK = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral DO l = 1, 3 - xii(1) = corSeg(l) - dPsi = self%dPsi(Xii) - detJ = self%detJac(Xii, dPsi) - invJ = self%invJac(Xii, dPsi) - localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - wSeg(l)/detJ + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + wSeg(l)/detJ END DO END FUNCTION elemKSegm - PURE FUNCTION elemFSegm(self, source) RESULT(localF) + !Compute the local source vector for a force f + PURE FUNCTION elemFSegm(self, nNodes, 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:2) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ, f - REAL(8):: Xii(1:3) INTEGER:: l - ALLOCATE(localF(1:2)) localF = 0.D0 - Xii = 0.D0 + Xi = 0.D0 + !Start 1D Gauss Quad Integral DO l = 1, 3 - xii(1) = corSeg(l) - detJ = self%detJac(Xii) - fPsi = self%fPsi(Xii) + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wSeg(l)*detJ @@ -357,151 +376,98 @@ MODULE moduleMesh1DCart END FUNCTION elemFSegm - PURE FUNCTION insideSegm(xi) RESULT(ins) + PURE FUNCTION insideSegm(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >=-1.D0 .AND. & - xi(1) <= 1.D0 + ins = Xi(1) >=-1.D0 .AND. & + Xi(1) <= 1.D0 END FUNCTION insideSegm - !Gathers EF at position Xii - PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF) + PURE FUNCTION phy2logSegm(self, r) RESULT(Xi) 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) /) - - fPsi = self%fPsi(xi) - 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(:) - - ALLOCATE(n(1:2)) - n = (/ self%n1%n, self%n2%n /) - - END FUNCTION getNodesSegm - - 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) + REAL(8):: Xi(1:3) - xN = 0.D0 - xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0 + Xi = 0.D0 + + Xi(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0 END FUNCTION phy2logSegm - !Get next element for a logical position xi - SUBROUTINE nextElementSegm(self, xi, nextElement) + !Get the next element for a logical position Xi + SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshCell1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - NULLIFY(nextElement) - IF (xi(1) < -1.D0) THEN - nextElement => self%e2 + NULLIFY(neighbourElement) + IF (Xi(1) < -1.D0) THEN + neighbourElement => self%e2 - ELSEIF (xi(1) > 1.D0) THEN - nextElement => self%e1 + ELSEIF (Xi(1) > 1.D0) THEN + neighbourElement => self%e1 END IF - END SUBROUTINE nextElementSegm + END SUBROUTINE neighbourElementSegm + + !Compute element volume + PURE SUBROUTINE volumeSegm(self) + IMPLICIT NONE + + CLASS(meshCell1DCartSegm), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:2) + + self%volume = 0.D0 + !1D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + !Compute total volume of the cell + self%volume = detJ*2.D0 + !Compute volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + + END SUBROUTINE volumeSegm !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Calculates a random position in 1D volume - !Computes the element Jacobian determinant - PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ) + !Compute element Jacobian determinant + PURE FUNCTION detJ1DCart(pDer) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol1DCart), 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):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dx(1) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - dJ = dx(1) + dJ = pDer(1, 1) END FUNCTION detJ1DCart - !Computes the invers Jacobian - PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ) + !Compute element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ1DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol1DCart), 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):: dx(1) - REAL(8):: invJ + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in + invJ = 0.D0 - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - invJ = 1.D0/dx(1) + invJ(1, 1) = 1.D0/pDer(1, 1) + invJ(2, 2) = 1.D0 + invJ(3, 3) = 1.D0 END FUNCTION invJ1DCart @@ -511,11 +477,11 @@ MODULE moduleMesh1DCart 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 @@ -523,9 +489,9 @@ MODULE moduleMesh1DCart 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 @@ -535,29 +501,29 @@ MODULE moduleMesh1DCart END SUBROUTINE connectMesh1DCart - 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(meshVol1DCartSegm) + TYPE IS(meshCell1DCartSegm) SELECT TYPE(elemB) - TYPE IS(meshVol1DCartSegm) + TYPE IS(meshCell1DCartSegm) CALL connectSegmSegm(elemA, elemB) END SELECT END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell 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 @@ -576,14 +542,14 @@ MODULE moduleMesh1DCart END SUBROUTINE connectSegmSegm - 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(elemA) - TYPE IS (meshVol1DCartSegm) + TYPE IS (meshCell1DCartSegm) SELECT TYPE(elemB) CLASS IS(meshEdge1DCart) CALL connectSegmEdge(elemA, elemB) @@ -592,12 +558,12 @@ MODULE moduleMesh1DCart END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge 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. & @@ -606,7 +572,7 @@ MODULE moduleMesh1DCart elemA%e1 => elemB elemB%e2 => elemA - !Revers the normal to point inside the domain + !Rever the normal to point inside the domain elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index b8edfdd..d998267 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -8,13 +8,14 @@ MODULE moduleMesh1DRad 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):: meshNode1DRad !Element coordinates REAL(8):: r = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode1DRad + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode1DRad PROCEDURE, PASS:: getCoordinates => getCoord1DRad END TYPE meshNode1DRad @@ -25,6 +26,7 @@ MODULE moduleMesh1DRad !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge1DRad PROCEDURE, PASS:: getNodes => getNodes1DRad PROCEDURE, PASS:: intersection => intersection1DRad @@ -32,58 +34,34 @@ MODULE moduleMesh1DRad END TYPE meshEdge1DRad - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad - CONTAINS - PROCEDURE, PASS:: detJac => detJ1DRad - PROCEDURE, PASS:: invJac => invJ1DRad - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol1DRad - - 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 - - CLASS(meshVol1DRad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1):: dx - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshVol1DRad):: meshVol1DRadSegm + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DRadSegm !Element coordinates REAL(8):: r(1:2) !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() - REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initVol1DRadSegm - PROCEDURE, PASS:: randPos => randPos1DRadSeg - PROCEDURE, PASS:: area => areaRad - PROCEDURE, NOPASS:: fPsi => fPsiRad - PROCEDURE, NOPASS:: dPsi => dPsiRad - PROCEDURE, PASS:: partialDer => partialDerRad - PROCEDURE, PASS:: elemK => elemKRad - PROCEDURE, PASS:: elemF => elemFRad - PROCEDURE, NOPASS:: inside => insideRad - PROCEDURE, PASS:: gatherEF => gatherEFRad - PROCEDURE, PASS:: gatherMF => gatherMFRad - PROCEDURE, PASS:: getNodes => getNodesRad - PROCEDURE, PASS:: phy2log => phy2logRad - PROCEDURE, PASS:: nextElement => nextElementRad + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCell1DRadSegm + PROCEDURE, PASS:: getNodes => getNodesSegm + PROCEDURE, PASS:: randPos => randPos1DRadSegm + PROCEDURE, NOPASS:: fPsi => fPsiSegm + PROCEDURE, NOPASS:: dPsi => dPsiSegm + PROCEDURE, PASS:: partialDer => partialDerSegm + PROCEDURE, NOPASS:: detJac => detJ1DRad + PROCEDURE, NOPASS:: invJac => invJ1DRad + PROCEDURE, PASS:: gatherElectricField => gatherEFSegm + PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm + PROCEDURE, PASS:: elemK => elemKSegm + PROCEDURE, PASS:: elemF => elemFSegm + PROCEDURE, NOPASS:: inside => insideSegm + PROCEDURE, PASS:: phy2log => phy2logSegm + PROCEDURE, PASS:: neighbourElement => neighbourElementSegm + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm - END TYPE meshVol1DRadSegm + END TYPE meshCell1DRadSegm CONTAINS !NODE FUNCTIONS @@ -103,7 +81,7 @@ MODULE moduleMesh1DRad !Node volume, to be determined in mesh self%v = 0.D0 - !Allocates output + !Allocate output ALLOCATE(self%output(1:nSpecies)) CALL OMP_INIT_LOCK(self%lock) @@ -121,7 +99,7 @@ MODULE moduleMesh1DRad END FUNCTION getCoord1DRad !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -137,6 +115,7 @@ MODULE moduleMesh1DRad INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj !Get element coordinates r1 = self%n1%getCoordinates() @@ -144,7 +123,6 @@ MODULE moduleMesh1DRad self%r = r1(1) self%normal = (/ 1.D0, 0.D0, 0.D0 /) - self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -161,13 +139,13 @@ MODULE moduleMesh1DRad END SUBROUTINE initEdge1DRad !Get nodes from edge - PURE FUNCTION getNodes1DRad(self) RESULT(n) + PURE FUNCTION getNodes1DRad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1)) n = (/ self%n1%n /) END FUNCTION getNodes1DRad @@ -183,7 +161,7 @@ MODULE moduleMesh1DRad END FUNCTION intersection1DRad - !Calculates a 'random' position in edge + !Calculate a 'random' position in edge FUNCTION randPos1DRad(self) RESULT(r) CLASS(meshEdge1DRad), INTENT(in):: self REAL(8):: r(1:3) @@ -194,18 +172,19 @@ MODULE moduleMesh1DRad !VOLUME FUNCTIONS !SEGMENT FUNCTIONS - !Init segment element - SUBROUTINE initVol1DRadSegm(self, n, p, nodes) + !Init element + 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 @@ -214,312 +193,296 @@ MODULE moduleMesh1DRad self%r = (/ r1(1), r2(1) /) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) 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) + !Get nodes from 1D volume + PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/ self%n1%n, self%n2%n /) + + END FUNCTION getNodesSegm + + !Random position in 1D volume + FUNCTION randPos1DRadSegm(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self + CLASS(meshCell1DRadSegm), 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:2) - xii(1) = random(-1.D0, 1.D0) - xii(2:3) = 0.D0 + Xi = 0.D0 + Xi(1) = random(-1.D0, 1.D0) - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 2) + + r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%r) - END FUNCTION randPos1DRadSeg + END FUNCTION randPos1DRadSegm - !Computes element area - PURE SUBROUTINE areaRad(self) + !Compute element functions at point Xi + PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(inout):: self - REAL(8):: l !element length - REAL(8):: fPsi(1:2) - REAL(8):: r - REAL(8):: detJ - REAL(8):: Xii(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - self%volume = 0.D0 - self%arNodes = 0.D0 - !1 point Gauss integral - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - !Computes total volume of the cell - r = DOT_PRODUCT(fPsi, self%r) - l = 2.D0*detJ - self%volume = r*l - !Computes volume per node - Xii = (/-5.D-1, 0.D0, 0.D0/) - r = DOT_PRODUCT(self%fPsi(Xii),self%r) - self%arNodes(1) = fPsi(1)*r*l - Xii = (/ 5.D-1, 0.D0, 0.D0/) - r = DOT_PRODUCT(self%fPsi(Xii),self%r) - self%arNodes(2) = fPsi(2)*r*l + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) - END SUBROUTINE areaRad + fPsi = fPsi * 0.50D0 - !Computes element functions at point xii - PURE FUNCTION fPsiRad(xi) RESULT(fPsi) + END FUNCTION fPsiSegm + + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - ALLOCATE(fPsi(1:2)) + dPsi = 0.D0 - fPsi(1) = 1.D0 - xi(1) - fPsi(2) = 1.D0 + xi(1) - fPsi = fPsi * 5.D-1 + dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /) - END FUNCTION fPsiRad + END FUNCTION dPsiSegm - !Computes element derivative shape function at Xii - PURE FUNCTION dPsiRad(xi) RESULT(dPsi) + !Partial derivative in global coordinates + PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - ALLOCATE(dPsi(1:1, 1:2)) + pDer = 0.D0 - dPsi(1, 1) = -5.D-1 - dPsi(1, 2) = 5.D-1 + pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%r(1:2)) + pDer(2,2) = 1.D0 + pDer(3,3) = 1.D0 - END FUNCTION dPsiRad + END FUNCTION partialDerSegm - !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerRad(self, dPsi, dx) + PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array) IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1):: dx - - dx(1) = DOT_PRODUCT(dPsi(1,:), self%r) - - END SUBROUTINE partialDerRad - - !Computes local stiffness matrix - PURE FUNCTION elemKRad(self) RESULT(localK) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: Xii(1:3) - REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ(1), detJ - REAL(8):: r, fPsi(1:2) - INTEGER:: l - - ALLOCATE(localK(1:2, 1:2)) - localK = 0.D0 - Xii = 0.D0 - DO l = 1, 3 - xii(1) = corSeg(l) - dPsi = self%dPsi(Xii) - detJ = self%detJac(Xii, dPsi) - invJ = self%invJac(Xii, dPsi) - fPsi = self%fPsi(Xii) - r = DOT_PRODUCT(fPsi, self%r) - localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - r*wSeg(l)/detJ - - END DO - - localK = localK*PI2 - - END FUNCTION elemKRad - - PURE FUNCTION elemFRad(self, source) RESULT(localF) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: fPsi(1:2) - REAL(8):: detJ, f, r - REAL(8):: Xii(1:3) - INTEGER:: l - - ALLOCATE(localF(1:2)) - localF = 0.D0 - Xii = 0.D0 - - DO l = 1, 3 - xii(1) = corSeg(l) - detJ = self%detJac(Xii) - fPsi = self%fPsi(Xii) - r = DOT_PRODUCT(fPsi, self%r) - f = DOT_PRODUCT(fPsi, source) - localF = localF + f*fPsi*r*wSeg(l)*detJ - - END DO - - END FUNCTION elemFRad - - PURE FUNCTION insideRad(xi) RESULT(ins) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = xi(1) >=-1.D0 .AND. & - xi(1) <= 1.D0 - - END FUNCTION insideRad - - !Gathers EF at position Xii - 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) + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: array(1:3) 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 + array = -self%gatherDF(Xi, 2, phi) - END FUNCTION gatherEFRad + END FUNCTION gatherEFSegm - PURE FUNCTION gatherMFRad(self, xi) RESULT(MF) + PURE FUNCTION gatherMFSegm(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, 2, B) + + END FUNCTION gatherMFSegm + + !Compute element local stiffness matrix + PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) + USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: r + REAL(8):: invJ(1:3, 1:3), detJ + INTEGER:: l + + localK = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + r = self%gatherF(Xi, 2, self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + r*wSeg(l)/detJ + + END DO + localK = localK*PI2 + + END FUNCTION elemKSegm + + !Compute the local source vector for a force f + PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:2) - REAL(8):: MF_Nodes(1:2, 1:3) - REAL(8):: MF(1:3) - REAL(8):: invJ + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) + REAL(8):: detJ, f + REAL(8):: r + INTEGER:: l - 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) /) + localF = 0.D0 - fPsi = self%fPsi(xi) - MF = MATMUL(fPsi, MF_Nodes) + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) + localF = localF + r*f*fPsi*wSeg(l)*detJ - END FUNCTION gatherMFRad + END DO + localF = localF*PI2 - !Get nodes from 1D volume - PURE FUNCTION getNodesRad(self) RESULT(n) + END FUNCTION elemFSegm + + PURE FUNCTION insideSegm(Xi) RESULT(ins) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + REAL(8), INTENT(in):: Xi(1:3) + LOGICAL:: ins - ALLOCATE(n(1:2)) - n = (/ self%n1%n, self%n2%n /) + ins = Xi(1) >=-1.D0 .AND. & + Xi(1) <= 1.D0 - END FUNCTION getNodesRad + END FUNCTION insideSegm - PURE FUNCTION phy2logRad(self, r) RESULT(xN) + PURE FUNCTION phy2logSegm(self, r) RESULT(Xi) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self + CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) + REAL(8):: Xi(1:3) - xN = 0.D0 - xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0 + Xi = 0.D0 - END FUNCTION phy2logRad + Xi(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0 - !Get next element for a logical position xi - SUBROUTINE nextElementRad(self, xi, nextElement) + END FUNCTION phy2logSegm + + !Get the next element for a logical position Xi + SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshCell1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - NULLIFY(nextElement) - IF (xi(1) < -1.D0) THEN - nextElement => self%e2 + NULLIFY(neighbourElement) + IF (Xi(1) < -1.D0) THEN + neighbourElement => self%e2 - ELSEIF (xi(1) > 1.D0) THEN - nextElement => self%e1 + ELSEIF (Xi(1) > 1.D0) THEN + neighbourElement => self%e1 END IF - END SUBROUTINE nextElementRad + END SUBROUTINE neighbourElementSegm + + !Compute element volume + PURE SUBROUTINE volumeSegm(self) + USE moduleConstParam, ONLY: PI4 + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:2) + REAL(8):: r + + self%volume = 0.D0 + !1D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + r = DOT_PRODUCT(fPsi, self%r) + !Compute total volume of the cell + self%volume = r*detJ*PI4 !2*2PI + !Compute volume per node + Xi = (/-5.D-1, 0.D0, 0.D0/) + r = self%gatherF(Xi, 2, self%r) + self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI4 + Xi = (/ 5.D-1, 0.D0, 0.D0/) + r = self%gatherF(Xi, 2, self%r) + self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI4 + + END SUBROUTINE volumeSegm !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Computes the element Jacobian determinant - PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ) + !Compute element Jacobian determinant + PURE FUNCTION detJ1DRad(pDer) RESULT(dJ) IMPLICIT NONE - CLASS(meshVol1DRad), 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):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dx(1) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - dJ = dx(1) + dJ = pDer(1, 1) END FUNCTION detJ1DRad - !Computes the invers Jacobian - PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ) + !Compute element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ1DRad(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol1DRad), 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):: dx(1) - REAL(8):: invJ + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in + invJ = 0.D0 - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - invJ = 1.D0/dx(1) + invJ(1, 1) = 1.D0/pDer(1, 1) + invJ(2, 2) = 1.D0 + invJ(3, 3) = 1.D0 END FUNCTION invJ1DRad @@ -529,11 +492,11 @@ MODULE moduleMesh1DRad 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 @@ -541,9 +504,9 @@ MODULE moduleMesh1DRad 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 @@ -553,29 +516,29 @@ MODULE moduleMesh1DRad END SUBROUTINE connectMesh1DRad - 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(meshVol1DRadSegm) + TYPE IS(meshCell1DRadSegm) SELECT TYPE(elemB) - TYPE IS(meshVol1DRadSegm) + TYPE IS(meshCell1DRadSegm) CALL connectSegmSegm(elemA, elemB) END SELECT END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell 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 @@ -594,14 +557,14 @@ MODULE moduleMesh1DRad END SUBROUTINE connectSegmSegm - 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(elemA) - TYPE IS (meshVol1DRadSegm) + TYPE IS (meshCell1DRadSegm) SELECT TYPE(elemB) CLASS IS(meshEdge1DRad) CALL connectSegmEdge(elemA, elemB) @@ -610,12 +573,12 @@ MODULE moduleMesh1DRad END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge 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. & @@ -624,7 +587,7 @@ MODULE moduleMesh1DRad elemA%e1 => elemB elemB%e2 => elemA - !Revers the normal to point inside the domain + !Rever the normal to point inside the domain elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index d01c0e4..c02078a 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -11,15 +11,16 @@ MODULE moduleMesh2DCart 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):: meshNode2DCart !Element coordinates REAL(8):: x = 0.D0, y = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode2DCart + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode2DCart PROCEDURE, PASS:: getCoordinates => getCoord2DCart END TYPE meshNode2DCart @@ -30,97 +31,75 @@ MODULE moduleMesh2DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdge2DCart - PROCEDURE, PASS:: getNodes => getNodes2DCart + !meshEdge DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initEdge2DCart + PROCEDURE, PASS:: getNodes => getNodes2DCart PROCEDURE, PASS:: intersection => intersection2DCartEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdge2DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ2DCart - PROCEDURE, PASS:: invJac => invJ2DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol2DCart - - 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:) - REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy - - END SUBROUTINE partialDer_interface - - END INTERFACE - !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartQuad + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartQuad !Element coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0 !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() - REAL(8):: arNodes(1:4) = 0.D0 - CONTAINS - PROCEDURE, PASS:: init => initVolQuad2DCart - 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:: elemK => elemKQuad - PROCEDURE, PASS:: elemF => elemFQuad - 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 + CONTAINS + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellQuad2DCart + PROCEDURE, PASS:: getNodes => getNodesQuad + PROCEDURE, PASS:: randPos => randPosCellQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, NOPASS:: detJac => detJ2DCart + PROCEDURE, NOPASS:: invJac => invJ2DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFQuad + PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad + PROCEDURE, PASS:: elemK => elemKQuad + PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: neighbourElement => neighbourElementQuad + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad + + END TYPE meshCell2DCartQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartTria + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartTria !Element coordinates REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0 !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() - REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolTria2DCart - PROCEDURE, PASS:: randPos => randPosVolTria - PROCEDURE, PASS:: area => areaTria - PROCEDURE, NOPASS:: fPsi => fPsiTria - PROCEDURE, NOPASS:: dPsi => dPsiTria - PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 - PROCEDURE, PASS:: partialDer => partialDerTria - PROCEDURE, PASS:: elemK => elemKTria - PROCEDURE, PASS:: elemF => elemFTria - PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: gatherEF => gatherEFTria - PROCEDURE, PASS:: gatherMF => gatherMFTria - PROCEDURE, PASS:: getNodes => getNodesTria - PROCEDURE, PASS:: phy2log => phy2logTria - PROCEDURE, PASS:: nextElement => nextElementTria + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellTria2DCart + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: randPos => randPosCellTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, NOPASS:: detJac => detJ2DCart + PROCEDURE, NOPASS:: invJac => invJ2DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFTria + PROCEDURE, PASS:: gatherMagneticField => gatherMFTria + PROCEDURE, PASS:: elemK => elemKTria + PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, NOPASS:: inside => insideTria + PROCEDURE, PASS:: phy2log => phy2logTria + PROCEDURE, PASS:: neighbourElement => neighbourElementTria + !PARTICULAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria - END TYPE meshVol2DCartTria + END TYPE meshCell2DCartTria CONTAINS !NODE FUNCTIONS @@ -160,7 +139,7 @@ MODULE moduleMesh2DCart END FUNCTION getCoord2DCart !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -176,6 +155,7 @@ MODULE moduleMesh2DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj !Get element coordinates @@ -183,6 +163,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) + self%weight = 1.D0 !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -203,40 +184,19 @@ MODULE moduleMesh2DCart END SUBROUTINE initEdge2DCart - !Random position in quadrilateral volume - FUNCTION randPosVolQuad(self) RESULT(r) - USE moduleRandom - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8):: r(1:3) - REAL(8):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - xii(1) = random(-1.D0, 1.D0) - xii(2) = random(-1.D0, 1.D0) - xii(3) = 0.D0 - - fPsi = self%fPsi(xii) - - r(1) = DOT_PRODUCT(fPsi, self%x) - r(2) = DOT_PRODUCT(fPsi, self%y) - r(3) = 0.D0 - - END FUNCTION randposVolQuad - !Get nodes from edge - PURE FUNCTION getNodes2DCart(self) RESULT(n) + PURE FUNCTION getNodes2DCart(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) END FUNCTION getNodes2DCart + !Calculate intersection between position and edge PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r) IMPLICIT NONE @@ -255,7 +215,7 @@ MODULE moduleMesh2DCart END FUNCTION intersection2DCartEdge - !Calculates a random position in edge + !Calculate a random position in edge FUNCTION randPosEdge(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -276,18 +236,24 @@ MODULE moduleMesh2DCart !VOLUME FUNCTIONS !QUAD FUNCTIONS - !Inits quadrilateral element - SUBROUTINE initVolQuad2DCart(self, n, p, nodes) + !Init element + 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 + !Assign node index self%n = n + + !Assign number of nodes of cell + self%nNodes = SIZE(p) + + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -301,139 +267,179 @@ MODULE moduleMesh2DCart self%y = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - self%n4%v = self%n4%v + self%arNodes(4) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolQuad2DCart + END SUBROUTINE initCellQuad2DCart - !Computes element area - PURE SUBROUTINE areaQuad(self) + !Get nodes from quadrilateral element + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(inout):: self - REAL(8):: xi(1:3) - REAL(8):: detJ + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/ self%n1%n, self%n2%n, self%n3%n, self%n4%n /) + + END FUNCTION getNodesQuad + + !Random position in quadrilateral volume + FUNCTION randPosCellQuad(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) - self%volume = 0.D0 - self%arNodes = 0.D0 - !2D 1 point Gauss Quad Integral - xi = 0.D0 - detJ = self%detJac(xi)*4.D0 !4 - fPsi = self%fPsi(xi) - self%volume = detJ - self%arNodes = fPsi*detJ + Xi(1) = random(-1.D0, 1.D0) + Xi(2) = random(-1.D0, 1.D0) + Xi(3) = 0.D0 - END SUBROUTINE areaQuad + fPsi = self%fPsi(Xi, 4) - !Computes element functions in point xi - PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) + r(1) = DOT_PRODUCT(fPsi, self%x) + r(2) = DOT_PRODUCT(fPsi, self%y) + r(3) = 0.D0 + + END FUNCTION randPosCellQuad + + !Compute element functions in point Xi + PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - REAL(8),INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - ALLOCATE(fPsi(1:4)) + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & + (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /) - fPsi(1) = (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 + fPsi = fPsi * 0.25D0 END FUNCTION fPsiQuad - !Derivative element function at coordinates xi - PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - ALLOCATE(dPsi(1:2,1:4)) + dPsi = 0.D0 - dPsi(1,:) = dPsiQuadXi1(xi(2)) - dPsi(2,:) = dPsiQuadXi2(xi(1)) + dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), & + (1.D0 - Xi(2)), & + (1.D0 + Xi(2)), & + -(1.D0 + Xi(2)) /) + + dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), & + -(1.D0 + Xi(1)), & + (1.D0 + Xi(1)), & + (1.D0 - Xi(1)) /) + + dPsi = dPsi * 0.25D0 END FUNCTION dPsiQuad - !Derivative element function (xi1) - PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) + !Partial derivative in global coordinates + PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - REAL(8),INTENT(in):: xi2 - REAL(8):: dPsiXi1(1:4) + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - 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 + pDer = 0.D0 - !Derivative element function (xi2) - PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /) + pDer(3,3) = 1.D0 + + END FUNCTION partialDerQuad + + 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, 4, 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, 4, B) + + END FUNCTION gatherMFQuad + + !Compute element local stiffness matrix + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) IMPLICIT NONE - 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 = dPsiXi2*0.25D0 - - END FUNCTION dPsiQuadXi2 - - PURE SUBROUTINE partialDerQuad(self, dPsi, dx, dy) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy - - dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) - dx(2) = DOT_PRODUCT(dPsi(2,:),self%x) - dy(1) = DOT_PRODUCT(dPsi(1,:),self%y) - dy(2) = DOT_PRODUCT(dPsi(2,:),self%y) - - END SUBROUTINE partialDerQuad - - !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(localK) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xi(1:3) - REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m - ALLOCATE(localK(1:4, 1:4)) - localK=0.D0 - xi=0.D0 + localK = 0.D0 + + Xi = 0.D0 !Start 2D Gauss Quad Integral - DO l=1, 3 - xi(2) = corQuad(l) - dPsi(1,:) = self%dPsiXi1(xi(2)) + DO l = 1, 3 + Xi(2) = corQuad(l) DO m = 1, 3 - xi(1) = corQuad(m) - dPsi(2,:) = self%dPsiXi2(xi(1)) - 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 + Xi(1) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + wQuad(l)*wQuad(m)/detJ END DO END DO @@ -441,27 +447,32 @@ MODULE moduleMesh2DCart END FUNCTION elemKQuad !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, source) RESULT(localF) + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: xi(1:3) + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: detJ, f INTEGER:: l, m - ALLOCATE(localF(1:4)) localF = 0.D0 - xi = 0.D0 - DO l=1, 3 - xi(1) = corQuad(l) + + Xi = 0.D0 + DO l = 1, 3 + Xi(1) = corQuad(l) DO m = 1, 3 - xi(2) = corQuad(m) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) - f = DOT_PRODUCT(fPsi,source) + Xi(2) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ END DO @@ -469,148 +480,113 @@ MODULE moduleMesh2DCart END FUNCTION elemFQuad - !Checks if a particle is inside a quad element - PURE FUNCTION insideQuad(xi) RESULT(ins) + !Check if Xi is inside the element + PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + 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) + 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) + !Transform physical coordinates to element coordinates + PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 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) /) - - fPsi = self%fPsi(xi) - 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(:) - - 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(xN) - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self + CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) - REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: Xi(1:3) + REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: dPsi(1:3,1:4), fPsi(1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv=1.D0 - xO=0.D0 + conv = 1.D0 + XiO = 0.D0 - DO WHILE(conv>1.D-4) - dPsi = self%dPsi(xO) - invJ = self%invJac(xO, dPsi) - fPsi = self%fPsi(xO) - f(1) = DOT_PRODUCT(fPsi,self%x)-r(1) - f(2) = DOT_PRODUCT(fPsi,self%y)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN + DO WHILE(conv > 1.D-2) + dPsi = self%dPsi(XiO, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(XiO, 4) + f = (/ DOT_PRODUCT(fPsi,self%x), & + DOT_PRODUCT(fPsi,self%y), & + 0.D0 /) + f = f - r + Xi = XiO - 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) + !Get the neighbour element for a logical position Xi + SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + CLASS(meshCell2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + REAL(8):: XiArray(1:4) INTEGER:: nextInt - 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) + XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) + nextInt = MAXLOC(XiArray,1) + !Select the higher value of directions and searches in that direction + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 CASE (4) - nextElement => self%e4 + neighbourElement => self%e4 END SELECT - END SUBROUTINE nextElementQuad + END SUBROUTINE neighbourElementQuad - !TRIA ELEMENT - !Init tria element - SUBROUTINE initVolTria2DCart(self, n, p, nodes) + !Compute element volume + PURE SUBROUTINE volumeQuad(self) + IMPLICIT NONE + + CLASS(meshCell2DCartQuad), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) + + self%volume = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + !Compute total volume of the cell + self%volume = detJ*4.D0 + !Compute volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + self%n4%v = self%n4%v + fPsi(4)*self%volume + + END SUBROUTINE volumeQuad + + !TRIA FUNCTIONS + !Init element + 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(:) @@ -619,6 +595,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 @@ -630,350 +609,308 @@ MODULE moduleMesh2DCart self%x = (/r1(1), r2(1), r3(1)/) self%y = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) 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) + !Random position in cell + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/self%n1%n, self%n2%n, self%n3%n /) + + END FUNCTION getNodesTria + + !Random position in cell + 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):: 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 - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 4) 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) + !Compute element functions in point Xi + PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(inout):: self - REAL(8):: xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - 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)/2.D0 - fPsi = self%fPsi(xi) - self%volume = detJ - self%arNodes = fPsi*detJ - - END SUBROUTINE areaTria - - !Shape functions for triangular element - PURE FUNCTION fPsiTria(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:3)) - - 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 FUNCTION fPsiTria - !Derivative element function at coordinates xi - PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + !Compute element derivative functions in point Xi + PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - ALLOCATE(dPsi(1:2,1:3)) + dPsi = 0.D0 - dPsi(1,:) = dPsiTriaXi1(xi(2)) - dPsi(2,:) = dPsiTriaXi2(xi(1)) + dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + !Compute the derivatives in global coordinates + PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - REAL(8), INTENT(in):: xi2 - REAL(8):: dPsiXi1(1:3) + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - dPsiXi1(1) = -1.D0 - dPsiXi1(2) = 1.D0 - dPsiXi1(3) = 0.D0 + pDer = 0.D0 - END FUNCTION dPsiTriaXi1 + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%x(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%x(1:3)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /) - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + END FUNCTION partialDerTria + + !Gather electric field at position Xi + 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, 3, phi) + + END FUNCTION gatherEFTria + + !Gather magnetic field at position Xi + 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 + + !Compute cell local stiffness matrix + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) IMPLICIT NONE - REAL(8), INTENT(in):: xi1 - REAL(8):: dPsiXi2(1:3) - - dPsiXi2(1) = -1.D0 - dPsiXi2(2) = 0.D0 - dPsiXi2(3) = 1.D0 - - END FUNCTION dPsiTriaXi2 - - PURE SUBROUTINE partialDerTria(self, dPsi, dx, dy) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy - - dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) - dx(2) = DOT_PRODUCT(dPsi(2,:),self%x) - dy(1) = DOT_PRODUCT(dPsi(1,:),self%y) - dy(2) = DOT_PRODUCT(dPsi(2,:),self%y) - - END SUBROUTINE partialDerTria - - !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(localK) - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xi(1:3) - REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: invJ(1:2,1:2), detJ + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3,1:3) + REAL(8):: pDer(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 + localK = 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) - fPsi = self%fPsi(xi) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO END FUNCTION elemKTria - !Computes element local source vector - PURE FUNCTION elemFTria(self, source) RESULT(localF) + !Compute element local source vector + PURE FUNCTION elemFTria(self, nNodes, 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:3) - REAL(8):: xi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) + REAL(8):: 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) - fPsi = self%fPsi(xi) - f = DOT_PRODUCT(fPsi,source) + DO l = 1, 4 + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) + f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wTria(l)*detJ END DO END FUNCTION elemFTria - PURE FUNCTION insideTria(xi) RESULT(ins) + !Check if Xi is inside the element + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + 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 + 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) + !Transform physical coordinates to element coordinates + PURE FUNCTION phy2logTria(self,r) RESULT(Xi) 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) /) - - fPsi = self%fPsi(xi) - 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(:) - - 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) - 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):: Xi(1:3) + REAL(8):: deltaR(1:3) + REAL(8):: dPsi(1:3, 1:3) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ !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, 3) + pDer = self%partialDer(3, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) + Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria - SUBROUTINE nextElementTria(self, xi, nextElement) + !Get the neighbour cell for a logical position Xi + SUBROUTINE neighbourElementTria(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:3) + CLASS(meshCell2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + REAL(8):: XiArray(1:3) INTEGER:: nextInt - xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) - nextInt = MINLOC(xiArray,1) - NULLIFY(nextElement) + XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray,1) + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 END SELECT - END SUBROUTINE nextElementTria + END SUBROUTINE neighbourElementTria - !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D - !Computes element Jacobian determinant - PURE FUNCTION detJ2DCart(self, xi, dPsi_in) RESULT(dJ) + !Calculate volume for triangular element + PURE SUBROUTINE volumeTria(self) IMPLICIT NONE - CLASS(meshVol2DCart), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) + CLASS(meshCell2DCartTria), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:3) + + self%volume = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) + !Computes total volume of the cell + self%volume = detJ + !Computes volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + + END SUBROUTINE volumeTria + + !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D + !Compute element Jacobian determinant + PURE FUNCTION detJ2DCart(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dx(1:2), dy(1:2) - IF(PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx, dy) - dJ = dx(1)*dy(2)-dx(2)*dy(1) + dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1) END FUNCTION detJ2DCart - !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCart(self,xi,dPsi_in) RESULT(invJ) + !Compute element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ2DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol2DCart), 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):: dx(1:2), dy(1:2) - REAL(8):: invJ(1:2,1:2) + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in + invJ = 0.D0 - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx, dy) - invJ(1,:) = (/ dy(2), -dx(2) /) - invJ(2,:) = (/ -dy(1), dx(1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(3, 3) = 1.D0 END FUNCTION invJ2DCart @@ -983,11 +920,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 @@ -995,9 +932,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 @@ -1008,34 +945,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) @@ -1043,22 +980,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) @@ -1066,13 +1003,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. & @@ -1115,8 +1052,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 @@ -1207,8 +1144,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 @@ -1280,7 +1217,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 @@ -1364,7 +1301,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 @@ -1381,7 +1318,7 @@ MODULE moduleMesh2DCart !Revers the normal to point inside the domain elemB%normal = - elemB%normal - + END IF END IF diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index c576c71..307f71c 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,15 +11,16 @@ 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 !Element coordinates REAL(8):: r = 0.D0, z = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode2DCyl + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode2DCyl PROCEDURE, PASS:: getCoordinates => getCoord2DCyl END TYPE meshNode2DCyl @@ -30,102 +31,79 @@ MODULE moduleMesh2DCyl !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdge2DCyl - PROCEDURE, PASS:: getNodes => getNodes2DCyl + !meshEdge DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initEdge2DCyl + PROCEDURE, PASS:: getNodes => getNodes2DCyl PROCEDURE, PASS:: intersection => intersection2DCylEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdge2DCyl - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl - CONTAINS - PROCEDURE, PASS:: detJac => detJ2DCyl - PROCEDURE, PASS:: invJac => invJ2DCyl - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol2DCyl - - 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:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - - END SUBROUTINE partialDer_interface - - END INTERFACE - !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylQuad + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylQuad !Element coordinates REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() - REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolQuad2DCyl - 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:: elemK => elemKQuad - PROCEDURE, PASS:: elemF => elemFQuad - PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: gatherEF => gatherEFQuad - PROCEDURE, PASS:: gatherMF => gatherMFQuad - PROCEDURE, PASS:: getNodes => getNodesQuad - PROCEDURE, PASS:: phy2log => phy2logQuad - PROCEDURE, PASS:: nextElement => nextElementQuad + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellQuad2DCyl + PROCEDURE, PASS:: getNodes => getNodesQuad + PROCEDURE, PASS:: randPos => randPosCellQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, NOPASS:: detJac => detJ2DCyl + PROCEDURE, NOPASS:: invJac => invJ2DCyl + PROCEDURE, PASS:: gatherElectricField => gatherEFQuad + PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad + PROCEDURE, PASS:: elemK => elemKQuad + PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: neighbourElement => neighbourElementQuad + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad - END TYPE meshVol2DCylQuad + END TYPE meshCell2DCylQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylTria + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylTria !Element coordinates REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() - REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolTria2DCyl - PROCEDURE, PASS:: randPos => randPosVolTria - PROCEDURE, PASS:: area => areaTria - PROCEDURE, NOPASS:: fPsi => fPsiTria - PROCEDURE, NOPASS:: dPsi => dPsiTria - PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 - PROCEDURE, PASS:: partialDer => partialDerTria - PROCEDURE, PASS:: elemK => elemKTria - PROCEDURE, PASS:: elemF => elemFTria - PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: gatherEF => gatherEFTria - PROCEDURE, PASS:: gatherMF => gatherMFTria - PROCEDURE, PASS:: getNodes => getNodesTria - PROCEDURE, PASS:: phy2log => phy2logTria - PROCEDURE, PASS:: nextElement => nextElementTria + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellTria2DCyl + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: randPos => randPosCellTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, NOPASS:: detJac => detJ2DCyl + PROCEDURE, NOPASS:: invJac => invJ2DCyl + PROCEDURE, PASS:: gatherElectricField => gatherEFTria + PROCEDURE, PASS:: gatherMagneticField => gatherMFTria + PROCEDURE, PASS:: elemK => elemKTria + PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, NOPASS:: inside => insideTria + PROCEDURE, PASS:: phy2log => phy2logTria + PROCEDURE, PASS:: neighbourElement => neighbourElementTria + !PARTICULAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria - END TYPE meshVol2DCylTria + END TYPE meshCell2DCylTria CONTAINS !NODE FUNCTIONS - !Inits node element + !Init node element SUBROUTINE initNode2DCyl(self, n, r) USE moduleSpecies USE moduleRefParam @@ -161,7 +139,7 @@ MODULE moduleMesh2DCyl END FUNCTION getCoord2DCyl !EDGE FUNCTIONS - !Inits edge element + !Init edge element SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -177,6 +155,7 @@ MODULE moduleMesh2DCyl INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj !Get element coordinates @@ -190,6 +169,7 @@ MODULE moduleMesh2DCyl self%z(2)-self%z(1) , & 0.D0 /) self%normal = self%normal/NORM2(self%normal) + !Boundary index self%boundary => boundary(bt) ALLOCATE(self%fboundary(1:nSpecies)) @@ -205,19 +185,19 @@ MODULE moduleMesh2DCyl END SUBROUTINE initEdge2DCyl !Get nodes from edge - PURE FUNCTION getNodes2DCyl(self) RESULT(n) + PURE FUNCTION getNodes2DCyl(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) END FUNCTION getNodes2DCyl + !Calculate intersection between position and edge PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r) - USE moduleMath IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self @@ -235,15 +215,15 @@ MODULE moduleMesh2DCyl END FUNCTION intersection2DCylEdge - !Calculates a random position in edge + !Calculate a random position in edge FUNCTION randPosEdge(self) RESULT(r) USE moduleRandom IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd - REAL(8):: dr, dz REAL(8):: r(1:3) + REAL(8):: dr, dz rnd = random() dr = self%r(2) - self%r(1) @@ -264,18 +244,24 @@ MODULE moduleMesh2DCyl !VOLUME FUNCTIONS !QUAD FUNCTIONS - !Inits quadrilateral element - SUBROUTINE initVolQuad2DCyl(self, n, p, nodes) + !Init element + 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 + !Assign node index self%n = n + + !Assign number of nodes of cell + self%nNodes = SIZE(p) + + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -289,181 +275,182 @@ MODULE moduleMesh2DCyl self%r = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - self%n4%v = self%n4%v + self%arNodes(4) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) 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 + !Get nodes from quadrilateral element + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(inout):: self - REAL(8):: r, xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:4) + CLASS(meshCell2DCylQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - 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 - 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/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) - self%arNodes(1) = fPsi(1)*r*detJ - xi = (/ 5.D-1, -5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) - self%arNodes(2) = fPsi(2)*r*detJ - xi = (/ 5.D-1, 5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) - self%arNodes(3) = fPsi(3)*r*detJ - xi = (/-5.D-1, 5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) - self%arNodes(4) = fPsi(4)*r*detJ + n = (/ self%n1%n, self%n2%n, self%n3%n, self%n4%n /) - END SUBROUTINE areaQuad - - !Computes element functions in point xi - PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8),INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:4)) - - 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 FUNCTION fPsiQuad - - !Derivative element function at coordinates xi - PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - ALLOCATE(dPsi(1:2,1:4)) - - dPsi(1,:) = dPsiQuadXi1(xi(2)) - dPsi(2,:) = dPsiQuadXi2(xi(1)) - - END FUNCTION dPsiQuad - - !Derivative element function (xi1) - PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) - IMPLICIT NONE - - 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 = dPsiXi1*0.25D0 - - END FUNCTION dPsiQuadXi1 - - !Derivative element function (xi2) - PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) - IMPLICIT NONE - - 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 = dPsiXi2*0.25D0 - - END FUNCTION dPsiQuadXi2 - - !Partial derivative in global coordinates - PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - - dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) - - END SUBROUTINE partialDerQuad + END FUNCTION getNodesQuad !Random position in quadrilateral volume - FUNCTION randPosVolQuad(self) RESULT(r) + FUNCTION randPosCellQuad(self) RESULT(r) 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 - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 4) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) r(3) = 0.D0 - END FUNCTION randposVolQuad + END FUNCTION randPosCellQuad + + !Computes element functions in point Xi + PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) + + fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 - Xi(2)), & + (1.D0 + Xi(1)) * (1.D0 + Xi(2)), & + (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /) + + fPsi = fPsi * 0.25D0 + + END FUNCTION fPsiQuad + + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) + + dPsi = 0.D0 + + dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), & + (1.D0 - Xi(2)), & + (1.D0 + Xi(2)), & + -(1.D0 + Xi(2)) /) + + dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), & + -(1.D0 + Xi(1)), & + (1.D0 + Xi(1)), & + (1.D0 - Xi(1)) /) + + dPsi = dPsi * 0.25D0 + + END FUNCTION dPsiQuad + + !Partial derivative in global coordinates + PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer) + IMPLICIT NONE + + CLASS(meshCell2DCylQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) + + pDer = 0.D0 + + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /) + pDer(3,3) = 1.D0 + + END FUNCTION partialDerQuad + + PURE FUNCTION gatherEFQuad(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):: phi(1:4) + + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) + + array = -self%gatherDF(Xi, 4, phi) + + END FUNCTION gatherEFQuad + + 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) + + 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, 4, B) + + END FUNCTION gatherMFQuad !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(localK) + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) 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 + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: r + REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l, m - ALLOCATE(localK(1:4, 1:4)) - localK=0.D0 - xi=0.D0 + localK = 0.D0 + + Xi = 0.D0 !Start 2D Gauss Quad Integral - DO l=1, 3 - xi(2) = corQuad(l) - dPsi(1,:) = self%dPsiXi1(xi(2)) + DO l = 1, 3 + Xi(2) = corQuad(l) DO m = 1, 3 - 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))* & - r*wQuad(l)*wQuad(m)/detJ + Xi(1) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + r = self%gatherF(Xi, 4, self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + r*wQuad(l)*wQuad(m)/detJ END DO END DO @@ -472,29 +459,35 @@ MODULE moduleMesh2DCyl END FUNCTION elemKQuad !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, source) RESULT(localF) + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: r REAL(8):: detJ, f INTEGER:: l, m - ALLOCATE(localF(1:4)) localF = 0.D0 - xi = 0.D0 - DO l=1, 3 - xi(1) = corQuad(l) + + Xi = 0.D0 + DO l = 1, 3 + Xi(1) = corQuad(l) DO m = 1, 3 - xi(2) = corQuad(m) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) + Xi(2) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ END DO @@ -503,148 +496,124 @@ MODULE moduleMesh2DCyl END FUNCTION elemFQuad - !Checks if a particle is inside a quad element - PURE FUNCTION insideQuad(xi) RESULT(ins) + !Checks if Xi is inside the element + PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + 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) + 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) + !Transform physical coordinates to element coordinates + PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 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 - 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(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) - - 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) /) - - fPsi = self%fPsi(xi) - MF = MATMUL(fPsi(:), MF_Nodes) - - END FUNCTION gatherMFQuad - - !Gets nodes from quadrilateral element - PURE FUNCTION getNodesQuad(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - 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(xN) - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(in):: self + CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) - REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: Xi(1:3) + REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: dPsi(1:3,1:4), fPsi(1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv=1.D0 - xO=0.D0 + conv = 1.D0 + XiO = 0.D0 - DO WHILE(conv>1.D-4) - dPsi = self%dPsi(xO) - invJ = self%invJac(xO, dPsi) - fPsi = self%fPsi(xO) - f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) - f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN + DO WHILE(conv > 1.D-2) + dPsi = self%dPsi(XiO, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(XiO, 4) + f = (/ DOT_PRODUCT(fPsi,self%z), & + DOT_PRODUCT(fPsi,self%r), & + 0.D0 /) + f = f - r + Xi = XiO - 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) + !Get the neighbour element for a logical position Xi + SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + CLASS(meshCell2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + 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) + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 CASE (4) - nextElement => self%e4 + neighbourElement => self%e4 END SELECT - END SUBROUTINE nextElementQuad + END SUBROUTINE neighbourElementQuad - !TRIA ELEMENT - !Init tria element - SUBROUTINE initVolTria2DCyl(self, n, p, nodes) + !Compute element volume + PURE SUBROUTINE volumeQuad(self) + USE moduleConstParam, ONLY: PI8 + IMPLICIT NONE + + CLASS(meshCell2DCylQuad), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: r + REAL(8):: detJ + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) + + self%volume = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + r = DOT_PRODUCT(fPsi,self%r) + !Computes total volume of the cell + self%volume = r*detJ*PI8 !4*2*pi + !Computes volume per node + Xi = (/-5.D-1, -5.D-1, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 + Xi = (/ 5.D-1, -5.D-1, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 + Xi = (/ 5.D-1, 5.D-1, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 + Xi = (/-5.D-1, 5.D-1, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 + + END SUBROUTINE volumeQuad + + !TRIA FUNCTIONS + !Init element + 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(:) @@ -653,6 +622,9 @@ MODULE moduleMesh2DCyl !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 @@ -664,157 +636,163 @@ MODULE moduleMesh2DCyl self%z = (/r1(1), r2(1), r3(1)/) self%r = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) 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) + !Random position in cell + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell2DCylTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/self%n1%n, self%n2%n, self%n3%n /) + + END FUNCTION getNodesTria + + !Random position in cell + FUNCTION randPosCellTria(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 - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 4) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) r(3) = 0.D0 - END FUNCTION randposVolTria + END FUNCTION randPosCellTria - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) - USE moduleConstParam, ONLY: PI + !Compute element functions in point Xi + PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(inout):: self - REAL(8):: r, xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - 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 - fPsi = self%fPsi(xi) - !Computes total volume of the cell - r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ - !Computes volume per node - self%arNodes = fPsi*r*detJ - - END SUBROUTINE areaTria - - !Shape functions for triangular element - PURE FUNCTION fPsiTria(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:3)) - - 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 FUNCTION fPsiTria - !Derivative element function at coordinates xi - PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + !Compute element derivative functions in point Xi + PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) - ALLOCATE(dPsi(1:2,1:3)) + dPsi = 0.D0 - dPsi(1,:) = dPsiTriaXi1(xi(2)) - dPsi(2,:) = dPsiTriaXi2(xi(1)) + dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + !Compute the derivatives in global coordinates + PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - REAL(8), INTENT(in):: xi2 - REAL(8):: dPsiXi1(1:3) + CLASS(meshCell2DCylTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - dPsiXi1(1) = -1.D0 - dPsiXi1(2) = 1.D0 - dPsiXi1(3) = 0.D0 + pDer = 0.D0 - END FUNCTION dPsiTriaXi1 + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%z(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%z(1:3)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%r(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%r(1:3)) /) - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + END FUNCTION partialDerTria + + !Gather electric field at position Xi + PURE FUNCTION gatherEFTria(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):: phi(1:3) - REAL(8), INTENT(in):: xi1 - REAL(8):: dPsiXi2(1:3) + phi = (/ self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi /) - dPsiXi2(1) = -1.D0 - dPsiXi2(2) = 0.D0 - dPsiXi2(3) = 1.D0 + array = -self%gatherDF(Xi, 3, phi) - END FUNCTION dPsiTriaXi2 + END FUNCTION gatherEFTria - PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) + !Gather magnetic field at position Xi + 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):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + B(:,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1) /) - dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) + B(:,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2) /) - END SUBROUTINE partialDerTria + B(:,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3) /) - !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(localK) + array = self%gatherF(Xi, 3, B) + + END FUNCTION gatherMFTria + + !Compute cell local stiffness matrix + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) 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 + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3,1:3) + REAL(8):: r + REAL(8):: pDer(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 + localK = 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) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + r = self%gatherF(Xi, 3, self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ END DO @@ -822,30 +800,35 @@ MODULE moduleMesh2DCyl END FUNCTION elemKTria - !Computes element local source vector - PURE FUNCTION elemFTria(self, source) RESULT(localF) + !Compute element local source vector + PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:3) - REAL(8):: r, xi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ, f + REAL(8):: r 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) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) + DO l = 1, 4 + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) localF = localF + r*f*fPsi*wTria(l)*detJ END DO @@ -853,171 +836,119 @@ MODULE moduleMesh2DCyl END FUNCTION elemFTria - PURE FUNCTION insideTria(xi) RESULT(ins) + !Check if Xi is inside the element + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + 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 + 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) + !Transform physical coordinates to element coordinates + PURE FUNCTION phy2logTria(self,r) RESULT(Xi) 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 - 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(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) - - 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) /) - - fPsi = self%fPsi(xi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFTria - - !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshVol2DCylTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - 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) - 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):: invJ(1:2,1:2), detJ - REAL(8):: deltaR(1:2) - REAL(8):: dPsi(1:2,1:3) + REAL(8):: Xi(1:3) + REAL(8):: deltaR(1:3) + REAL(8):: dPsi(1:3, 1:3) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - 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 + Xi = 0.D0 + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) + Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria - SUBROUTINE nextElementTria(self, xi, nextElement) + !Get the neighbour cell for a logical position Xi + SUBROUTINE neighbourElementTria(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:3) + CLASS(meshCell2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + REAL(8):: XiArray(1:3) INTEGER:: nextInt - xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) - nextInt = MINLOC(xiArray,1) - NULLIFY(nextElement) + XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray,1) + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 END SELECT - END SUBROUTINE nextElementTria + END SUBROUTINE neighbourElementTria - !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS - !Computes element Jacobian determinant - PURE FUNCTION detJ2DCyl(self, xi, dPsi_in) RESULT(dJ) + !Calculate volume for triangular element + PURE SUBROUTINE volumeTria(self) + USE moduleConstParam, ONLY: PI 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(meshCell2DCylTria), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:3) + REAL(8):: r + + self%volume = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) + r = DOT_PRODUCT(fPsi, self%r) + !Computes total volume of the cell + self%volume = r*detJ*PI !2PI*1/2 + !Computes volume per node + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + + END SUBROUTINE volumeTria + + !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS + !Compute element Jacobian determinant + PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dz(1:2), dr(1:2) - IF(PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dz, dr) - dJ = dz(1)*dr(2)-dz(2)*dr(1) + dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1) END FUNCTION detJ2DCyl - !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCyl(self,xi,dPsi_in) RESULT(invJ) + !Compute element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ2DCyl(pDer) 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(:,:) - REAL(8):: dz(1:2), dr(1:2) - REAL(8):: invJ(1:2,1:2) + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in + invJ = 0.D0 - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dz, dr) - invJ(1,:) = (/ dr(2), -dz(2) /) - invJ(2,:) = (/ -dr(1), dz(1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl @@ -1027,11 +958,11 @@ MODULE moduleMesh2DCyl 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 @@ -1039,9 +970,9 @@ MODULE moduleMesh2DCyl 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 @@ -1052,34 +983,34 @@ MODULE moduleMesh2DCyl END SUBROUTINE connectMesh2DCyl !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(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) @@ -1087,22 +1018,22 @@ MODULE moduleMesh2DCyl 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(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) @@ -1110,13 +1041,13 @@ MODULE moduleMesh2DCyl END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge 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. & @@ -1159,8 +1090,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 @@ -1251,8 +1182,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 @@ -1324,7 +1255,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 @@ -1408,7 +1339,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 0add18a..c451689 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -11,6 +11,7 @@ MODULE moduleMesh3DCart !Element coordinates REAL(8):: x, y, z CONTAINS + !meshNode DEFERRED PROCEDURES PROCEDURE, PASS:: init => initNode3DCart PROCEDURE, PASS:: getCoordinates => getCoord3DCart @@ -23,42 +24,18 @@ MODULE moduleMesh3DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge3DCartTria PROCEDURE, PASS:: getNodes => getNodes3DCartTria PROCEDURE, PASS:: intersection => intersection3DCartTria PROCEDURE, PASS:: randPos => randPosEdgeTria - PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria + !PARTICULAR PROCEDURES + PROCEDURE, NOPASS, PRIVATE:: fPsi => fPsiEdgeTria END TYPE meshEdge3DCartTria - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol3DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ3DCart - PROCEDURE, PASS:: invJac => invJ3DCart - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol3DCart - - ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xii) RESULT(dPsi) - REAL(8), INTENT(in):: xii(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:) - REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz - - END SUBROUTINE partialDer_interface - - END INTERFACE - !Tetrahedron volume element - TYPE, PUBLIC, EXTENDS(meshVol3DCart):: meshVol3DCartTetra + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell3DCartTetra !Element Coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -66,28 +43,30 @@ 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 => dPsiTetraXii1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXii2 - PROCEDURE, PASS:: partialDer => partialDerTetra - PROCEDURE, PASS:: elemK => elemKTetra - PROCEDURE, PASS:: elemF => elemFTetra - PROCEDURE, NOPASS:: inside => insideTetra - PROCEDURE, PASS:: gatherEF => gatherEFTetra - PROCEDURE, PASS:: gatherMF => gatherMFTetra - PROCEDURE, PASS:: getNodes => getNodesTetra - PROCEDURE, PASS:: phy2log => phy2logTetra - PROCEDURE, PASS:: nextElement => nextElementTetra + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellTetra + PROCEDURE, PASS:: getNodes => getNodesTetra + PROCEDURE, PASS:: randPos => randPosCellTetra + PROCEDURE, NOPASS:: fPsi => fPsiTetra + PROCEDURE, NOPASS:: dPsi => dPsiTetra + PROCEDURE, PASS:: partialDer => partialDerTetra + PROCEDURE, NOPASS:: detJac => detJ3DCart + PROCEDURE, NOPASS:: invJac => invJ3DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFTetra + PROCEDURE, PASS:: gatherMagneticField => gatherMFTetra + PROCEDURE, PASS:: elemK => elemKTetra + PROCEDURE, PASS:: elemF => elemFTetra + PROCEDURE, NOPASS:: inside => insideTetra + PROCEDURE, PASS:: phy2log => phy2logTetra + PROCEDURE, PASS:: neighbourElement => neighbourElementTetra + !PARTICULAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTetra - END TYPE meshVol3DCartTetra + END TYPE meshCell3DCartTetra CONTAINS !NODE FUNCTIONS - !Inits node element + !Init node element SUBROUTINE initNode3DCart(self, n, r) USE moduleSpecies USE moduleRefParam @@ -123,8 +102,8 @@ MODULE moduleMesh3DCart END FUNCTION getCoord3DCart - !SURFACE FUNCTIONS - !Inits surface element + !EDGE FUNCTIONS + !Init surface element SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary @@ -142,6 +121,7 @@ MODULE moduleMesh3DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj self%n3 => mesh%nodes(p(3))%obj @@ -177,17 +157,18 @@ MODULE moduleMesh3DCart END SUBROUTINE initEdge3DCartTria !Get nodes from surface - PURE FUNCTION getNodes3DCartTria(self) RESULT(n) + PURE FUNCTION getNodes3DCartTria(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:3)) n = (/self%n1%n, self%n2%n, self%n3%n/) END FUNCTION getNodes3DCartTria + !Calculate intersection between position and edge PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r) IMPLICIT NONE @@ -206,21 +187,21 @@ MODULE moduleMesh3DCart END FUNCTION intersection3DCartTria - !Calculates a random position in the surface + !Calculate a random position in the surface FUNCTION randPosEdgeTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) + 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 - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi) r = (/DOT_PRODUCT(fPsi, self%x), & DOT_PRODUCT(fPsi, self%y), & DOT_PRODUCT(fPsi, self%z)/) @@ -228,35 +209,38 @@ MODULE moduleMesh3DCart END FUNCTION randPosEdgeTria !Shape functions for triangular surface - PURE FUNCTION fPsiEdgeTria(xii) RESULT(fPsi) + PURE FUNCTION fPsiEdgeTria(Xi) RESULT(fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + REAL(8):: fPsi(1:3) - ALLOCATE(fPsi(1:3)) - - fPsi(1) = 1.D0 - xii(1) - xii(2) - fPsi(2) = xii(1) - fPsi(3) = xii(2) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) END FUNCTION fPsiEdgeTria !VOLUME FUNCTIONS !TETRA FUNCTIONS - !Inits tetrahedron element - SUBROUTINE initVolTetra3DCart(self, n, p, nodes) + !Init element + SUBROUTINE initCellTetra(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):: volNodes(1:4) !Volume of each node + !Assign node index self%n = n + + !Assign number of nodes of cell + self%nNodes = SIZE(p) + + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -271,431 +255,338 @@ MODULE moduleMesh3DCart self%z = (/r1(3), r2(3), r3(3), r4(3)/) !Computes the element volume - CALL self%calcVol() - - !Assign proportional volume to each node - !TODO: Review this to apply to all elements in the future - volNodes = self%fPsi((/0.25D0, 0.25D0, 0.25D0/))*self%volume - self%n1%v = self%n1%v + volNodes(1) - self%n2%v = self%n2%v + volNodes(2) - self%n3%v = self%n3%v + volNodes(3) - self%n4%v = self%n4%v + volNodes(4) + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initVolTetra3DCart + END SUBROUTINE initCellTetra - !Random position in volume tetrahedron - FUNCTION randPosVolTetra(self) RESULT(r) + !Gets node indexes from cell + PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell3DCartTetra), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) + + END FUNCTION getNodesTetra + + !Random position in cell + 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):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:4) - xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0 - xii(1)) - xii(3) = random( 0.D0, 1.D0 - xii(1) - xii(2)) + Xi(1) = random( 0.D0, 1.D0) + Xi(2) = random( 0.D0, 1.D0 - Xi(1)) + Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2)) - ALLOCATE(fPsi(1:4)) - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi, 4) - r(1) = DOT_PRODUCT(fPsi, self%x) - r(2) = DOT_PRODUCT(fPsi, self%y) - r(3) = DOT_PRODUCT(fPsi, self%z) + 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) + !Compute element functions in point Xi + PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(inout):: self - REAL(8):: xii(1:3) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - self%volume = 0.D0 - xii = (/0.25D0, 0.25D0, 0.25D0/) - self%volume = self%detJac(xii) - - END SUBROUTINE volumeTetra - - !Computes element functions in point xii - PURE FUNCTION fPsiTetra(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:4)) - - fPsi(1) = 1.D0 - xi(1) - xi(2) - xi(3) - fPsi(2) = xi(1) - fPsi(3) = xi(2) - fPsi(4) = xi(3) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) + fPsi(4) = Xi(3) END FUNCTION fPsiTetra - !Derivative element function at coordinates xii - PURE FUNCTION dPsiTetra(xii) RESULT(dPsi) + !Compute element derivative functions in point Xi + PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xii(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3, 1:nNodes) - ALLOCATE(dPsi(1:3,1:4)) + dPsi = 0.D0 - dPsi(1,:) = dPsiTetraXii1(xii(2), xii(3)) - dPsi(2,:) = dPsiTetraXii2(xii(1), xii(3)) - dPsi(3,:) = dPsiTetraXii3(xii(1), xii(2)) + dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /) + dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /) + dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTetra - !Derivative element function respect to xii1 - PURE FUNCTION dPsiTetraXii1(xii2, xii3) RESULT(dPsiXii1) - IMPLICIT NONE - REAL(8), INTENT(in):: xii2, xii3 - REAL(8):: dPsiXii1(1:4) - - dPsiXii1(1) = -1.D0 - dPsiXii1(2) = 1.D0 - dPsiXii1(3) = 0.D0 - dPsiXii1(4) = 0.D0 - - END FUNCTION dPsiTetraXii1 - - !Derivative element function respect to xii2 - PURE FUNCTION dPsiTetraXii2(xii1, xii3) RESULT(dPsiXii2) - IMPLICIT NONE - REAL(8), INTENT(in):: xii1, xii3 - REAL(8):: dPsiXii2(1:4) - - dPsiXii2(1) = -1.D0 - dPsiXii2(2) = 0.D0 - dPsiXii2(3) = 1.D0 - dPsiXii2(4) = 0.D0 - - END FUNCTION dPsiTetraXii2 - - !Derivative element function respect to xii3 - PURE FUNCTION dPsiTetraXii3(xii1, xii2) RESULT(dPsiXii3) - IMPLICIT NONE - REAL(8), INTENT(in):: xii1, xii2 - REAL(8):: dPsiXii3(1:4) - - dPsiXii3(1) = -1.D0 - dPsiXii3(2) = 0.D0 - dPsiXii3(3) = 0.D0 - dPsiXii3(4) = 1.D0 - - END FUNCTION dPsiTetraXii3 - - !Computes the derivatives in global coordinates - PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz) + !Compute the derivatives in global coordinates + PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:, 1:) - REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz + CLASS(meshCell3DCartTetra), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3, 1:3) - dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) - dx(2) = DOT_PRODUCT(dPsi(2,:), self%x) - dx(3) = DOT_PRODUCT(dPsi(3,:), self%x) + pDer = 0.D0 - dy(1) = DOT_PRODUCT(dPsi(1,:), self%y) - dy(2) = DOT_PRODUCT(dPsi(2,:), self%y) - dy(3) = DOT_PRODUCT(dPsi(3,:), self%y) + pDer(1, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%x(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%x(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%x(1:4)) /) - dz(1) = DOT_PRODUCT(dPsi(1,:), self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:), self%z) - dz(3) = DOT_PRODUCT(dPsi(3,:), self%z) + pDer(2, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%y(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%y(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%y(1:4)) /) - END SUBROUTINE partialDerTetra + pDer(3, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%z(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%z(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%z(1:4)) /) - PURE FUNCTION elemKTetra(self) RESULT(localK) + END FUNCTION partialDerTetra + + !Gather electric field at position Xi + 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, 4, phi) + + END FUNCTION gatherEFTetra + + !Gather magnetic field at position Xi + 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, 4, B) + + END FUNCTION gatherMFTetra + + !Compute cell local stiffness matrix + PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xii(1:3) + CLASS(meshCell3DCartTetra), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3), detJ - ALLOCATE(localK(1:4,1:4)) localK = 0.D0 - xii = 0.D0 + Xi = 0.D0 !TODO: One point Gauss integral. Upgrade when possible - xii = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(xii) - detJ = self%detJac(xii, dPsi) - invJ = self%invJac(xii, dPsi) - fPsi = self%fPsi(xii) + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(Xi, 4) localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra - PURE FUNCTION elemFTetra(self, source) RESULT(localF) + !Compute element local source vector + PURE FUNCTION elemFTetra(self, nNodes, 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 + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) - REAL(8):: xii(1:3) + REAL(8):: pDer(1:3, 1:3) REAL(8):: detJ, f - ALLOCATE(localF(1:4)) localF = 0.D0 - xii = 0.D0 - !TODO: One point Gauss integral. Upgrade when possible - xii = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(xii) - detJ = self%detJac(xii, dPsi) - fPsi = self%fPsi(xii) + Xi = 0.D0 + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi, source) localF = f*fPsi*1.D0*detJ END FUNCTION elemFTetra - PURE FUNCTION insideTetra(xi) RESULT(ins) + !Check if Xi is inside the element + PURE FUNCTION insideTetra(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >= 0.D0 .AND. & - xi(2) >= 0.D0 .AND. & - xi(3) >= 0.D0 .AND. & - 1.D0 - xi(1) - xi(2) - xi(3) >= 0.D0 + ins = Xi(1) >= 0.D0 .AND. & + Xi(2) >= 0.D0 .AND. & + Xi(3) >= 0.D0 .AND. & + 1.D0 - Xi(1) - Xi(2) - Xi(3) >= 0.D0 END FUNCTION insideTetra - PURE FUNCTION gatherEFTetra(self, xi) RESULT(EF) + !Transform physical coordinates to element coordinates + PURE FUNCTION phy2logTetra(self,r) RESULT(Xi) 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) /) - - fPsi = self%fPsi(xi) - MF = MATMUL(fPsi, MF_Nodes) - - END FUNCTION gatherMFTetra - - PURE FUNCTION getNodesTetra(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:4)) - n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) - - END FUNCTION getNodesTetra - - 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):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: invJ(1:3, 1:3), detJ REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:4) - xi = 0.D0 + !Direct method to convert coordinates + Xi = 0.D0 deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /) - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - detJ = self%detJac(xi, dPsi) - xi = MATMUL(invJ, deltaR)/detJ + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) + Xi = MATMUL(invJ, deltaR)/detJ END FUNCTION phy2logTetra - SUBROUTINE nextElementTetra(self, xi, nextElement) + !Get the neighbour cell for a logical position Xi + SUBROUTINE neighbourElementTetra(self, Xi, neighbourElement) IMPLICIT NONE - CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + CLASS(meshCell3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + REAL(8):: XiArray(1:4) INTEGER:: nextInt !TODO: Review when connectivity - xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /) - nextInt = MINLOC(xiArray, 1) - NULLIFY(nextElement) + XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray, 1) + NULLIFY(neighbourElement) SELECT CASE(nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 CASE (4) - nextElement => self%e4 + neighbourElement => self%e4 END SELECT - END SUBROUTINE nextElementTetra + END SUBROUTINE neighbourElementTetra - !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D - !Computes element Jacobian determinant - PURE FUNCTION detJ3DCart(self, xi, dPsi_in) RESULT(dJ) + !Calculate volume for triangular element + PURE SUBROUTINE volumeTetra(self) IMPLICIT NONE - CLASS(meshVol3DCart), INTENT(in)::self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:, 1:) + CLASS(meshCell3DCartTetra), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) + + self%volume = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = (/0.25D0, 0.25D0, 0.25D0/) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + !Computes total volume of the cell + self%volume = detJ + !Computes volume per node + fPsi = self%fPsi(Xi, 4) + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + self%n4%v = self%n4%v + fPsi(4)*self%volume + + END SUBROUTINE volumeTetra + + !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D + !Compute element Jacobian determinant + PURE FUNCTION detJ3DCart(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8), ALLOCATABLE:: dPsi(:,:) - REAL(8):: dx(1:3), dy(1:3), dz(1:3) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx, dy, dz) - dJ = dx(1)*(dy(2)*dz(3) - dy(3)*dz(2)) & - - dx(2)*(dy(1)*dz(3) - dy(3)*dz(1)) & - + dx(3)*(dy(1)*dz(2) - dy(2)*dz(1)) + dJ = pDer(1,1)*(pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)) & + - pDer(1,2)*(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)) & + + pDer(1,3)*(pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1)) END FUNCTION detJ3DCart - PURE FUNCTION invJ3DCart(self,xi,dPsi_in) RESULT(invJ) + !Compute element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ3DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshVol3DCart), 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), DIMENSION(1:3):: dx, dy, dz + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3) - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in + invJ(1,1:3) = (/ (pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)), & + -(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)), & + (pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1)) /) - ELSE - dPsi = self%dPsi(xi) + invJ(2,1:3) = (/ -(pDer(1,2)*pDer(3,3) - pDer(1,3)*pDer(3,2)), & + (pDer(1,1)*pDer(3,3) - pDer(1,3)*pDer(3,1)), & + -(pDer(1,1)*pDer(3,2) - pDer(1,2)*pDer(3,1)) /) - END IF - - CALL self%partialDer(dPsi, dx, dy, dz) - invJ(1,1) = (dy(2)*dz(3) - dy(3)*dz(2)) - invJ(1,2) = -(dy(1)*dz(3) - dy(3)*dz(1)) - invJ(1,3) = (dy(1)*dz(2) - dy(2)*dz(1)) - - invJ(2,1) = -(dx(2)*dz(3) - dx(3)*dz(2)) - invJ(2,2) = (dx(1)*dz(3) - dx(3)*dz(1)) - invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1)) - - invJ(3,1) = (dx(2)*dy(3) - dx(3)*dy(2)) - invJ(3,2) = -(dx(1)*dy(3) - dx(3)*dy(1)) - invJ(3,3) = (dx(1)*dy(2) - dx(2)*dy(1)) + invJ(3,1:3) = (/ (pDer(1,2)*pDer(2,3) - pDer(1,3)*pDer(2,2)), & + -(pDer(1,1)*pDer(2,3) - pDer(1,3)*pDer(2,1)), & + (pDer(1,1)*pDer(2,2) - pDer(1,2)*pDer(2,1)) /) invJ = TRANSPOSE(invJ) END FUNCTION invJ3DCart - !Selects type of elements to build connection - SUBROUTINE connectVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) - !Element A is a tetrahedron - SELECT TYPE(elemB) - TYPE IS(meshVol3DCartTetra) - !Element B is a tetrahedron - CALL connectTetraTetra(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectVolVol - - SUBROUTINE connectVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemB) - CLASS IS(meshEdge3DCartTria) - SELECT TYPE(elemA) - TYPE IS(meshVol3DCartTetra) - !Element A is a tetrahedron - CALL connectTetraEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectVolEdge - SUBROUTINE connectMesh3DCart(self) IMPLICIT NONE 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 @@ -703,9 +594,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 @@ -715,6 +606,46 @@ MODULE moduleMesh3DCart END SUBROUTINE connectMesh3DCart + !Select type of elements to build connection + SUBROUTINE connectCellCell(elemA, elemB) + IMPLICIT NONE + + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshCell), INTENT(inout):: elemB + + SELECT TYPE(elemA) + TYPE IS(meshCell3DCartTetra) + !Element A is a tetrahedron + SELECT TYPE(elemB) + TYPE IS(meshCell3DCartTetra) + !Element B is a tetrahedron + CALL connectTetraTetra(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectCellCell + + SUBROUTINE connectCellEdge(elemA, elemB) + IMPLICIT NONE + + CLASS(meshCell), INTENT(inout):: elemA + CLASS(meshEdge), INTENT(inout):: elemB + + SELECT TYPE(elemB) + CLASS IS(meshEdge3DCartTria) + SELECT TYPE(elemA) + TYPE IS(meshCell3DCartTetra) + !Element A is a tetrahedron + CALL connectTetraEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectCellEdge + !Checks if two sets of nodes are coincidend in any order PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) IMPLICIT NONE @@ -741,8 +672,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 @@ -870,11 +801,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 /) @@ -889,10 +820,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 @@ -922,10 +853,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 @@ -955,10 +886,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 @@ -988,10 +919,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..d968ba3 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,15 +59,14 @@ 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 - SUBROUTINE readInitial0D(sp, filename, density, velocity, temperature) + SUBROUTINE readInitial0D(filename, density, velocity, temperature) IMPLICIT NONE - INTEGER, INTENT(in):: sp CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity diff --git a/src/modules/mesh/inout/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..ae1fe05 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 @@ -321,10 +321,9 @@ MODULE moduleMeshInputGmsh2 END SUBROUTINE readGmsh2 !Reads the initial information from an output file for an species - SUBROUTINE readInitialGmsh2(sp, filename, density, velocity, temperature) + SUBROUTINE readInitialGmsh2(filename, density, velocity, temperature) IMPLICIT NONE - INTEGER, INTENT(in):: sp CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity diff --git a/src/modules/mesh/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 05c4b85..7482842 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -24,8 +24,10 @@ MODULE moduleMesh !Lock indicator for scattering INTEGER(KIND=OMP_LOCK_KIND):: lock CONTAINS + !DEFERED PROCEDURES PROCEDURE(initNode_interface), DEFERRED, PASS:: init PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates + !GENERIC PROCEDURES PROCEDURE, PASS:: resetOutput END TYPE meshNode @@ -66,10 +68,12 @@ 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() + !Nomber of nodes in the edge + INTEGER:: nNodes + !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 @@ -81,6 +85,7 @@ MODULE moduleMesh !Physical surface for the edge INTEGER:: physicalSurface CONTAINS + !DEFERED PROCEDURES PROCEDURE(initEdge_interface), DEFERRED, PASS:: init PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection @@ -102,10 +107,11 @@ MODULE moduleMesh END SUBROUTINE initEdge_interface !Get nodes index from node - PURE FUNCTION getNodesEdge_interface(self) RESULT(n) + PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) END FUNCTION getNodesEdge_interface @@ -146,8 +152,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 +169,160 @@ 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 - PROCEDURE, PASS:: findCell - PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log - PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside - PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement + !DEFERRED PROCEDURES + !Init the cell + PROCEDURE(initCell_interface), DEFERRED, PASS:: init + !Get the index of the nodes + PROCEDURE(getNodesCell_interface), DEFERRED, PASS:: getNodes + !Calculate random position on the cell + PROCEDURE(randPosCell_interface), DEFERRED, PASS:: randPos + !Obtain functions and values of cell natural functions + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi + PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + PROCEDURE(detJac_interface), DEFERRED, NOPASS:: detJac + PROCEDURE(invJac_interface), DEFERRED, NOPASS:: invJac + !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 + !Check if particle is inside the cell + PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside + !Convert physical coordinates (r) into logical coordinates (Xi) + PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log + !Returns the neighbour element based on particle position outside the cell + PROCEDURE(neighbourElement_interface), DEFERRED, PASS:: neighbourElement + !Scatter properties of particles on cell nodes + PROCEDURE, PASS:: scatter + !Subroutine to find in which cell a particle is located + PROCEDURE, PASS:: findCell + !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 - 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 + END SUBROUTINE initCell_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) + PURE FUNCTION getNodesCell_interface(self, nNodes) RESULT(n) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - END FUNCTION gatherEF_interface + END FUNCTION getNodesCell_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) + FUNCTION randPosCell_interface(self) RESULT(r) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8):: r(1:3) - END FUNCTION gatherMF_interface + END FUNCTION randPosCell_interface - PURE FUNCTION getNodesVol_interface(self) RESULT(n) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - END FUNCTION getNodesVol_interface - - PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) END FUNCTION fPsi_interface - PURE FUNCTION elemK_interface(self) RESULT(localK) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), ALLOCATABLE:: localK(:,:) + PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3, 1:nNodes) + + END FUNCTION dPsi_interface + + PURE FUNCTION partialDer_interface(self, nNodes, dPsi) RESULT(pDer) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3, 1:3) + + END FUNCTION partialDer_interface + + PURE FUNCTION detJac_interface(pDer) RESULT(dJ) + REAL(8), INTENT(in):: pDer(1:3,1:3) + REAL(8):: dJ + + END FUNCTION detJac_interface + + PURE FUNCTION invJac_interface(pDer) RESULT(invJ) + REAL(8), INTENT(in):: pDer(1:3,1:3) + 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, nNodes) RESULT(localK) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1: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(:) + PURE FUNCTION elemF_interface(self, nNodes, source) RESULT(localF) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1: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) - 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 - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(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 - REAL(8):: r(1:3) + PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) + IMPORT:: meshCell + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: Xi(1:3) - END FUNCTION randPosVol_interface + END FUNCTION phy2log_interface + + SUBROUTINE neighbourElement_interface(self, Xi, neighbourElement) + IMPORT:: meshCell, meshElement + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement + + END SUBROUTINE neighbourElement_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,16 +331,18 @@ 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(:) + !PROCEDURES SPECIFIC OF FILE TYPE PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() CONTAINS + !GENERIC PROCEDURES PROCEDURE, PASS:: doCollisions END TYPE meshGeneric @@ -295,14 +351,12 @@ MODULE moduleMesh !Reads the mesh from a file SUBROUTINE readMesh_interface(self, filename) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename END SUBROUTINE readMesh_interface - SUBROUTINE readInitial_interface(sp, filename, density, velocity, temperature) - INTEGER, INTENT(in):: sp + SUBROUTINE readInitial_interface(filename, density, velocity, temperature) CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity @@ -310,18 +364,16 @@ 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 - CLASS(meshGeneric), INTENT(inout):: self 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 - CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t @@ -338,28 +390,21 @@ MODULE moduleMesh REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV + !PROCEDURES SPECIFIC OF FILE TYPE PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS + !GENERIC PROCEDURES PROCEDURE, PASS:: constructGlobalK END TYPE meshParticles ABSTRACT INTERFACE - !Perform Coulomb Scattering - SUBROUTINE doCoulomb_interface(self) - IMPORT meshParticles - - CLASS(meshParticles), INTENT(inout):: self - - END SUBROUTINE doCoulomb_interface - !Prints Species data SUBROUTINE printOutput_interface(self, t) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t @@ -368,21 +413,25 @@ MODULE moduleMesh !Prints EM info SUBROUTINE printEM_interface(self, t) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface + !Perform Coulomb Scattering + SUBROUTINE doCoulomb_interface(self) + IMPORT meshParticles + CLASS(meshParticles), INTENT(inout):: self + + END SUBROUTINE doCoulomb_interface + !Prints average values SUBROUTINE printAverage_interface(self) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self END SUBROUTINE printAverage_interface - END INTERFACE TYPE(meshParticles), TARGET:: mesh @@ -390,6 +439,7 @@ MODULE moduleMesh !Collision (MCC) mesh TYPE, EXTENDS(meshGeneric):: meshCollisions CONTAINS + !GENERIC PROCEDURES END TYPE meshCollisions @@ -398,7 +448,6 @@ MODULE moduleMesh ABSTRACT INTERFACE SUBROUTINE readMeshColl_interface(self, filename) IMPORT meshCollisions - CLASS(meshCollisions), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename @@ -416,7 +465,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 @@ -431,24 +480,29 @@ MODULE moduleMesh !Logical to indicate if an specific mesh for MC Collisions is used LOGICAL:: doubleMesh + !Logical to indicate if MCC collisions are performed + LOGICAL:: doMCC !Complete path for the two meshes CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle CONTAINS !Constructs the global K matrix - SUBROUTINE constructGlobalK(self) + PURE SUBROUTINE constructGlobalK(self) IMPLICIT NONE CLASS(meshParticles), INTENT(inout):: self INTEGER:: e + INTEGER:: nNodes INTEGER, ALLOCATABLE:: n(:) REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: nNodes, i, j + INTEGER:: i, j - DO e = 1, self%numVols - n = self%vols(e)%obj%getNodes() - localK = self%vols(e)%obj%elemK() - nNodes = SIZE(n) + DO e = 1, self%numCells + nNodes = self%cells(e)%obj%nNodes + ALLOCATE(n(1:nNodes)) + ALLOCATE(localK(1:nNodes, 1:nNodes)) + n = self%cells(e)%obj%getNodes(nNodes) + localK = self%cells(e)%obj%elemK(nNodes) DO i = 1, nNodes DO j = 1, nNodes @@ -457,6 +511,8 @@ MODULE moduleMesh END DO END DO + + DEALLOCATE(n, localK) END DO @@ -480,30 +536,89 @@ MODULE moduleMesh END SUBROUTINE resetOutput - !Scatters particle properties into vol nodes - SUBROUTINE scatter(self, part) + !Gather the value of valNodes (scalar) at position Xi + PURE FUNCTION gatherF_scalar(self, Xi, nNodes, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) + REAL(8):: f + REAL(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = DOT_PRODUCT(fPsi, valNodes) + + END FUNCTION gatherF_scalar + + !Gather the value of valNodes (array) at position Xi + PURE FUNCTION gatherF_array(self, Xi, nNodes, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3) + REAL(8):: f(1:3) + REAL(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = MATMUL(fPsi, valNodes) + + END FUNCTION gatherF_array + + !Gather the spatial derivative of valNodes (scalar) at position Xi + PURE FUNCTION gatherDF_scalar(self, Xi, nNodes, valNodes) RESULT(df) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) + REAL(8):: df(1:3) + REAL(8):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3,1:3) + REAL(8):: dPsiR(1:3, 1:nNodes) + REAL(8):: invJ(1:3, 1:3), detJ + + dPsi = self%dPsi(Xi, nNodes) + pDer = self%partialDer(nNodes, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + 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, nNodes, part) USE moduleMath USE moduleSpecies USE OMP_LIB IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: self + CLASS(meshCell), INTENT(inout):: self + INTEGER, INTENT(in):: nNodes CLASS(particle), INTENT(in):: part - REAL(8), ALLOCATABLE:: fPsi(:) - INTEGER, ALLOCATABLE:: volNodes(:) + REAL(8):: fPsi(1:nNodes) + INTEGER:: cellNodes(1:nNodes) REAL(8):: tensorS(1:3, 1:3) INTEGER:: sp - INTEGER:: i, nNodes + INTEGER:: i CLASS(meshNode), POINTER:: node - fPsi = self%fPsi(part%xi) + cellNodes = self%getNodes(nNodes) + fPsi = self%fPsi(part%Xi, nNodes) + tensorS = outerProduct(part%v, part%v) + sp = part%species%n - volNodes = self%getNodes() - nNodes = SIZE(volNodes) DO i = 1, nNodes - node => mesh%nodes(volNodes(i))%obj + 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(:) @@ -521,38 +636,41 @@ 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(meshElement), POINTER:: nextElement + CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell + REAL(8):: Xi(1:3) + CLASS(meshElement), POINTER:: neighbourElement INTEGER:: sp - xi = self%phy2log(part%r) + Xi = self%phy2log(part%r) !Checks if particle is inside 'self' cell - IF (self%inside(xi)) THEN - part%vol = self%n - part%xi = xi + IF (self%inside(Xi)) THEN + part%cell = self%n + part%Xi = Xi part%n_in = .TRUE. !Assign particle to listPart_in - CALL OMP_SET_LOCK(self%lock) - sp = part%species%n - CALL self%listPart_in(sp)%add(part) - self%totalWeight(sp) = self%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(self%lock) + IF (doMCC) THEN + CALL OMP_SET_LOCK(self%lock) + sp = part%species%n + CALL self%listPart_in(sp)%add(part) + self%totalWeight(sp) = self%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(self%lock) + + END IF ELSE !If not, searches for a neighbour and repeats the process. - CALL self%nextElement(xi, nextElement) + CALL self%neighbourElement(Xi, neighbourElement) !Defines the next step - SELECT TYPE(nextElement) - CLASS IS(meshVol) + SELECT TYPE(neighbourElement) + CLASS IS(meshCell) !Particle moved to new cell, repeat find procedure - CALL nextElement%findCell(part, self) + CALL neighbourElement%findCell(part, self) CLASS IS (meshEdge) !Particle encountered a surface, apply boundary - CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part) + CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part) !If particle is still inside the domain, call findCell IF (part%n_in) THEN @@ -575,14 +693,14 @@ MODULE moduleMesh END SUBROUTINE findCell - !If Coll and Particle are the same, simply copy the part%vol into part%volColl + !If Coll and Particle are the same, simply copy the part%cell into part%cellColl SUBROUTINE findCellSameMesh(part) USE moduleSpecies IMPLICIT NONE TYPE(particle), INTENT(inout):: part - part%volColl = part%vol + part%cellColl = part%cell END SUBROUTINE findCellSameMesh @@ -595,31 +713,34 @@ MODULE moduleMesh TYPE(particle), INTENT(inout):: part LOGICAL:: found - CLASS(meshVol), POINTER:: vol - REAL(8), DIMENSION(1:3):: xii - CLASS(meshElement), POINTER:: nextElement + CLASS(meshCell), POINTER:: cell + REAL(8), DIMENSION(1:3):: Xi + CLASS(meshElement), POINTER:: neighbourElement INTEGER:: sp found = .FALSE. - vol => meshColl%vols(part%volColl)%obj + cell => meshColl%cells(part%cellColl)%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) - 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) + Xi = cell%phy2log(part%r) + IF (cell%inside(Xi)) THEN + part%cellColl = cell%n + IF (doMCC) THEN + CALL OMP_SET_LOCK(cell%lock) + sp = part%species%n + CALL cell%listPart_in(sp)%add(part) + cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(cell%lock) + + END IF found = .TRUE. ELSE - CALL vol%nextElement(xii, nextElement) - SELECT TYPE(nextElement) - CLASS IS(meshVol) + CALL cell%neighbourElement(Xi, neighbourElement) + SELECT TYPE(neighbourElement) + CLASS IS(meshCell) !Try next element - vol => nextElement + cell => neighbourElement CLASS DEFAULT !Should never happend, but just in case, stops loops @@ -644,15 +765,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 @@ -675,7 +796,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 @@ -686,21 +807,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 @@ -710,8 +832,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 @@ -721,15 +843,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 @@ -737,10 +859,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 @@ -764,32 +886,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..6120111 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 @@ -55,10 +55,10 @@ MODULE moduleMeshBoundary !Scatter particle in associated volume IF (ASSOCIATED(edge%e1)) THEN - CALL edge%e1%scatter(part) + CALL edge%e1%scatter(edge%e1%nNodes, part) ELSE - CALL edge%e2%scatter(part) + CALL edge%e2%scatter(edge%e2%nNodes, part) END IF @@ -156,11 +156,11 @@ MODULE moduleMeshBoundary newElectron%r = edge%randPos() newIon%r = newElectron%r - newElectron%vol = part%vol - newIon%vol = part%vol + newElectron%cell = part%cell + newIon%cell = part%cell - newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) - newIon%xi = newElectron%xi + newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r) + newIon%Xi = newElectron%Xi newElectron%weight = part%weight newIon%weight = newElectron%weight diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index ba1edfd..ccca930 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -111,13 +111,13 @@ MODULE moduleCollisions IMPLICIT NONE REAL(8):: n(1:3) - REAL(8):: cosXii, sinXii, eps + REAL(8):: cosXi, sinXi, eps - cosXii = random(-1.D0, 1.D0) - sinXii = DSQRT(1.D0 - cosXii**2) + cosXi = random(-1.D0, 1.D0) + sinXi = DSQRT(1.D0 - cosXi**2) eps = random(0.D0, PI2) - n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /) + n = (/ cosXi, sinXi*DCOS(eps), sinXi*DSIN(eps) /) END FUNCTION randomDirectionVHS @@ -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..daa3846 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") @@ -285,16 +285,16 @@ MODULE moduleInject partInj(n)%r = randomEdge%randPos() !Volume associated to the edge: IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%vol = randomEdge%e1%n + partInj(n)%cell = randomEdge%e1%n ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%vol = randomEdge%e2%n + partInj(n)%cell = randomEdge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - partInj(n)%volColl = randomEdge%eColl%n + partInj(n)%cellColl = randomEdge%eColl%n sp = self%species%n !Assign particle type @@ -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)%cell)%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/moduleList.f90 b/src/modules/moduleList.f90 index b040c80..b1dafdc 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -92,9 +92,12 @@ MODULE moduleList DEALLOCATE(current) current => next END DO + IF (ASSOCIATED(self%head)) NULLIFY(self%head) IF (ASSOCIATED(self%tail)) NULLIFY(self%tail) + self%amount = 0 + END SUBROUTINE eraseList SUBROUTINE setLock(self) diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index d19ff28..ab08f08 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -38,9 +38,9 @@ MODULE moduleSpecies REAL(8):: r(1:3) !Position REAL(8):: v(1:3) !Velocity CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle - INTEGER:: vol !Index of element in which the particle is located - INTEGER:: volColl !Index of element in which the particle is located in the Collision Mesh - REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p. + INTEGER:: cell !Index of element in which the particle is located + INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh + REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p. LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index ae73ebc..bdf6b03 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -26,12 +26,14 @@ MODULE moduleEM CLASS(boundaryEM), INTENT(in):: self CLASS(meshEdge):: edge + INTEGER:: nNodes INTEGER, ALLOCATABLE:: nodes(:) INTEGER:: n - nodes = edge%getNodes() + nNodes = edge%nNodes + nodes = edge%getNodes(nNodes) - DO n = 1, SIZE(nodes) + DO n = 1, nNodes SELECT CASE(self%typeEM) CASE ("dirichlet") mesh%K(nodes(n), :) = 0.D0 @@ -46,40 +48,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,9 +67,9 @@ MODULE moduleEM !$OMP END SINGLE !$OMP DO REDUCTION(+:vectorF) - DO e = 1, mesh%numVols - nodes = mesh%vols(e)%obj%getNodes() - nNodes = SIZE(nodes) + DO e = 1, mesh%numCells + nNodes = mesh%cells(e)%obj%nNodes + nodes = mesh%cells(e)%obj%getNodes(nNodes) !Calculates charge density (rho) in element nodes ALLOCATE(rho(1:nNodes)) rho = 0.D0 @@ -113,7 +81,7 @@ MODULE moduleEM END DO !Calculates local F vector - localF = mesh%vols(e)%obj%elemF(rho) + localF = mesh%cells(e)%obj%elemF(nNodes, 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..e557495 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -43,14 +43,14 @@ MODULE moduleSolver END SUBROUTINE solveEM_interface !Apply nonAnalogue scheme to a particle - SUBROUTINE weightingScheme_interface(part, volOld, volNew) + SUBROUTINE weightingScheme_interface(part, cellOld, cellNew) USE moduleSpecies USE moduleMesh IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: cellOld + CLASS(meshCell), POINTER, INTENT(inout):: cellNew END SUBROUTINE weightingScheme_interface @@ -204,7 +204,10 @@ MODULE moduleSolver DO n = 1, partList%amount partNext => partCurr%next partArray(nStart + n) = partCurr%part - CALL doProbes(partArray(nStart+n)) + IF (nProbes > 0) THEN + CALL doProbes(partArray(nStart+n)) + + END IF DEALLOCATE(partCurr) partCurr => partNext @@ -270,7 +273,10 @@ MODULE moduleSolver IF (partInj(n)%n_in) THEN nn = nn + 1 partOld(nn) = partInj(n) - CALL doProbes(partOld(nn)) + IF (nProbes > 0) THEN + CALL doProbes(partOld(nn)) + + END IF END IF @@ -283,7 +289,10 @@ MODULE moduleSolver IF (partTemp(n)%n_in) THEN nn = nn + 1 partOld(nn) = partTemp(n) - CALL doProbes(partOld(nn)) + IF (nProbes > 0) THEN + CALL doProbes(partOld(nn)) + + END IF END IF @@ -313,31 +322,37 @@ 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 - IF (solver%pusher(s)%pushSpecies) THEN - CALL mesh%vols(e)%obj%listPart_in(s)%erase() - mesh%vols(e)%obj%totalWeight(s) = 0.D0 + IF (doMCC) THEN + DO s = 1, nSpecies + DO e = 1, mesh%numCells + IF (solver%pusher(s)%pushSpecies) THEN + CALL mesh%cells(e)%obj%listPart_in(s)%erase() + mesh%cells(e)%obj%totalWeight(s) = 0.D0 - END IF + END IF + + END DO END DO - END DO + END IF !$OMP SECTION !Erase the list of particles inside the cell in coll mesh - DO s = 1, nSpecies - DO e = 1, meshColl%numVols - IF (solver%pusher(s)%pushSpecies) THEN - CALL meshColl%vols(e)%obj%listPart_in(s)%erase() - meshColl%vols(e)%obj%totalWeight(s) = 0.D0 + IF (doubleMesh) THEN + DO s = 1, nSpecies + DO e = 1, meshColl%numCells + IF (solver%pusher(s)%pushSpecies) THEN + CALL meshColl%cells(e)%obj%listPart_in(s)%erase() + meshColl%cells(e)%obj%totalWeight(s) = 0.D0 - END IF + END IF + + END DO END DO - - END DO + + END IF !$OMP END SECTIONS @@ -354,11 +369,13 @@ MODULE moduleSolver IMPLICIT NONE INTEGER:: n + CLASS(meshCell), POINTER:: cell !Loops over the particles to scatter them !$OMP DO DO n = 1, nPartOld - CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) + cell => mesh%cells(partOld(n)%cell)%obj + CALL cell%scatter(cell%nNodes, partOld(n)) END DO !$OMP END DO @@ -376,28 +393,28 @@ MODULE moduleSolver END SUBROUTINE doEMField !Split particles as a function of cell volume and splits particle - SUBROUTINE volumeWScheme(part, volOld, volNew) + SUBROUTINE volumeWScheme(part, cellOld, cellNew) USE moduleSpecies USE moduleMesh USE moduleRandom IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: cellOld + CLASS(meshCell), POINTER, INTENT(inout):: cellNew REAL(8):: fractionVolume, pSplit !If particle changes volume to smaller cell - IF (volOld%volume > volNew%volume .AND. & + IF (cellOld%volume > cellNew%volume .AND. & part%weight >= part%species%weight*1.0D-1) THEN - fractionVolume = volOld%volume/volNew%volume + fractionVolume = cellOld%volume/cellNew%volume !Calculate probability of splitting particle pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1) IF (random() < pSplit) THEN !Split particle in two - CALL splitParticle(part, 2, volNew) + CALL splitParticle(part, 2, cellNew) END IF @@ -407,7 +424,7 @@ MODULE moduleSolver !Subroutine to split the particle 'part' into a number 'nSplit' of particles. !'nSplit-1' particles are added to the partNAScheme list - SUBROUTINE splitParticle(part, nSplit, vol) + SUBROUTINE splitParticle(part, nSplit, cell) USE moduleSpecies USE moduleList USE moduleMesh @@ -416,7 +433,7 @@ MODULE moduleSolver TYPE(particle), INTENT(inout):: part INTEGER, INTENT(in):: nSplit - CLASS(meshVol), INTENT(inout):: vol + CLASS(meshCell), INTENT(inout):: cell REAL(8):: newWeight TYPE(particle), POINTER:: newPart INTEGER:: p @@ -438,10 +455,13 @@ MODULE moduleSolver CALL partWScheme%add(newPart) CALL partWScheme%unsetLock() !Add particle to cell list - CALL OMP_SET_lock(vol%lock) sp = part%species%n - CALL vol%listPart_in(sp)%add(newPart) - CALL OMP_UNSET_lock(vol%lock) + IF (doMCC) THEN + CALL OMP_SET_lock(cell%lock) + CALL cell%listPart_in(sp)%add(newPart) + CALL OMP_UNSET_lock(cell%lock) + + END IF END DO @@ -454,18 +474,18 @@ MODULE moduleSolver CLASS(solverGeneric), INTENT(in):: self TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER:: volOld, volNew + CLASS(meshCell), POINTER:: cellOld, cellNew !Assume that particle is outside the domain part%n_in = .FALSE. - volOld => mesh%vols(part%vol)%obj - CALL volOld%findCell(part) + cellOld => mesh%cells(part%cell)%obj + CALL cellOld%findCell(part) CALL findCellColl(part) - volNew => mesh%vols(part%vol)%obj !Call the NA shcme IF (ASSOCIATED(self%weightingScheme)) THEN - CALL self%weightingScheme(part, volOld, volNew) + cellNew => mesh%cells(part%cell)%obj + CALL self%weightingScheme(part, cellOld, cellNew) END IF diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90 index 69fcaab..bc02912 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%cell)%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%cell)%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%cell)%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%cell)%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%cell)%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