JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'improvement/performance' into 'development'

Browse source 
Performance improvement.

See merge request JorgeGonz/fpakc!33
This commit is contained in:
Jorge Gonzalez 2023-01-07 17:47:29 +00:00
commit d8eb2c2c58
23 changed files with 2522 additions and 2623 deletions
Download patch file Download diff file Expand all files Collapse all files

2
runs/1D_Cathode/inputRad.json
View file

@ -34,7 +34,7 @@
{"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} {"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
], ],
"inject": [ "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} "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}
], ],
"solver": { "solver": {

10
src/fpakc.f90
View file

@ -74,7 +74,10 @@ PROGRAM fpakc
tColl = omp_get_wtime() tColl = omp_get_wtime()
!$OMP END SINGLE !$OMP END SINGLE
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions(t) IF (doMCC) THEN
CALL meshForMCC%doCollisions(t)
END IF
!$OMP SINGLE !$OMP SINGLE
tColl = omp_get_wtime() - tColl tColl = omp_get_wtime() - tColl
@ -83,7 +86,10 @@ PROGRAM fpakc
tCoul = omp_get_wTime() tCoul = omp_get_wTime()
!$OMP END SINGLE !$OMP END SINGLE
IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb() IF (ASSOCIATED(mesh%doCoulomb)) THEN
CALL mesh%doCoulomb()
END IF
!$OMP SINGLE !$OMP SINGLE
tCoul = omp_get_wTime() - tCoul tCoul = omp_get_wTime() - tCoul

3
src/modules/common/moduleConstParam.f90
View file

@ -6,7 +6,8 @@ MODULE moduleConstParam
REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi
REAL(8), PARAMETER:: PI2 = 2.D0*PI !2*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:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1
REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge
REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI

106
src/modules/init/moduleInput.f90
View file

@ -331,13 +331,14 @@ MODULE moduleInput
REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
INTEGER, ALLOCATABLE, DIMENSION(:):: nodes INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
INTEGER:: nNodes INTEGER:: nNodes
REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi REAL(8), ALLOCATABLE, DIMENSION(:):: sourceScalar
REAL(8), ALLOCATABLE, DIMENSION(:,:):: sourceArray
!Density at the volume centroid !Density at the volume centroid
REAL(8):: densityCen REAL(8):: densityCen
!Mean velocity and temperature at particle position !Mean velocity and temperature at particle position
REAL(8):: velocityXi(1:3), temperatureXi REAL(8):: velocityXi(1:3), temperatureXi
INTEGER:: nNewPart = 0.D0 INTEGER:: nNewPart = 0.D0
CLASS(meshVol), POINTER:: vol CLASS(meshCell), POINTER:: cell
TYPE(particle), POINTER:: partNew TYPE(particle), POINTER:: partNew
REAL(8):: vTh REAL(8):: vTh
TYPE(lNode), POINTER:: partCurr, partNext TYPE(lNode), POINTER:: partCurr, partNext
@ -354,87 +355,74 @@ MODULE moduleInput
CALL config%get(object // '.file', spFile, found) CALL config%get(object // '.file', spFile, found)
!Reads node values at the nodes !Reads node values at the nodes
filename = path // spFile 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 !For each volume in the node, create corresponding particles
DO e = 1, mesh%numVols DO e = 1, mesh%numCells
!Scale variables !Scale variables
!Density at centroid of cell !Density at centroid of cell
nodes = mesh%vols(e)%obj%getNodes() nNodes = mesh%cells(e)%obj%nNodes
nNodes = SIZE(nodes) nodes = mesh%cells(e)%obj%getNodes(nNodes)
fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) ALLOCATE(sourceScalar(1:nNodes))
ALLOCATE(source(1:nNodes)) ALLOCATE(sourceArray(1:nNodes, 1:3))
DO j = 1, nNodes sourceScalar = density(nodes)
source(j) = density(nodes(j)) densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)
END DO
densityCen = DOT_PRODUCT(fPsi, source)
DEALLOCATE(fPsi)
!Calculate number of particles !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 !Allocate new particles
DO p = 1, nNewPart DO p = 1, nNewPart
ALLOCATE(partNew) ALLOCATE(partNew)
partNew%species => species(sp)%obj partNew%species => species(sp)%obj
partNew%r = mesh%vols(e)%obj%randPos() partNew%r = mesh%cells(e)%obj%randPos()
partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r)
!Get mean velocity at particle position !Get mean velocity at particle position
fPsi = mesh%vols(e)%obj%fPsi(partNew%xi) sourceArray(:,1) = velocity((nodes), 1)
DO j = 1, nNodes sourceArray(:,2) = velocity((nodes), 2)
source(j) = velocity(nodes(j), 1) sourceArray(:,3) = velocity((nodes), 3)
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)
velocityXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceArray)
velocityXi = velocityXi / v_ref velocityXi = velocityXi / v_ref
!Get temperature at particle position !Get temperature at particle position
DO j = 1, nNodes sourceScalar = temperature(nodes)
source(j) = temperature(nodes(j)) temperatureXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceScalar)
END DO
temperatureXi = DOT_PRODUCT(fPsi, source)
temperatureXi = temperatureXi / T_ref temperatureXi = temperatureXi / T_ref
vTh = DSQRT(temperatureXi / species(sp)%obj%m) vTh = DSQRT(temperatureXi / species(sp)%obj%m)
partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian() partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian()
partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian()
partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian()
partNew%vol = e
partNew%cell = e
IF (doubleMesh) THEN IF (doubleMesh) THEN
partNew%volColl = findCellBrute(meshColl, partNew%r) partNew%cellColl = findCellBrute(meshColl, partNew%r)
ELSE ELSE
partNew%volColl = partNew%vol partNew%cellColl = partNew%cell
END IF END IF
partNew%n_in = .TRUE. partNew%n_in = .TRUE.
partNew%weight = species(sp)%obj%weight partNew%weight = species(sp)%obj%weight
!Assign particle to temporal list of particles !Assign particle to temporal list of particles
CALL partInitial%add(partNew) CALL partInitial%add(partNew)
!Assign particle to list in volume !Assign particle to list in volume
vol => meshforMCC%vols(partNew%volColl)%obj IF (doMCC) THEN
CALL OMP_SET_LOCK(vol%lock) cell => meshforMCC%cells(partNew%cellColl)%obj
CALL vol%listPart_in(sp)%add(partNew) CALL OMP_SET_LOCK(cell%lock)
vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight CALL cell%listPart_in(sp)%add(partNew)
CALL OMP_UNSET_LOCK(vol%lock) cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight
CALL OMP_UNSET_LOCK(cell%lock)
END IF
END DO END DO
DEALLOCATE(source) DEALLOCATE(sourceScalar, sourceArray)
END DO END DO
@ -643,9 +631,9 @@ MODULE moduleInput
REAL(8):: energyThreshold, energyBinding REAL(8):: energyThreshold, energyBinding
CHARACTER(:), ALLOCATABLE:: electron CHARACTER(:), ALLOCATABLE:: electron
INTEGER:: e 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) CALL config%info('interactions', found)
IF (found) THEN IF (found) THEN
!Checks if MC collisions have been defined !Checks if MC collisions have been defined
@ -739,18 +727,18 @@ MODULE moduleInput
END DO END DO
!Init the required arrays in each volume to account for MCC. !Init the required arrays in each volume to account for MCC.
DO e = 1, meshForMCC%numVols DO e = 1, meshForMCC%numCells
vol => meshForMCC%vols(e)%obj cell => meshForMCC%cells(e)%obj
!Allocate Maximum cross section per collision pair and assign the initial collision rate !Allocate Maximum cross section per collision pair and assign the initial collision rate
ALLOCATE(vol%sigmaVrelMax(1:nCollPairs)) ALLOCATE(cell%sigmaVrelMax(1:nCollPairs))
vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
IF (collOutput) THEN IF (collOutput) THEN
ALLOCATE(vol%tallyColl(1:nCollPairs)) ALLOCATE(cell%tallyColl(1:nCollPairs))
DO k = 1, nCollPairs DO k = 1, nCollPairs
ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
vol%tallyColl(k)%tally = 0 cell%tallyColl(k)%tally = 0
END DO END DO
@ -907,6 +895,8 @@ MODULE moduleInput
END IF END IF
doMCC = ASSOCIATED(meshForMCC)
!Get the dimension of the geometry !Get the dimension of the geometry
CALL config%get(object // '.dimension', mesh%dimen, found) CALL config%get(object // '.dimension', mesh%dimen, found)
IF (.NOT. found) THEN IF (.NOT. found) THEN
@ -930,8 +920,8 @@ MODULE moduleInput
CALL config%get(object // '.volume', volume, found) CALL config%get(object // '.volume', volume, found)
!Rescale the volumne !Rescale the volumne
IF (found) THEN IF (found) THEN
mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref
mesh%nodes(1)%obj%v = mesh%vols(1)%obj%volume mesh%nodes(1)%obj%v = mesh%cells(1)%obj%volume
END IF END IF

160
src/modules/mesh/0D/moduleMesh0D.f90
View file

@ -6,27 +6,32 @@ MODULE moduleMesh0D
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D
INTEGER:: n1 INTEGER:: n1
CONTAINS CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode0D PROCEDURE, PASS:: init => initNode0D
PROCEDURE, PASS:: getCoordinates => getCoord0D PROCEDURE, PASS:: getCoordinates => getCoord0D
END TYPE meshNode0D END TYPE meshNode0D
TYPE, PUBLIC, EXTENDS(meshVol):: meshVol0D TYPE, PUBLIC, EXTENDS(meshCell):: meshCell0D
CLASS(meshNode), POINTER:: n1 CLASS(meshNode), POINTER:: n1
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initVol0D PROCEDURE, PASS:: init => initCell0D
PROCEDURE, PASS:: getNodes => getNodes0D PROCEDURE, PASS:: getNodes => getNodes0D
PROCEDURE, PASS:: randPos => randPos0D PROCEDURE, PASS:: randPos => randPos0D
PROCEDURE, NOPASS:: fPsi => fPsi0D PROCEDURE, NOPASS:: fPsi => fPsi0D
PROCEDURE, PASS:: gatherEF => gatherEF0D PROCEDURE, NOPASS:: dPsi => dPsi0D
PROCEDURE, PASS:: gatherMF => gatherMF0D 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:: elemK => elemK0D
PROCEDURE, PASS:: elemF => elemF0D PROCEDURE, PASS:: elemF => elemF0D
PROCEDURE, PASS:: phy2log => phy2log0D PROCEDURE, PASS:: phy2log => phy2log0D
PROCEDURE, NOPASS:: inside => inside0D PROCEDURE, NOPASS:: inside => inside0D
PROCEDURE, PASS:: nextElement => nextElement0D PROCEDURE, PASS:: neighbourElement => neighbourElement0D
END TYPE meshVol0D END TYPE meshCell0D
CONTAINS CONTAINS
!NODE FUNCTIONS !NODE FUNCTIONS
@ -61,18 +66,20 @@ MODULE moduleMesh0D
!VOLUME FUNCTIONS !VOLUME FUNCTIONS
!Inits dummy 0D volume !Inits dummy 0D volume
SUBROUTINE initVol0D(self, n, p, nodes) SUBROUTINE initCell0D(self, n, p, nodes)
USE moduleRefParam USE moduleRefParam
USE moduleSpecies USE moduleSpecies
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(out):: self CLASS(meshCell0D), INTENT(out):: self
INTEGER, INTENT(in):: n INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => nodes(p(1))%obj self%n1 => nodes(p(1))%obj
self%volume = 1.D0 self%volume = 1.D0
self%n1%v = 1.D0 self%n1%v = 1.D0
@ -82,88 +89,122 @@ MODULE moduleMesh0D
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(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 IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
ALLOCATE(n(1:1))
n = self%n1%n n = self%n1%n
END FUNCTION getNodes0D END FUNCTION getNodes0D
!Calculate random position inside the volume
FUNCTION randPos0D(self) RESULT(r) FUNCTION randPos0D(self) RESULT(r)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
r = 0.D0 r = 0.D0
END FUNCTION randPos0D END FUNCTION randPos0D
PURE FUNCTION fPsi0D(xi) RESULT(fPsi) PURE FUNCTION fPsi0D(Xi, nNodes) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3) IMPLICIT NONE
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
ALLOCATE(fPsi(1:1))
fPsi = 1.D0 fPsi = 1.D0
END FUNCTION fPsi0D END FUNCTION fPsi0D
PURE FUNCTION gatherEF0D(self, xi) RESULT(EF) PURE FUNCTION dPsi0D(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in):: xi(1:3) INTEGER, INTENT(in):: nNodes
REAL(8):: EF(1:3) 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 IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) INTEGER, INTENT(in):: nNodes
REAL(8):: MF(1:3) 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 IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:) INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0 localK = 0.D0
END FUNCTION elemK0D END FUNCTION elemK0D
PURE FUNCTION elemF0D(self, source) RESULT(localF) PURE FUNCTION elemF0D(self, nNodes, source) RESULT(localF)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: source(1:) INTEGER, INTENT(in):: nNodes
REAL(8), ALLOCATABLE:: localF(:) REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
ALLOCATE(localF(1:1))
localF = 0.D0 localF = 0.D0
END FUNCTION elemF0D 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) PURE FUNCTION phy2log0D(self,r) RESULT(xN)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3) REAL(8):: xN(1:3)
@ -171,25 +212,46 @@ MODULE moduleMesh0D
END FUNCTION phy2log0D END FUNCTION phy2log0D
PURE FUNCTION inside0D(xi) RESULT(ins) PURE FUNCTION inside0D(Xi) RESULT(ins)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins LOGICAL:: ins
ins = .TRUE. ins = .TRUE.
END FUNCTION inside0D END FUNCTION inside0D
SUBROUTINE nextElement0D(self, xi, nextElement) SUBROUTINE neighbourElement0D(self, Xi, neighbourElement)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol0D), INTENT(in):: self CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement 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 END MODULE moduleMesh0D

480
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -14,6 +14,7 @@ MODULE moduleMesh1DCart
!Element coordinates !Element coordinates
REAL(8):: x = 0.D0 REAL(8):: x = 0.D0
CONTAINS CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode1DCart PROCEDURE, PASS:: init => initNode1DCart
PROCEDURE, PASS:: getCoordinates => getCoord1DCart PROCEDURE, PASS:: getCoordinates => getCoord1DCart
@ -25,6 +26,7 @@ MODULE moduleMesh1DCart
!Connectivity to nodes !Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL() CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge1DCart PROCEDURE, PASS:: init => initEdge1DCart
PROCEDURE, PASS:: getNodes => getNodes1DCart PROCEDURE, PASS:: getNodes => getNodes1DCart
PROCEDURE, PASS:: intersection => intersection1DCart PROCEDURE, PASS:: intersection => intersection1DCart
@ -32,57 +34,34 @@ MODULE moduleMesh1DCart
END TYPE meshEdge1DCart END TYPE meshEdge1DCart
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DCartSegm
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
!Element coordinates !Element coordinates
REAL(8):: x(1:2) REAL(8):: x(1:2)
!Connectivity to nodes !Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
!Connectivity to adjacent elements !Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initVol1DCartSegm !meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: randPos => randPos1DCartSeg PROCEDURE, PASS:: init => initCell1DCartSegm
PROCEDURE, PASS:: area => areaSegm PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: randPos => randPos1DCartSegm
PROCEDURE, NOPASS:: fPsi => fPsiSegm PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, NOPASS:: dPsi => dPsiSegm PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: partialDer => partialDerSegm 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:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: inside => insideSegm PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: gatherEF => gatherEFSegm
PROCEDURE, PASS:: gatherMF => gatherMFSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: nextElement => nextElementSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm
!PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm
END TYPE meshVol1DCartSegm END TYPE meshCell1DCartSegm
CONTAINS CONTAINS
!NODE FUNCTIONS !NODE FUNCTIONS
@ -120,7 +99,7 @@ MODULE moduleMesh1DCart
END FUNCTION getCoord1DCart END FUNCTION getCoord1DCart
!EDGE FUNCTIONS !EDGE FUNCTIONS
!Inits edge element !Init edge element
SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface) SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface)
USE moduleSpecies USE moduleSpecies
USE moduleBoundary USE moduleBoundary
@ -136,6 +115,7 @@ MODULE moduleMesh1DCart
INTEGER:: s INTEGER:: s
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => mesh%nodes(p(1))%obj self%n1 => mesh%nodes(p(1))%obj
!Get element coordinates !Get element coordinates
r1 = self%n1%getCoordinates() r1 = self%n1%getCoordinates()
@ -159,13 +139,13 @@ MODULE moduleMesh1DCart
END SUBROUTINE initEdge1DCart END SUBROUTINE initEdge1DCart
!Get nodes from edge !Get nodes from edge
PURE FUNCTION getNodes1DCart(self) RESULT(n) PURE FUNCTION getNodes1DCart(self, nNodes) RESULT(n)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge1DCart), INTENT(in):: self CLASS(meshEdge1DCart), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
ALLOCATE(n(1))
n = (/ self%n1%n /) n = (/ self%n1%n /)
END FUNCTION getNodes1DCart END FUNCTION getNodes1DCart
@ -181,7 +161,7 @@ MODULE moduleMesh1DCart
END FUNCTION intersection1DCart END FUNCTION intersection1DCart
!Calculates a 'random' position in edge !Calculate a 'random' position in edge
FUNCTION randPosEdge(self) RESULT(r) FUNCTION randPosEdge(self) RESULT(r)
CLASS(meshEdge1DCart), INTENT(in):: self CLASS(meshEdge1DCart), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
@ -192,18 +172,19 @@ MODULE moduleMesh1DCart
!VOLUME FUNCTIONS !VOLUME FUNCTIONS
!SEGMENT FUNCTIONS !SEGMENT FUNCTIONS
!Init segment element !Init element
SUBROUTINE initVol1DCartSegm(self, n, p, nodes) SUBROUTINE initCell1DCartSegm(self, n, p, nodes)
USE moduleRefParam USE moduleRefParam
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(out):: self CLASS(meshCell1DCartSegm), INTENT(out):: self
INTEGER, INTENT(in):: n INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2 REAL(8), DIMENSION(1:3):: r1, r2
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => nodes(p(1))%obj self%n1 => nodes(p(1))%obj
self%n2 => nodes(p(2))%obj self%n2 => nodes(p(2))%obj
!Get element coordinates !Get element coordinates
@ -212,144 +193,182 @@ MODULE moduleMesh1DCart
self%x = (/ r1(1), r2(1) /) self%x = (/ r1(1), r2(1) /)
!Assign node volume !Assign node volume
CALL self%area() CALL self%calculateVolume()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol1DCartSegm END SUBROUTINE initCell1DCartSegm
!Calculates a random position in 1D volume !Get nodes from 1D volume
FUNCTION randPos1DCartSeg(self) RESULT(r) 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 USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
REAL(8):: xii(1:3) REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8):: fPsi(1:2)
xii(1) = random(-1.D0, 1.D0) Xi = 0.D0
xii(2:3) = 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) r(1) = DOT_PRODUCT(fPsi, self%x)
END FUNCTION randPos1DCartSeg END FUNCTION randPos1DCartSegm
!Computes element area !Compute element functions at point Xi
PURE SUBROUTINE areaSegm(self) PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(inout):: self REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: l !element length INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:2) REAL(8):: fPsi(1:nNodes)
REAL(8):: detJ
REAL(8):: Xii(1:3)
self%volume = 0.D0 fPsi = (/ 1.D0 - Xi(1), &
self%arNodes = 0.D0 1.D0 + Xi(1) /)
!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
END SUBROUTINE areaSegm fPsi = fPsi * 0.50D0
!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
END FUNCTION fPsiSegm END FUNCTION fPsiSegm
!Computes element derivative shape function at Xii !Derivative element function at coordinates Xi
PURE FUNCTION dPsiSegm(xi) RESULT(dPsi) PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) 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, 1:2) = (/ -5.D-1, 5.D-1 /)
dPsi(1, 2) = 5.D-1
END FUNCTION dPsiSegm END FUNCTION dPsiSegm
!Computes partial derivatives of coordinates !Partial derivative in global coordinates
PURE SUBROUTINE partialDerSegm(self, dPsi, dx) PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:) INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(out), DIMENSION(1):: dx 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 END FUNCTION partialDerSegm
PURE FUNCTION elemKSegm(self) RESULT(localK)
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 IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:) INTEGER, INTENT(in):: nNodes
REAL(8):: Xii(1:3) REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: dPsi(1:1, 1:2) REAL(8):: Xi(1:3)
REAL(8):: invJ(1), detJ REAL(8):: dPsi(1:3, 1:2)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3, 1:3), detJ
INTEGER:: l INTEGER:: l
ALLOCATE(localK(1:2,1:2))
localK = 0.D0 localK = 0.D0
Xii = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3 DO l = 1, 3
xii(1) = corSeg(l) Xi(1) = corSeg(l)
dPsi = self%dPsi(Xii) dPsi = self%dPsi(Xi, 2)
detJ = self%detJac(Xii, dPsi) pDer = self%partialDer(2, dPsi)
invJ = self%invJac(Xii, dPsi) detJ = self%detJac(pDer)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & invJ = self%invJac(pDer)
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
wSeg(l)/detJ wSeg(l)/detJ
END DO END DO
END FUNCTION elemKSegm 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 IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: source(1:) INTEGER, INTENT(in):: nNodes
REAL(8), ALLOCATABLE:: localF(:) REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:2) 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):: detJ, f
REAL(8):: Xii(1:3)
INTEGER:: l INTEGER:: l
ALLOCATE(localF(1:2))
localF = 0.D0 localF = 0.D0
Xii = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3 DO l = 1, 3
xii(1) = corSeg(l) Xi(1) = corSeg(l)
detJ = self%detJac(Xii) dPsi = self%dPsi(Xi, 2)
fPsi = self%fPsi(Xii) pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 2)
f = DOT_PRODUCT(fPsi, source) f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*wSeg(l)*detJ localF = localF + f*fPsi*wSeg(l)*detJ
@ -357,151 +376,98 @@ MODULE moduleMesh1DCart
END FUNCTION elemFSegm END FUNCTION elemFSegm
PURE FUNCTION insideSegm(xi) RESULT(ins) PURE FUNCTION insideSegm(Xi) RESULT(ins)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins LOGICAL:: ins
ins = xi(1) >=-1.D0 .AND. & ins = Xi(1) >=-1.D0 .AND. &
xi(1) <= 1.D0 Xi(1) <= 1.D0
END FUNCTION insideSegm END FUNCTION insideSegm
!Gathers EF at position Xii PURE FUNCTION phy2logSegm(self, r) RESULT(Xi)
PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), 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
REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3) REAL(8):: Xi(1:3)
xN = 0.D0 Xi = 0.D0
xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
Xi(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
END FUNCTION phy2logSegm END FUNCTION phy2logSegm
!Get next element for a logical position xi !Get the next element for a logical position Xi
SUBROUTINE nextElementSegm(self, xi, nextElement) SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(in):: self CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
NULLIFY(nextElement) NULLIFY(neighbourElement)
IF (xi(1) < -1.D0) THEN IF (Xi(1) < -1.D0) THEN
nextElement => self%e2 neighbourElement => self%e2
ELSEIF (xi(1) > 1.D0) THEN ELSEIF (Xi(1) > 1.D0) THEN
nextElement => self%e1 neighbourElement => self%e1
END IF 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 !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Calculates a random position in 1D volume !Compute element Jacobian determinant
!Computes the element Jacobian determinant PURE FUNCTION detJ1DCart(pDer) RESULT(dJ)
PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCart), INTENT(in):: self REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ REAL(8):: dJ
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN dJ = pDer(1, 1)
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
dJ = dx(1)
END FUNCTION detJ1DCart END FUNCTION detJ1DCart
!Computes the invers Jacobian !Compute element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ) PURE FUNCTION invJ1DCart(pDer) RESULT(invJ)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCart), INTENT(in):: self REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8), INTENT(in):: xi(1:3) REAL(8):: invJ(1:3,1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1)
REAL(8):: invJ
IF (PRESENT(dPsi_in)) THEN invJ = 0.D0
dPsi = dPsi_in
ELSE invJ(1, 1) = 1.D0/pDer(1, 1)
dPsi = self%dPsi(xi) invJ(2, 2) = 1.D0
invJ(3, 3) = 1.D0
END IF
CALL self%partialDer(dPsi, dx)
invJ = 1.D0/dx(1)
END FUNCTION invJ1DCart END FUNCTION invJ1DCart
@ -511,11 +477,11 @@ MODULE moduleMesh1DCart
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER:: e, et INTEGER:: e, et
DO e = 1, self%numVols DO e = 1, self%numCells
!Connect Vol-Vol !Connect Cell-Cell
DO et = 1, self%numVols DO et = 1, self%numCells
IF (e /= et) THEN 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 END IF
@ -523,9 +489,9 @@ MODULE moduleMesh1DCart
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
!Connect Vol-Edge !Connect Cell-Edge
DO et = 1, self%numEdges 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 END DO
@ -535,29 +501,29 @@ MODULE moduleMesh1DCart
END SUBROUTINE connectMesh1DCart END SUBROUTINE connectMesh1DCart
SUBROUTINE connectVolVol(elemA, elemB) SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB CLASS(meshCell), INTENT(inout):: elemB
SELECT TYPE(elemA) SELECT TYPE(elemA)
TYPE IS(meshVol1DCartSegm) TYPE IS(meshCell1DCartSegm)
SELECT TYPE(elemB) SELECT TYPE(elemB)
TYPE IS(meshVol1DCartSegm) TYPE IS(meshCell1DCartSegm)
CALL connectSegmSegm(elemA, elemB) CALL connectSegmSegm(elemA, elemB)
END SELECT END SELECT
END SELECT END SELECT
END SUBROUTINE connectVolVol END SUBROUTINE connectCellCell
SUBROUTINE connectSegmSegm(elemA, elemB) SUBROUTINE connectSegmSegm(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemA
CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemB CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. & IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN elemA%n2%n == elemB%n1%n) THEN
@ -576,14 +542,14 @@ MODULE moduleMesh1DCart
END SUBROUTINE connectSegmSegm END SUBROUTINE connectSegmSegm
SUBROUTINE connectVolEdge(elemA, elemB) SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA) SELECT TYPE(elemA)
TYPE IS (meshVol1DCartSegm) TYPE IS (meshCell1DCartSegm)
SELECT TYPE(elemB) SELECT TYPE(elemB)
CLASS IS(meshEdge1DCart) CLASS IS(meshEdge1DCart)
CALL connectSegmEdge(elemA, elemB) CALL connectSegmEdge(elemA, elemB)
@ -592,12 +558,12 @@ MODULE moduleMesh1DCart
END SELECT END SELECT
END SUBROUTINE connectVolEdge END SUBROUTINE connectCellEdge
SUBROUTINE connectSegmEdge(elemA, elemB) SUBROUTINE connectSegmEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DCartSegm), INTENT(inout), TARGET:: elemA CLASS(meshCell1DCartSegm), INTENT(inout), TARGET:: elemA
CLASS(meshEdge1DCart), INTENT(inout), TARGET:: elemB CLASS(meshEdge1DCart), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. & IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
@ -606,7 +572,7 @@ MODULE moduleMesh1DCart
elemA%e1 => elemB elemA%e1 => elemB
elemB%e2 => elemA elemB%e2 => elemA
!Revers the normal to point inside the domain !Rever the normal to point inside the domain
elemB%normal = - elemB%normal elemB%normal = - elemB%normal
END IF END IF

563
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -14,6 +14,7 @@ MODULE moduleMesh1DRad
!Element coordinates !Element coordinates
REAL(8):: r = 0.D0 REAL(8):: r = 0.D0
CONTAINS CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode1DRad PROCEDURE, PASS:: init => initNode1DRad
PROCEDURE, PASS:: getCoordinates => getCoord1DRad PROCEDURE, PASS:: getCoordinates => getCoord1DRad
@ -25,6 +26,7 @@ MODULE moduleMesh1DRad
!Connectivity to nodes !Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL() CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge1DRad PROCEDURE, PASS:: init => initEdge1DRad
PROCEDURE, PASS:: getNodes => getNodes1DRad PROCEDURE, PASS:: getNodes => getNodes1DRad
PROCEDURE, PASS:: intersection => intersection1DRad PROCEDURE, PASS:: intersection => intersection1DRad
@ -32,58 +34,34 @@ MODULE moduleMesh1DRad
END TYPE meshEdge1DRad END TYPE meshEdge1DRad
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DRadSegm
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
!Element coordinates !Element coordinates
REAL(8):: r(1:2) REAL(8):: r(1:2)
!Connectivity to nodes !Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
!Connectivity to adjacent elements !Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initVol1DRadSegm !meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: randPos => randPos1DRadSeg PROCEDURE, PASS:: init => initCell1DRadSegm
PROCEDURE, PASS:: area => areaRad PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, NOPASS:: fPsi => fPsiRad PROCEDURE, PASS:: randPos => randPos1DRadSegm
PROCEDURE, NOPASS:: dPsi => dPsiRad PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, PASS:: partialDer => partialDerRad PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: elemK => elemKRad PROCEDURE, PASS:: partialDer => partialDerSegm
PROCEDURE, PASS:: elemF => elemFRad PROCEDURE, NOPASS:: detJac => detJ1DRad
PROCEDURE, NOPASS:: inside => insideRad PROCEDURE, NOPASS:: invJac => invJ1DRad
PROCEDURE, PASS:: gatherEF => gatherEFRad PROCEDURE, PASS:: gatherElectricField => gatherEFSegm
PROCEDURE, PASS:: gatherMF => gatherMFRad PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm
PROCEDURE, PASS:: getNodes => getNodesRad PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: phy2log => phy2logRad PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, PASS:: nextElement => nextElementRad 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 CONTAINS
!NODE FUNCTIONS !NODE FUNCTIONS
@ -103,7 +81,7 @@ MODULE moduleMesh1DRad
!Node volume, to be determined in mesh !Node volume, to be determined in mesh
self%v = 0.D0 self%v = 0.D0
!Allocates output !Allocate output
ALLOCATE(self%output(1:nSpecies)) ALLOCATE(self%output(1:nSpecies))
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
@ -121,7 +99,7 @@ MODULE moduleMesh1DRad
END FUNCTION getCoord1DRad END FUNCTION getCoord1DRad
!EDGE FUNCTIONS !EDGE FUNCTIONS
!Inits edge element !Init edge element
SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface) SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface)
USE moduleSpecies USE moduleSpecies
USE moduleBoundary USE moduleBoundary
@ -137,6 +115,7 @@ MODULE moduleMesh1DRad
INTEGER:: s INTEGER:: s
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => mesh%nodes(p(1))%obj self%n1 => mesh%nodes(p(1))%obj
!Get element coordinates !Get element coordinates
r1 = self%n1%getCoordinates() r1 = self%n1%getCoordinates()
@ -144,7 +123,6 @@ MODULE moduleMesh1DRad
self%r = r1(1) self%r = r1(1)
self%normal = (/ 1.D0, 0.D0, 0.D0 /) self%normal = (/ 1.D0, 0.D0, 0.D0 /)
self%normal = self%normal/NORM2(self%normal)
!Boundary index !Boundary index
self%boundary => boundary(bt) self%boundary => boundary(bt)
@ -161,13 +139,13 @@ MODULE moduleMesh1DRad
END SUBROUTINE initEdge1DRad END SUBROUTINE initEdge1DRad
!Get nodes from edge !Get nodes from edge
PURE FUNCTION getNodes1DRad(self) RESULT(n) PURE FUNCTION getNodes1DRad(self, nNodes) RESULT(n)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge1DRad), INTENT(in):: self CLASS(meshEdge1DRad), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
ALLOCATE(n(1))
n = (/ self%n1%n /) n = (/ self%n1%n /)
END FUNCTION getNodes1DRad END FUNCTION getNodes1DRad
@ -183,7 +161,7 @@ MODULE moduleMesh1DRad
END FUNCTION intersection1DRad END FUNCTION intersection1DRad
!Calculates a 'random' position in edge !Calculate a 'random' position in edge
FUNCTION randPos1DRad(self) RESULT(r) FUNCTION randPos1DRad(self) RESULT(r)
CLASS(meshEdge1DRad), INTENT(in):: self CLASS(meshEdge1DRad), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
@ -194,18 +172,19 @@ MODULE moduleMesh1DRad
!VOLUME FUNCTIONS !VOLUME FUNCTIONS
!SEGMENT FUNCTIONS !SEGMENT FUNCTIONS
!Init segment element !Init element
SUBROUTINE initVol1DRadSegm(self, n, p, nodes) SUBROUTINE initCell1DRadSegm(self, n, p, nodes)
USE moduleRefParam USE moduleRefParam
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(out):: self CLASS(meshCell1DRadSegm), INTENT(out):: self
INTEGER, INTENT(in):: n INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2 REAL(8), DIMENSION(1:3):: r1, r2
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => nodes(p(1))%obj self%n1 => nodes(p(1))%obj
self%n2 => nodes(p(2))%obj self%n2 => nodes(p(2))%obj
!Get element coordinates !Get element coordinates
@ -214,312 +193,296 @@ MODULE moduleMesh1DRad
self%r = (/ r1(1), r2(1) /) self%r = (/ r1(1), r2(1) /)
!Assign node volume !Assign node volume
CALL self%area() CALL self%calculateVolume()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol1DRadSegm END SUBROUTINE initCell1DRadSegm
!Calculates a random position in 1D volume !Get nodes from 1D volume
FUNCTION randPos1DRadSeg(self) RESULT(r) 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 USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
REAL(8):: xii(1:3) REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8):: fPsi(1:2)
xii(1) = random(-1.D0, 1.D0) Xi = 0.D0
xii(2:3) = 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) r(1) = DOT_PRODUCT(fPsi, self%r)
END FUNCTION randPos1DRadSeg END FUNCTION randPos1DRadSegm
!Computes element area !Compute element functions at point Xi
PURE SUBROUTINE areaRad(self) PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout):: self REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: l !element length INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:2) REAL(8):: fPsi(1:nNodes)
REAL(8):: r
REAL(8):: detJ
REAL(8):: Xii(1:3)
self%volume = 0.D0 fPsi = (/ 1.D0 - Xi(1), &
self%arNodes = 0.D0 1.D0 + Xi(1) /)
!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
END SUBROUTINE areaRad fPsi = fPsi * 0.50D0
!Computes element functions at point xii END FUNCTION fPsiSegm
PURE FUNCTION fPsiRad(xi) RESULT(fPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
ALLOCATE(fPsi(1:2)) dPsi = 0.D0
fPsi(1) = 1.D0 - xi(1) dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiRad END FUNCTION dPsiSegm
!Computes element derivative shape function at Xii !Partial derivative in global coordinates
PURE FUNCTION dPsiRad(xi) RESULT(dPsi) PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), ALLOCATABLE:: dPsi(:,:) 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 pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%r(1:2))
dPsi(1, 2) = 5.D-1 pDer(2,2) = 1.D0
pDer(3,3) = 1.D0
END FUNCTION dPsiRad END FUNCTION partialDerSegm
!Computes partial derivatives of coordinates PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array)
PURE SUBROUTINE partialDerRad(self, dPsi, dx)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
CLASS(meshVol1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in):: dPsi(1:,1:) REAL(8):: array(1:3)
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)
REAL(8):: phi(1:2) REAL(8):: phi(1:2)
REAL(8):: EF(1:3)
REAL(8):: invJ
phi = (/ self%n1%emData%phi, & phi = (/ self%n1%emData%phi, &
self%n2%emData%phi /) self%n2%emData%phi /)
dPsi = self%dPsi(xi) array = -self%gatherDF(Xi, 2, phi)
invJ = self%invJac(xi, dPsi)
EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
EF(2) = 0.D0
EF(3) = 0.D0
END FUNCTION gatherEFRad END FUNCTION gatherEFSegm
PURE FUNCTION gatherMFRad(self, xi) RESULT(MF) PURE FUNCTION gatherMFSegm(self, Xi) RESULT(array)
IMPLICIT NONE 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)
CLASS(meshVol1DRadSegm), INTENT(in):: self B(:,1) = (/ self%n1%emData%B(1), &
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) /) self%n2%emData%B(1) /)
MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
B(:,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2) /) self%n2%emData%B(2) /)
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
B(:,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /) self%n2%emData%B(3) /)
fPsi = self%fPsi(xi) array = self%gatherF(Xi, 2, B)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFRad END FUNCTION gatherMFSegm
!Get nodes from 1D volume !Compute element local stiffness matrix
PURE FUNCTION getNodesRad(self) RESULT(n) PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) 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
ALLOCATE(n(1:2)) localK = 0.D0
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesRad 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
PURE FUNCTION phy2logRad(self, r) RESULT(xN) 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 IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self 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):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
REAL(8):: r
INTEGER:: l
localF = 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)
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 DO
localF = localF*PI2
END FUNCTION elemFSegm
PURE FUNCTION insideSegm(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 insideSegm
PURE FUNCTION phy2logSegm(self, r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3) REAL(8):: Xi(1:3)
xN = 0.D0 Xi = 0.D0
xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.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 END FUNCTION phy2logSegm
SUBROUTINE nextElementRad(self, xi, nextElement)
!Get the next element for a logical position Xi
SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
NULLIFY(nextElement) NULLIFY(neighbourElement)
IF (xi(1) < -1.D0) THEN IF (Xi(1) < -1.D0) THEN
nextElement => self%e2 neighbourElement => self%e2
ELSEIF (xi(1) > 1.D0) THEN ELSEIF (Xi(1) > 1.D0) THEN
nextElement => self%e1 neighbourElement => self%e1
END IF 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 !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes the element Jacobian determinant !Compute element Jacobian determinant
PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ) PURE FUNCTION detJ1DRad(pDer) RESULT(dJ)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRad), INTENT(in):: self REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ REAL(8):: dJ
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN dJ = pDer(1, 1)
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
dJ = dx(1)
END FUNCTION detJ1DRad END FUNCTION detJ1DRad
!Computes the invers Jacobian !Compute element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ) PURE FUNCTION invJ1DRad(pDer) RESULT(invJ)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRad), INTENT(in):: self REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8), INTENT(in):: xi(1:3) REAL(8):: invJ(1:3,1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1)
REAL(8):: invJ
IF (PRESENT(dPsi_in)) THEN invJ = 0.D0
dPsi = dPsi_in
ELSE invJ(1, 1) = 1.D0/pDer(1, 1)
dPsi = self%dPsi(xi) invJ(2, 2) = 1.D0
invJ(3, 3) = 1.D0
END IF
CALL self%partialDer(dPsi, dx)
invJ = 1.D0/dx(1)
END FUNCTION invJ1DRad END FUNCTION invJ1DRad
@ -529,11 +492,11 @@ MODULE moduleMesh1DRad
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER:: e, et INTEGER:: e, et
DO e = 1, self%numVols DO e = 1, self%numCells
!Connect Vol-Vol !Connect Cell-Cell
DO et = 1, self%numVols DO et = 1, self%numCells
IF (e /= et) THEN 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 END IF
@ -541,9 +504,9 @@ MODULE moduleMesh1DRad
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
!Connect Vol-Edge !Connect Cell-Edge
DO et = 1, self%numEdges 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 END DO
@ -553,29 +516,29 @@ MODULE moduleMesh1DRad
END SUBROUTINE connectMesh1DRad END SUBROUTINE connectMesh1DRad
SUBROUTINE connectVolVol(elemA, elemB) SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB CLASS(meshCell), INTENT(inout):: elemB
SELECT TYPE(elemA) SELECT TYPE(elemA)
TYPE IS(meshVol1DRadSegm) TYPE IS(meshCell1DRadSegm)
SELECT TYPE(elemB) SELECT TYPE(elemB)
TYPE IS(meshVol1DRadSegm) TYPE IS(meshCell1DRadSegm)
CALL connectSegmSegm(elemA, elemB) CALL connectSegmSegm(elemA, elemB)
END SELECT END SELECT
END SELECT END SELECT
END SUBROUTINE connectVolVol END SUBROUTINE connectCellCell
SUBROUTINE connectSegmSegm(elemA, elemB) SUBROUTINE connectSegmSegm(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemA
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. & IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN elemA%n2%n == elemB%n1%n) THEN
@ -594,14 +557,14 @@ MODULE moduleMesh1DRad
END SUBROUTINE connectSegmSegm END SUBROUTINE connectSegmSegm
SUBROUTINE connectVolEdge(elemA, elemB) SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA) SELECT TYPE(elemA)
TYPE IS (meshVol1DRadSegm) TYPE IS (meshCell1DRadSegm)
SELECT TYPE(elemB) SELECT TYPE(elemB)
CLASS IS(meshEdge1DRad) CLASS IS(meshEdge1DRad)
CALL connectSegmEdge(elemA, elemB) CALL connectSegmEdge(elemA, elemB)
@ -610,12 +573,12 @@ MODULE moduleMesh1DRad
END SELECT END SELECT
END SUBROUTINE connectVolEdge END SUBROUTINE connectCellEdge
SUBROUTINE connectSegmEdge(elemA, elemB) SUBROUTINE connectSegmEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA CLASS(meshCell1DRadSegm), INTENT(inout), TARGET:: elemA
CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. & IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
@ -624,7 +587,7 @@ MODULE moduleMesh1DRad
elemA%e1 => elemB elemA%e1 => elemB
elemB%e2 => elemA elemB%e2 => elemA
!Revers the normal to point inside the domain !Rever the normal to point inside the domain
elemB%normal = - elemB%normal elemB%normal = - elemB%normal
END IF END IF

1079
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

File diff suppressed because it is too large Load diff

1117
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

File diff suppressed because it is too large Load diff

715
src/modules/mesh/3DCart/moduleMesh3DCart.f90
View file

@ -11,6 +11,7 @@ MODULE moduleMesh3DCart
!Element coordinates !Element coordinates
REAL(8):: x, y, z REAL(8):: x, y, z
CONTAINS CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode3DCart PROCEDURE, PASS:: init => initNode3DCart
PROCEDURE, PASS:: getCoordinates => getCoord3DCart PROCEDURE, PASS:: getCoordinates => getCoord3DCart
@ -23,42 +24,18 @@ MODULE moduleMesh3DCart
!Connectivity to nodes !Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
CONTAINS CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge3DCartTria PROCEDURE, PASS:: init => initEdge3DCartTria
PROCEDURE, PASS:: getNodes => getNodes3DCartTria PROCEDURE, PASS:: getNodes => getNodes3DCartTria
PROCEDURE, PASS:: intersection => intersection3DCartTria PROCEDURE, PASS:: intersection => intersection3DCartTria
PROCEDURE, PASS:: randPos => randPosEdgeTria PROCEDURE, PASS:: randPos => randPosEdgeTria
PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria !PARTICULAR PROCEDURES
PROCEDURE, NOPASS, PRIVATE:: fPsi => fPsiEdgeTria
END TYPE meshEdge3DCartTria 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 !Tetrahedron volume element
TYPE, PUBLIC, EXTENDS(meshVol3DCart):: meshVol3DCartTetra TYPE, PUBLIC, EXTENDS(meshCell):: meshCell3DCartTetra
!Element Coordinates !Element Coordinates
REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0 REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0
!Connectivity to nodes !Connectivity to nodes
@ -66,28 +43,30 @@ MODULE moduleMesh3DCart
!Connectivity to adjacent elements !Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initVolTetra3DCart !meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: randPos => randPosVolTetra PROCEDURE, PASS:: init => initCellTetra
PROCEDURE, PASS:: calcVol => volumeTetra PROCEDURE, PASS:: getNodes => getNodesTetra
PROCEDURE, PASS:: randPos => randPosCellTetra
PROCEDURE, NOPASS:: fPsi => fPsiTetra PROCEDURE, NOPASS:: fPsi => fPsiTetra
PROCEDURE, NOPASS:: dPsi => dPsiTetra PROCEDURE, NOPASS:: dPsi => dPsiTetra
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXii1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXii2
PROCEDURE, PASS:: partialDer => partialDerTetra 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:: elemK => elemKTetra
PROCEDURE, PASS:: elemF => elemFTetra PROCEDURE, PASS:: elemF => elemFTetra
PROCEDURE, NOPASS:: inside => insideTetra PROCEDURE, NOPASS:: inside => insideTetra
PROCEDURE, PASS:: gatherEF => gatherEFTetra
PROCEDURE, PASS:: gatherMF => gatherMFTetra
PROCEDURE, PASS:: getNodes => getNodesTetra
PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: phy2log => phy2logTetra
PROCEDURE, PASS:: nextElement => nextElementTetra PROCEDURE, PASS:: neighbourElement => neighbourElementTetra
!PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTetra
END TYPE meshVol3DCartTetra END TYPE meshCell3DCartTetra
CONTAINS CONTAINS
!NODE FUNCTIONS !NODE FUNCTIONS
!Inits node element !Init node element
SUBROUTINE initNode3DCart(self, n, r) SUBROUTINE initNode3DCart(self, n, r)
USE moduleSpecies USE moduleSpecies
USE moduleRefParam USE moduleRefParam
@ -123,8 +102,8 @@ MODULE moduleMesh3DCart
END FUNCTION getCoord3DCart END FUNCTION getCoord3DCart
!SURFACE FUNCTIONS !EDGE FUNCTIONS
!Inits surface element !Init surface element
SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface) SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface)
USE moduleSpecies USE moduleSpecies
USE moduleBoundary USE moduleBoundary
@ -142,6 +121,7 @@ MODULE moduleMesh3DCart
INTEGER:: s INTEGER:: s
self%n = n self%n = n
self%nNodes = SIZE(p)
self%n1 => mesh%nodes(p(1))%obj self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj self%n2 => mesh%nodes(p(2))%obj
self%n3 => mesh%nodes(p(3))%obj self%n3 => mesh%nodes(p(3))%obj
@ -177,17 +157,18 @@ MODULE moduleMesh3DCart
END SUBROUTINE initEdge3DCartTria END SUBROUTINE initEdge3DCartTria
!Get nodes from surface !Get nodes from surface
PURE FUNCTION getNodes3DCartTria(self) RESULT(n) PURE FUNCTION getNodes3DCartTria(self, nNodes) RESULT(n)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge3DCartTria), INTENT(in):: self 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/) n = (/self%n1%n, self%n2%n, self%n3%n/)
END FUNCTION getNodes3DCartTria END FUNCTION getNodes3DCartTria
!Calculate intersection between position and edge
PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r) PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r)
IMPLICIT NONE IMPLICIT NONE
@ -206,21 +187,21 @@ MODULE moduleMesh3DCart
END FUNCTION intersection3DCartTria END FUNCTION intersection3DCartTria
!Calculates a random position in the surface !Calculate a random position in the surface
FUNCTION randPosEdgeTria(self) RESULT(r) FUNCTION randPosEdgeTria(self) RESULT(r)
USE moduleRandom USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge3DCartTria), INTENT(in):: self CLASS(meshEdge3DCartTria), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
REAL(8):: xii(1:3) REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3) REAL(8):: fPsi(1:3)
xii(1) = random( 0.D0, 1.D0) Xi(1) = random( 0.D0, 1.D0)
xii(2) = random( 0.D0, 1.D0 - xii(1)) Xi(2) = random( 0.D0, 1.D0 - Xi(1))
xii(3) = 0.D0 Xi(3) = 0.D0
fPsi = self%fPsi(xii) fPsi = self%fPsi(Xi)
r = (/DOT_PRODUCT(fPsi, self%x), & r = (/DOT_PRODUCT(fPsi, self%x), &
DOT_PRODUCT(fPsi, self%y), & DOT_PRODUCT(fPsi, self%y), &
DOT_PRODUCT(fPsi, self%z)/) DOT_PRODUCT(fPsi, self%z)/)
@ -228,35 +209,38 @@ MODULE moduleMesh3DCart
END FUNCTION randPosEdgeTria END FUNCTION randPosEdgeTria
!Shape functions for triangular surface !Shape functions for triangular surface
PURE FUNCTION fPsiEdgeTria(xii) RESULT(fPsi) PURE FUNCTION fPsiEdgeTria(Xi) RESULT(fPsi)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3) REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8):: fPsi(1:3)
ALLOCATE(fPsi(1:3)) fPsi(1) = 1.D0 - Xi(1) - Xi(2)
fPsi(2) = Xi(1)
fPsi(1) = 1.D0 - xii(1) - xii(2) fPsi(3) = Xi(2)
fPsi(2) = xii(1)
fPsi(3) = xii(2)
END FUNCTION fPsiEdgeTria END FUNCTION fPsiEdgeTria
!VOLUME FUNCTIONS !VOLUME FUNCTIONS
!TETRA FUNCTIONS !TETRA FUNCTIONS
!Inits tetrahedron element !Init element
SUBROUTINE initVolTetra3DCart(self, n, p, nodes) SUBROUTINE initCellTetra(self, n, p, nodes)
USE moduleRefParam USE moduleRefParam
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(out):: self CLASS(meshCell3DCartTetra), INTENT(out):: self
INTEGER, INTENT(in):: n INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node 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 self%n = n
!Assign number of nodes of cell
self%nNodes = SIZE(p)
!Assign nodes to element
self%n1 => nodes(p(1))%obj self%n1 => nodes(p(1))%obj
self%n2 => nodes(p(2))%obj self%n2 => nodes(p(2))%obj
self%n3 => nodes(p(3))%obj self%n3 => nodes(p(3))%obj
@ -271,431 +255,338 @@ MODULE moduleMesh3DCart
self%z = (/r1(3), r2(3), r3(3), r4(3)/) self%z = (/r1(3), r2(3), r3(3), r4(3)/)
!Computes the element volume !Computes the element volume
CALL self%calcVol() CALL self%calculateVolume()
!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 OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVolTetra3DCart END SUBROUTINE initCellTetra
!Random position in volume tetrahedron !Gets node indexes from cell
FUNCTION randPosVolTetra(self) RESULT(r) 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 USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), INTENT(in):: self
REAL(8):: r(1:3) REAL(8):: r(1:3)
REAL(8):: xii(1:3) REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8):: fPsi(1:4)
xii(1) = random( 0.D0, 1.D0) Xi(1) = random( 0.D0, 1.D0)
xii(2) = random( 0.D0, 1.D0 - xii(1)) Xi(2) = random( 0.D0, 1.D0 - Xi(1))
xii(3) = random( 0.D0, 1.D0 - xii(1) - xii(2)) Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2))
ALLOCATE(fPsi(1:4)) fPsi = self%fPsi(Xi, 4)
fPsi = self%fPsi(xii)
r(1) = DOT_PRODUCT(fPsi, self%x) r = (/ DOT_PRODUCT(fPsi, self%x), &
r(2) = DOT_PRODUCT(fPsi, self%y) DOT_PRODUCT(fPsi, self%y), &
r(3) = DOT_PRODUCT(fPsi, self%z) DOT_PRODUCT(fPsi, self%z) /)
END FUNCTION randPosVolTetra END FUNCTION randPosCellTetra
!Computes the element volume !Compute element functions in point Xi
PURE SUBROUTINE volumeTetra(self) PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout):: self REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: xii(1:3) INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
self%volume = 0.D0 fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3)
xii = (/0.25D0, 0.25D0, 0.25D0/) fPsi(2) = Xi(1)
self%volume = self%detJac(xii) fPsi(3) = Xi(2)
fPsi(4) = Xi(3)
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)
END FUNCTION fPsiTetra END FUNCTION fPsiTetra
!Derivative element function at coordinates xii !Compute element derivative functions in point Xi
PURE FUNCTION dPsiTetra(xii) RESULT(dPsi) PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3) REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) 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(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /)
dPsi(2,:) = dPsiTetraXii2(xii(1), xii(3)) dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /)
dPsi(3,:) = dPsiTetraXii3(xii(1), xii(2)) dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /)
END FUNCTION dPsiTetra END FUNCTION dPsiTetra
!Derivative element function respect to xii1 !Compute the derivatives in global coordinates
PURE FUNCTION dPsiTetraXii1(xii2, xii3) RESULT(dPsiXii1) PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer)
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)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:, 1:) INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz 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
dx(2) = DOT_PRODUCT(dPsi(2,:), self%x)
dx(3) = DOT_PRODUCT(dPsi(3,:), self%x)
dy(1) = DOT_PRODUCT(dPsi(1,:), self%y) pDer(1, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%x(1:4)), &
dy(2) = DOT_PRODUCT(dPsi(2,:), self%y) DOT_PRODUCT(dPsi(2,1:4), self%x(1:4)), &
dy(3) = DOT_PRODUCT(dPsi(3,:), self%y) DOT_PRODUCT(dPsi(3,1:4), self%x(1:4)) /)
dz(1) = DOT_PRODUCT(dPsi(1,:), self%z) pDer(2, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%y(1:4)), &
dz(2) = DOT_PRODUCT(dPsi(2,:), self%z) DOT_PRODUCT(dPsi(2,1:4), self%y(1:4)), &
dz(3) = DOT_PRODUCT(dPsi(3,:), self%z) 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 IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:) INTEGER, INTENT(in):: nNodes
REAL(8):: xii(1:3) REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) 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 REAL(8):: invJ(1:3,1:3), detJ
ALLOCATE(localK(1:4,1:4))
localK = 0.D0 localK = 0.D0
xii = 0.D0 Xi = 0.D0
!TODO: One point Gauss integral. Upgrade when possible !TODO: One point Gauss integral. Upgrade when possible
xii = (/ 0.25D0, 0.25D0, 0.25D0 /) Xi = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(xii) dPsi = self%dPsi(Xi, 4)
detJ = self%detJac(xii, dPsi) pDer = self%partialDer(4, dPsi)
invJ = self%invJac(xii, dPsi) detJ = self%detJac(pDer)
fPsi = self%fPsi(xii) invJ = self%invJac(pDer)
fPsi = self%fPsi(Xi, 4)
localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
END FUNCTION elemKTetra 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 IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: source(1:) INTEGER, INTENT(in):: nNodes
REAL(8), ALLOCATABLE:: localF(:) 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):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: xii(1:3) REAL(8):: pDer(1:3, 1:3)
REAL(8):: detJ, f REAL(8):: detJ, f
ALLOCATE(localF(1:4))
localF = 0.D0 localF = 0.D0
xii = 0.D0 Xi = 0.D0
!TODO: One point Gauss integral. Upgrade when possible Xi = (/ 0.25D0, 0.25D0, 0.25D0 /)
xii = (/ 0.25D0, 0.25D0, 0.25D0 /) dPsi = self%dPsi(Xi, 4)
dPsi = self%dPsi(xii) pDer = self%partialDer(4, dPsi)
detJ = self%detJac(xii, dPsi) detJ = self%detJac(pDer)
fPsi = self%fPsi(xii) fPsi = self%fPsi(Xi, 4)
f = DOT_PRODUCT(fPsi, source) f = DOT_PRODUCT(fPsi, source)
localF = f*fPsi*1.D0*detJ localF = f*fPsi*1.D0*detJ
END FUNCTION elemFTetra END FUNCTION elemFTetra
PURE FUNCTION insideTetra(xi) RESULT(ins) !Check if Xi is inside the element
PURE FUNCTION insideTetra(Xi) RESULT(ins)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins LOGICAL:: ins
ins = xi(1) >= 0.D0 .AND. & ins = Xi(1) >= 0.D0 .AND. &
xi(2) >= 0.D0 .AND. & Xi(2) >= 0.D0 .AND. &
xi(3) >= 0.D0 .AND. & Xi(3) >= 0.D0 .AND. &
1.D0 - xi(1) - xi(2) - xi(3) >= 0.D0 1.D0 - Xi(1) - Xi(2) - Xi(3) >= 0.D0
END FUNCTION insideTetra 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 IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), 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
REAL(8), INTENT(in):: r(1:3) 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):: invJ(1:3, 1:3), detJ
REAL(8):: deltaR(1:3) 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) /) deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /)
dPsi = self%dPsi(xi) dPsi = self%dPsi(Xi, 4)
invJ = self%invJac(xi, dPsi) pDer = self%partialDer(4, dPsi)
detJ = self%detJac(xi, dPsi) invJ = self%invJac(pDer)
xi = MATMUL(invJ, deltaR)/detJ detJ = self%detJac(pDer)
Xi = MATMUL(invJ, deltaR)/detJ
END FUNCTION phy2logTetra END FUNCTION phy2logTetra
SUBROUTINE nextElementTetra(self, xi, nextElement) !Get the neighbour cell for a logical position Xi
SUBROUTINE neighbourElementTetra(self, Xi, neighbourElement)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self CLASS(meshCell3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
REAL(8):: xiArray(1:4) REAL(8):: XiArray(1:4)
INTEGER:: nextInt INTEGER:: nextInt
!TODO: Review when connectivity !TODO: Review when connectivity
xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /) XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), Xi(2), Xi(1) /)
nextInt = MINLOC(xiArray, 1) nextInt = MINLOC(XiArray, 1)
NULLIFY(nextElement) NULLIFY(neighbourElement)
SELECT CASE(nextInt) SELECT CASE(nextInt)
CASE (1) CASE (1)
nextElement => self%e1 neighbourElement => self%e1
CASE (2) CASE (2)
nextElement => self%e2 neighbourElement => self%e2
CASE (3) CASE (3)
nextElement => self%e3 neighbourElement => self%e3
CASE (4) CASE (4)
nextElement => self%e4 neighbourElement => self%e4
END SELECT END SELECT
END SUBROUTINE nextElementTetra END SUBROUTINE neighbourElementTetra
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D !Calculate volume for triangular element
!Computes element Jacobian determinant PURE SUBROUTINE volumeTetra(self)
PURE FUNCTION detJ3DCart(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCart), INTENT(in)::self CLASS(meshCell3DCartTetra), INTENT(inout):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8):: Xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:, 1:) 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):: dJ
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1:3), dy(1:3), dz(1:3)
IF (PRESENT(dPsi_in)) THEN dJ = pDer(1,1)*(pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)) &
dPsi = dPsi_in - 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))
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))
END FUNCTION detJ3DCart 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 IMPLICIT NONE
CLASS(meshVol3DCart), INTENT(in):: self REAL(8), INTENT(in):: pDer(1:3, 1:3)
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):: invJ(1:3,1:3) REAL(8):: invJ(1:3,1:3)
IF(PRESENT(dPsi_in)) THEN invJ(1,1:3) = (/ (pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)), &
dPsi=dPsi_in -(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 invJ(2,1:3) = (/ -(pDer(1,2)*pDer(3,3) - pDer(1,3)*pDer(3,2)), &
dPsi = self%dPsi(xi) (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 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)), &
CALL self%partialDer(dPsi, dx, dy, dz) (pDer(1,1)*pDer(2,2) - pDer(1,2)*pDer(2,1)) /)
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 = TRANSPOSE(invJ) invJ = TRANSPOSE(invJ)
END FUNCTION invJ3DCart 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) SUBROUTINE connectMesh3DCart(self)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER:: e, et INTEGER:: e, et
DO e = 1, self%numVols DO e = 1, self%numCells
!Connect Vol-Vol !Connect Cell-Cell
DO et = 1, self%numVols DO et = 1, self%numCells
IF (e /= et) THEN 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 END IF
@ -703,9 +594,9 @@ MODULE moduleMesh3DCart
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
!Connect Vol-Edge !Connect Cell-Edge
DO et = 1, self%numEdges 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 END DO
@ -715,6 +606,46 @@ MODULE moduleMesh3DCart
END SUBROUTINE connectMesh3DCart 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 !Checks if two sets of nodes are coincidend in any order
PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident)
IMPLICIT NONE IMPLICIT NONE
@ -741,8 +672,8 @@ MODULE moduleMesh3DCart
SUBROUTINE connectTetraTetra(elemA, elemB) SUBROUTINE connectTetraTetra(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA CLASS(meshCell3DCartTetra), INTENT(inout), TARGET:: elemA
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB CLASS(meshCell3DCartTetra), INTENT(inout), TARGET:: elemB
!Check surface 1 !Check surface 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN IF (.NOT. ASSOCIATED(elemA%e1)) THEN
@ -870,11 +801,11 @@ MODULE moduleMesh3DCart
USE moduleMath USE moduleMath
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA CLASS(meshCell3DCartTetra), INTENT(inout), TARGET:: elemA
CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
INTEGER:: nodesEdge(1:3) INTEGER:: nodesEdge(1:3)
REAL(8), DIMENSION(1:3):: vec1, vec2 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 /) nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /)
@ -889,10 +820,10 @@ MODULE moduleMesh3DCart
vec2 = (/ elemA%x(3) - elemA%x(1), & vec2 = (/ elemA%x(3) - elemA%x(1), &
elemA%y(3) - elemA%y(1), & elemA%y(3) - elemA%y(1), &
elemA%z(3) - elemA%z(1) /) elemA%z(3) - elemA%z(1) /)
normVol = crossProduct(vec1, vec2) normCell = crossProduct(vec1, vec2)
normVol = normalize(normVol) normCell = normalize(normCell)
IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN
elemA%e1 => elemB elemA%e1 => elemB
elemB%e1 => elemA elemB%e1 => elemA
@ -922,10 +853,10 @@ MODULE moduleMesh3DCart
vec2 = (/ elemA%x(4) - elemA%x(2), & vec2 = (/ elemA%x(4) - elemA%x(2), &
elemA%y(4) - elemA%y(2), & elemA%y(4) - elemA%y(2), &
elemA%z(4) - elemA%z(2) /) elemA%z(4) - elemA%z(2) /)
normVol = crossProduct(vec1, vec2) normCell = crossProduct(vec1, vec2)
normVol = normalize(normVol) normCell = normalize(normCell)
IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN
elemA%e2 => elemB elemA%e2 => elemB
elemB%e1 => elemA elemB%e1 => elemA
@ -955,10 +886,10 @@ MODULE moduleMesh3DCart
vec2 = (/ elemA%x(4) - elemA%x(1), & vec2 = (/ elemA%x(4) - elemA%x(1), &
elemA%y(4) - elemA%y(1), & elemA%y(4) - elemA%y(1), &
elemA%z(4) - elemA%z(1) /) elemA%z(4) - elemA%z(1) /)
normVol = crossProduct(vec1, vec2) normCell = crossProduct(vec1, vec2)
normVol = normalize(normVol) normCell = normalize(normCell)
IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN
elemA%e3 => elemB elemA%e3 => elemB
elemB%e1 => elemA elemB%e1 => elemA
@ -988,10 +919,10 @@ MODULE moduleMesh3DCart
vec2 = (/ elemA%x(4) - elemA%x(1), & vec2 = (/ elemA%x(4) - elemA%x(1), &
elemA%y(4) - elemA%y(1), & elemA%y(4) - elemA%y(1), &
elemA%z(4) - elemA%z(1) /) elemA%z(4) - elemA%z(1) /)
normVol = crossProduct(vec1, vec2) normCell = crossProduct(vec1, vec2)
normVol = normalize(normVol) normCell = normalize(normCell)
IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN IF (DOT_PRODUCT(elemB%normal, normCell) == -1.D0) THEN
elemA%e4 => elemB elemA%e4 => elemB
elemB%e1 => elemA elemB%e1 => elemA

11
src/modules/mesh/inout/0D/moduleMeshInput0D.f90
View file

@ -41,8 +41,8 @@ MODULE moduleMeshInput0D
self%numNodes = 1 self%numNodes = 1
ALLOCATE(self%nodes(1:1)) ALLOCATE(self%nodes(1:1))
!Allocates one volume !Allocates one volume
self%numVols = 1 self%numCells = 1
ALLOCATE(self%vols(1:1)) ALLOCATE(self%cells(1:1))
!Allocates matrix K !Allocates matrix K
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
@ -59,15 +59,14 @@ MODULE moduleMeshInput0D
CALL self%nodes(1)%obj%init(1, r) CALL self%nodes(1)%obj%init(1, r)
!Creates the volume !Creates the volume
ALLOCATE(meshVol0D:: self%vols(1)%obj) ALLOCATE(meshCell0D:: self%cells(1)%obj)
CALL self%vols(1)%obj%init(1, (/ 1/), self%nodes) CALL self%cells(1)%obj%init(1, (/ 1/), self%nodes)
END SUBROUTINE read0D END SUBROUTINE read0D
SUBROUTINE readInitial0D(sp, filename, density, velocity, temperature) SUBROUTINE readInitial0D(filename, density, velocity, temperature)
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: sp
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity

2
src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
View file

@ -57,7 +57,7 @@ MODULE moduleMeshOutput0D
END IF END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') 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) CLOSE(20)
END SUBROUTINE printColl0D END SUBROUTINE printColl0D

27
src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
View file

@ -136,7 +136,7 @@ MODULE moduleMeshInputGmsh2
!Substract the number of edges to the total number of elements !Substract the number of edges to the total number of elements
!to obtain the number of volume elements !to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges self%numCells = TotalnumElem - self%numEdges
ALLOCATE(self%edges(1:self%numEdges)) ALLOCATE(self%edges(1:self%numEdges))
numEdges = self%numEdges numEdges = self%numEdges
@ -146,13 +146,13 @@ MODULE moduleMeshInputGmsh2
END DO END DO
TYPE IS(meshCollisions) TYPE IS(meshCollisions)
self%numVols = TotalnumElem self%numCells = TotalnumElem
numEdges = 0 numEdges = 0
END SELECT END SELECT
!Allocates arrays !Allocates arrays
ALLOCATE(self%vols(1:self%numVols)) ALLOCATE(self%cells(1:self%numCells))
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
@ -232,7 +232,7 @@ MODULE moduleMeshInputGmsh2
END SELECT END SELECT
!Read and initialize volumes !Read and initialize volumes
DO e = 1, self%numVols DO e = 1, self%numCells
!Reads the volume according to the geometry !Reads the volume according to the geometry
SELECT CASE(self%dimen) SELECT CASE(self%dimen)
CASE(3) CASE(3)
@ -244,7 +244,7 @@ MODULE moduleMeshInputGmsh2
!Tetrahedron element !Tetrahedron element
ALLOCATE(p(1:4)) ALLOCATE(p(1:4))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, 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 END SELECT
@ -259,13 +259,13 @@ MODULE moduleMeshInputGmsh2
!Triangular element !Triangular element
ALLOCATE(p(1:3)) ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, 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) CASE (3)
!Quadrilateral element !Quadrilateral element
ALLOCATE(p(1:4)) ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, 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 END SELECT
@ -278,13 +278,13 @@ MODULE moduleMeshInputGmsh2
!Triangular element !Triangular element
ALLOCATE(p(1:3)) ALLOCATE(p(1:3))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, 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) CASE (3)
!Quadrilateral element !Quadrilateral element
ALLOCATE(p(1:4)) ALLOCATE(p(1:4))
READ(10,*) n, elemType, eTemp, eTemp, eTemp, 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 END SELECT
@ -296,19 +296,19 @@ MODULE moduleMeshInputGmsh2
ALLOCATE(p(1:2)) ALLOCATE(p(1:2))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, 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") CASE("Cart")
ALLOCATE(p(1:2)) ALLOCATE(p(1:2))
READ(10, *) n, elemType, eTemp, eTemp, eTemp, 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
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) DEALLOCATE(p)
END DO END DO
@ -321,10 +321,9 @@ MODULE moduleMeshInputGmsh2
END SUBROUTINE readGmsh2 END SUBROUTINE readGmsh2
!Reads the initial information from an output file for an species !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 IMPLICIT NONE
INTEGER, INTENT(in):: sp
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity

16
src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
View file

@ -181,9 +181,9 @@ MODULE moduleMeshOutputGmsh2
DO c = 1, interactionMatrix(k)%amount DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c WRITE(cString, "(I2)") c
title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numVols) CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells)
DO n=1, self%numVols DO n=1, self%numCells
WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c) WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
END DO END DO
CALL writeGmsh2FooterElementData(60) CALL writeGmsh2FooterElementData(60)
@ -211,9 +211,9 @@ MODULE moduleMeshOutputGmsh2
REAL(8):: time REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring 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 IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref time = DBLE(t)*tauMin*ti_ref
@ -231,9 +231,9 @@ MODULE moduleMeshOutputGmsh2
END DO END DO
CALL writeGmsh2FooterNodeData(20) CALL writeGmsh2FooterNodeData(20)
CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numVols) CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells)
DO e=1, self%numVols DO e=1, self%numCells
WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
END DO END DO
CALL writeGmsh2FooterElementData(20) CALL writeGmsh2FooterElementData(20)

476
src/modules/mesh/moduleMesh.f90
View file

@ -24,8 +24,10 @@ MODULE moduleMesh
!Lock indicator for scattering !Lock indicator for scattering
INTEGER(KIND=OMP_LOCK_KIND):: lock INTEGER(KIND=OMP_LOCK_KIND):: lock
CONTAINS CONTAINS
!DEFERED PROCEDURES
PROCEDURE(initNode_interface), DEFERRED, PASS:: init PROCEDURE(initNode_interface), DEFERRED, PASS:: init
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
!GENERIC PROCEDURES
PROCEDURE, PASS:: resetOutput PROCEDURE, PASS:: resetOutput
END TYPE meshNode END TYPE meshNode
@ -66,10 +68,12 @@ MODULE moduleMesh
!Parent of Edge element !Parent of Edge element
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge
!Connectivity to vols !Nomber of nodes in the edge
CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL() INTEGER:: nNodes
!Connectivity to vols in meshColl !Connectivity to cells
CLASS(meshVol), POINTER:: eColl => NULL() CLASS(meshCell), POINTER:: e1 => NULL(), e2 => NULL()
!Connectivity to cells in meshColl
CLASS(meshCell), POINTER:: eColl => NULL()
!Normal vector !Normal vector
REAL(8):: normal(1:3) REAL(8):: normal(1:3)
!Weight for random injection of particles !Weight for random injection of particles
@ -81,6 +85,7 @@ MODULE moduleMesh
!Physical surface for the edge !Physical surface for the edge
INTEGER:: physicalSurface INTEGER:: physicalSurface
CONTAINS CONTAINS
!DEFERED PROCEDURES
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
@ -102,10 +107,11 @@ MODULE moduleMesh
END SUBROUTINE initEdge_interface END SUBROUTINE initEdge_interface
!Get nodes index from node !Get nodes index from node
PURE FUNCTION getNodesEdge_interface(self) RESULT(n) PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n)
IMPORT:: meshEdge IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self CLASS(meshEdge), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
END FUNCTION getNodesEdge_interface END FUNCTION getNodesEdge_interface
@ -146,8 +152,10 @@ MODULE moduleMesh
END TYPE meshEdgeCont END TYPE meshEdgeCont
!Parent of Volume element !Parent of cell element
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshCell
!Number of nodes in the cell
INTEGER:: nNodes
!Maximum collision rate !Maximum collision rate
REAL(8), ALLOCATABLE:: sigmaVrelMax(:) REAL(8), ALLOCATABLE:: sigmaVrelMax(:)
!Arrays for counting number of collisions !Arrays for counting number of collisions
@ -161,114 +169,160 @@ MODULE moduleMesh
!Total weight of particles inside cell !Total weight of particles inside cell
REAL(8), ALLOCATABLE:: totalWeight(:) REAL(8), ALLOCATABLE:: totalWeight(:)
CONTAINS CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init !DEFERRED PROCEDURES
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes !Init the cell
PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos 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(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE, PASS:: scatter PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
PROCEDURE(gatherMF_interface), DEFERRED, PASS:: gatherMF 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(elemK_interface), DEFERRED, PASS:: elemK
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
PROCEDURE, PASS:: findCell !Check if particle is inside the cell
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside
PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement !Convert physical coordinates (r) into logical coordinates (Xi)
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
!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 ABSTRACT INTERFACE
SUBROUTINE initVol_interface(self, n, p, nodes) SUBROUTINE initCell_interface(self, n, p, nodes)
IMPORT:: meshVol IMPORT:: meshCell
IMPORT meshNodeCont IMPORT meshNodeCont
CLASS(meshVol), INTENT(out):: self CLASS(meshCell), INTENT(out):: self
INTEGER, INTENT(in):: n INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
END SUBROUTINE initVol_interface END SUBROUTINE initCell_interface
PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF) PURE FUNCTION getNodesCell_interface(self, nNodes) RESULT(n)
IMPORT:: meshVol IMPORT:: meshCell
CLASS(meshVol), INTENT(in):: self CLASS(meshCell), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) INTEGER, INTENT(in):: nNodes
REAL(8):: EF(1:3) INTEGER:: n(1:nNodes)
END FUNCTION gatherEF_interface END FUNCTION getNodesCell_interface
PURE FUNCTION gatherMF_interface(self, xi) RESULT(MF) FUNCTION randPosCell_interface(self) RESULT(r)
IMPORT:: meshVol IMPORT:: meshCell
CLASS(meshVol), INTENT(in):: self CLASS(meshCell), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3) REAL(8):: r(1:3)
REAL(8):: MF(1:3)
END FUNCTION gatherMF_interface END FUNCTION randPosCell_interface
PURE FUNCTION getNodesVol_interface(self) RESULT(n) PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi)
IMPORT:: meshVol REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshVol), INTENT(in):: self INTEGER, INTENT(in):: nNodes
INTEGER, ALLOCATABLE:: n(:) REAL(8):: fPsi(1:nNodes)
END FUNCTION getNodesVol_interface
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface END FUNCTION fPsi_interface
PURE FUNCTION elemK_interface(self) RESULT(localK) PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi)
IMPORT:: meshVol REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshVol), INTENT(in):: self INTEGER, INTENT(in):: nNodes
REAL(8), ALLOCATABLE:: localK(:,:) 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 END FUNCTION elemK_interface
PURE FUNCTION elemF_interface(self, source) RESULT(localF) PURE FUNCTION elemF_interface(self, nNodes, source) RESULT(localF)
IMPORT:: meshVol IMPORT:: meshCell
CLASS(meshVol), INTENT(in):: self CLASS(meshCell), INTENT(in):: self
REAL(8), INTENT(in):: source(1:) INTEGER, INTENT(in):: nNodes
REAL(8), ALLOCATABLE:: localF(:) REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
END FUNCTION elemF_interface END FUNCTION elemF_interface
SUBROUTINE nextElement_interface(self, xi, nextElement) PURE FUNCTION inside_interface(Xi) RESULT(ins)
IMPORT:: meshVol, meshElement IMPORT:: meshCell
CLASS(meshVol), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3)
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)
LOGICAL:: ins LOGICAL:: ins
END FUNCTION inside_interface END FUNCTION inside_interface
FUNCTION randPosVol_interface(self) RESULT(r) PURE FUNCTION phy2log_interface(self,r) RESULT(Xi)
IMPORT:: meshVol IMPORT:: meshCell
CLASS(meshVol), INTENT(in):: self CLASS(meshCell), INTENT(in):: self
REAL(8):: r(1:3) 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 END INTERFACE
!Containers for volumes in the mesh !Containers for cells in the mesh
TYPE:: meshVolCont TYPE:: meshCellCont
CLASS(meshVol), ALLOCATABLE:: obj CLASS(meshCell), ALLOCATABLE:: obj
END TYPE meshVolCont END TYPE meshCellCont
!Generic mesh type !Generic mesh type
TYPE, ABSTRACT:: meshGeneric TYPE, ABSTRACT:: meshGeneric
@ -277,16 +331,18 @@ MODULE moduleMesh
!Geometry of the mesh !Geometry of the mesh
CHARACTER(:), ALLOCATABLE:: geometry CHARACTER(:), ALLOCATABLE:: geometry
!Number of elements !Number of elements
INTEGER:: numNodes, numVols INTEGER:: numNodes, numCells
!Array of nodes !Array of nodes
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of volume elements !Array of cell elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:) TYPE(meshCellCont), ALLOCATABLE:: cells(:)
!PROCEDURES SPECIFIC OF FILE TYPE
PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL()
PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL()
PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
CONTAINS CONTAINS
!GENERIC PROCEDURES
PROCEDURE, PASS:: doCollisions PROCEDURE, PASS:: doCollisions
END TYPE meshGeneric END TYPE meshGeneric
@ -295,14 +351,12 @@ MODULE moduleMesh
!Reads the mesh from a file !Reads the mesh from a file
SUBROUTINE readMesh_interface(self, filename) SUBROUTINE readMesh_interface(self, filename)
IMPORT meshGeneric IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
END SUBROUTINE readMesh_interface END SUBROUTINE readMesh_interface
SUBROUTINE readInitial_interface(sp, filename, density, velocity, temperature) SUBROUTINE readInitial_interface(filename, density, velocity, temperature)
INTEGER, INTENT(in):: sp
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
@ -310,18 +364,16 @@ MODULE moduleMesh
END SUBROUTINE readInitial_interface END SUBROUTINE readInitial_interface
!Connects volume and edges to the mesh !Connects cell and edges to the mesh
SUBROUTINE connectMesh_interface(self) SUBROUTINE connectMesh_interface(self)
IMPORT meshGeneric IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
END SUBROUTINE connectMesh_interface END SUBROUTINE connectMesh_interface
!Prints number of collisions in each volume !Prints number of collisions in each cell
SUBROUTINE printColl_interface(self, t) SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
@ -338,28 +390,21 @@ MODULE moduleMesh
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Permutation matrix for P L U factorization !Permutation matrix for P L U factorization
INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
!PROCEDURES SPECIFIC OF FILE TYPE
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL()
PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
CONTAINS CONTAINS
!GENERIC PROCEDURES
PROCEDURE, PASS:: constructGlobalK PROCEDURE, PASS:: constructGlobalK
END TYPE meshParticles END TYPE meshParticles
ABSTRACT INTERFACE ABSTRACT INTERFACE
!Perform Coulomb Scattering
SUBROUTINE doCoulomb_interface(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(inout):: self
END SUBROUTINE doCoulomb_interface
!Prints Species data !Prints Species data
SUBROUTINE printOutput_interface(self, t) SUBROUTINE printOutput_interface(self, t)
IMPORT meshParticles IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
@ -368,21 +413,25 @@ MODULE moduleMesh
!Prints EM info !Prints EM info
SUBROUTINE printEM_interface(self, t) SUBROUTINE printEM_interface(self, t)
IMPORT meshParticles IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface 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 !Prints average values
SUBROUTINE printAverage_interface(self) SUBROUTINE printAverage_interface(self)
IMPORT meshParticles IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
END SUBROUTINE printAverage_interface END SUBROUTINE printAverage_interface
END INTERFACE END INTERFACE
TYPE(meshParticles), TARGET:: mesh TYPE(meshParticles), TARGET:: mesh
@ -390,6 +439,7 @@ MODULE moduleMesh
!Collision (MCC) mesh !Collision (MCC) mesh
TYPE, EXTENDS(meshGeneric):: meshCollisions TYPE, EXTENDS(meshGeneric):: meshCollisions
CONTAINS CONTAINS
!GENERIC PROCEDURES
END TYPE meshCollisions END TYPE meshCollisions
@ -398,7 +448,6 @@ MODULE moduleMesh
ABSTRACT INTERFACE ABSTRACT INTERFACE
SUBROUTINE readMeshColl_interface(self, filename) SUBROUTINE readMeshColl_interface(self, filename)
IMPORT meshCollisions IMPORT meshCollisions
CLASS(meshCollisions), INTENT(inout):: self CLASS(meshCollisions), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
@ -416,7 +465,7 @@ MODULE moduleMesh
!Pointer to mesh used for MC collisions !Pointer to mesh used for MC collisions
CLASS(meshGeneric), POINTER:: meshForMCC => NULL() 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() PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL()
ABSTRACT INTERFACE ABSTRACT INTERFACE
@ -431,24 +480,29 @@ MODULE moduleMesh
!Logical to indicate if an specific mesh for MC Collisions is used !Logical to indicate if an specific mesh for MC Collisions is used
LOGICAL:: doubleMesh LOGICAL:: doubleMesh
!Logical to indicate if MCC collisions are performed
LOGICAL:: doMCC
!Complete path for the two meshes !Complete path for the two meshes
CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle
CONTAINS CONTAINS
!Constructs the global K matrix !Constructs the global K matrix
SUBROUTINE constructGlobalK(self) PURE SUBROUTINE constructGlobalK(self)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(inout):: self CLASS(meshParticles), INTENT(inout):: self
INTEGER:: e INTEGER:: e
INTEGER:: nNodes
INTEGER, ALLOCATABLE:: n(:) INTEGER, ALLOCATABLE:: n(:)
REAL(8), ALLOCATABLE:: localK(:,:) REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j INTEGER:: i, j
DO e = 1, self%numVols DO e = 1, self%numCells
n = self%vols(e)%obj%getNodes() nNodes = self%cells(e)%obj%nNodes
localK = self%vols(e)%obj%elemK() ALLOCATE(n(1:nNodes))
nNodes = SIZE(n) 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 i = 1, nNodes
DO j = 1, nNodes DO j = 1, nNodes
@ -458,6 +512,8 @@ MODULE moduleMesh
END DO END DO
DEALLOCATE(n, localK)
END DO END DO
END SUBROUTINE constructGlobalK END SUBROUTINE constructGlobalK
@ -480,30 +536,89 @@ MODULE moduleMesh
END SUBROUTINE resetOutput END SUBROUTINE resetOutput
!Scatters particle properties into vol nodes !Gather the value of valNodes (scalar) at position Xi
SUBROUTINE scatter(self, part) 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 moduleMath
USE moduleSpecies USE moduleSpecies
USE OMP_LIB USE OMP_LIB
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self CLASS(meshCell), INTENT(inout):: self
INTEGER, INTENT(in):: nNodes
CLASS(particle), INTENT(in):: part CLASS(particle), INTENT(in):: part
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8):: fPsi(1:nNodes)
INTEGER, ALLOCATABLE:: volNodes(:) INTEGER:: cellNodes(1:nNodes)
REAL(8):: tensorS(1:3, 1:3) REAL(8):: tensorS(1:3, 1:3)
INTEGER:: sp INTEGER:: sp
INTEGER:: i, nNodes INTEGER:: i
CLASS(meshNode), POINTER:: node 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) tensorS = outerProduct(part%v, part%v)
sp = part%species%n sp = part%species%n
volNodes = self%getNodes()
nNodes = SIZE(volNodes)
DO i = 1, nNodes DO i = 1, nNodes
node => mesh%nodes(volNodes(i))%obj node => mesh%nodes(cellNodes(i))%obj
CALL OMP_SET_LOCK(node%lock) CALL OMP_SET_LOCK(node%lock)
node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) 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(:) node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:)
@ -521,38 +636,41 @@ MODULE moduleMesh
USE OMP_LIB USE OMP_LIB
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self CLASS(meshCell), INTENT(inout):: self
CLASS(particle), INTENT(inout), TARGET:: part CLASS(particle), INTENT(inout), TARGET:: part
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell
REAL(8):: xi(1:3) REAL(8):: Xi(1:3)
CLASS(meshElement), POINTER:: nextElement CLASS(meshElement), POINTER:: neighbourElement
INTEGER:: sp INTEGER:: sp
xi = self%phy2log(part%r) Xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell !Checks if particle is inside 'self' cell
IF (self%inside(xi)) THEN IF (self%inside(Xi)) THEN
part%vol = self%n part%cell = self%n
part%xi = xi part%Xi = Xi
part%n_in = .TRUE. part%n_in = .TRUE.
!Assign particle to listPart_in !Assign particle to listPart_in
IF (doMCC) THEN
CALL OMP_SET_LOCK(self%lock) CALL OMP_SET_LOCK(self%lock)
sp = part%species%n sp = part%species%n
CALL self%listPart_in(sp)%add(part) CALL self%listPart_in(sp)%add(part)
self%totalWeight(sp) = self%totalWeight(sp) + part%weight self%totalWeight(sp) = self%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(self%lock) CALL OMP_UNSET_LOCK(self%lock)
END IF
ELSE ELSE
!If not, searches for a neighbour and repeats the process. !If not, searches for a neighbour and repeats the process.
CALL self%nextElement(xi, nextElement) CALL self%neighbourElement(Xi, neighbourElement)
!Defines the next step !Defines the next step
SELECT TYPE(nextElement) SELECT TYPE(neighbourElement)
CLASS IS(meshVol) CLASS IS(meshCell)
!Particle moved to new cell, repeat find procedure !Particle moved to new cell, repeat find procedure
CALL nextElement%findCell(part, self) CALL neighbourElement%findCell(part, self)
CLASS IS (meshEdge) CLASS IS (meshEdge)
!Particle encountered a surface, apply boundary !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 particle is still inside the domain, call findCell
IF (part%n_in) THEN IF (part%n_in) THEN
@ -575,14 +693,14 @@ MODULE moduleMesh
END SUBROUTINE findCell 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) SUBROUTINE findCellSameMesh(part)
USE moduleSpecies USE moduleSpecies
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
part%volColl = part%vol part%cellColl = part%cell
END SUBROUTINE findCellSameMesh END SUBROUTINE findCellSameMesh
@ -595,31 +713,34 @@ MODULE moduleMesh
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
LOGICAL:: found LOGICAL:: found
CLASS(meshVol), POINTER:: vol CLASS(meshCell), POINTER:: cell
REAL(8), DIMENSION(1:3):: xii REAL(8), DIMENSION(1:3):: Xi
CLASS(meshElement), POINTER:: nextElement CLASS(meshElement), POINTER:: neighbourElement
INTEGER:: sp INTEGER:: sp
found = .FALSE. found = .FALSE.
vol => meshColl%vols(part%volColl)%obj cell => meshColl%cells(part%cellColl)%obj
DO WHILE(.NOT. found) DO WHILE(.NOT. found)
xii = vol%phy2log(part%r) Xi = cell%phy2log(part%r)
IF (vol%inside(xii)) THEN IF (cell%inside(Xi)) THEN
part%volColl = vol%n part%cellColl = cell%n
CALL OMP_SET_LOCK(vol%lock) IF (doMCC) THEN
CALL OMP_SET_LOCK(cell%lock)
sp = part%species%n sp = part%species%n
CALL vol%listPart_in(sp)%add(part) CALL cell%listPart_in(sp)%add(part)
vol%totalWeight(sp) = vol%totalWeight(sp) + part%weight cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(vol%lock) CALL OMP_UNSET_LOCK(cell%lock)
END IF
found = .TRUE. found = .TRUE.
ELSE ELSE
CALL vol%nextElement(xii, nextElement) CALL cell%neighbourElement(Xi, neighbourElement)
SELECT TYPE(nextElement) SELECT TYPE(neighbourElement)
CLASS IS(meshVol) CLASS IS(meshCell)
!Try next element !Try next element
vol => nextElement cell => neighbourElement
CLASS DEFAULT CLASS DEFAULT
!Should never happend, but just in case, stops loops !Should never happend, but just in case, stops loops
@ -644,15 +765,15 @@ MODULE moduleMesh
REAL(8), DIMENSION(1:3), INTENT(in):: r REAL(8), DIMENSION(1:3), INTENT(in):: r
INTEGER:: nVol INTEGER:: nVol
INTEGER:: e INTEGER:: e
REAL(8), DIMENSION(1:3):: xii REAL(8), DIMENSION(1:3):: Xi
!Inits RESULT !Inits RESULT
nVol = 0 nVol = 0
DO e = 1, self%numVols DO e = 1, self%numCells
xii = self%vols(e)%obj%phy2log(r) Xi = self%cells(e)%obj%phy2log(r)
IF(self%vols(e)%obj%inside(xii)) THEN IF(self%cells(e)%obj%inside(Xi)) THEN
nVol = self%vols(e)%obj%n nVol = self%cells(e)%obj%n
EXIT EXIT
END IF END IF
@ -675,7 +796,7 @@ MODULE moduleMesh
CLASS(meshGeneric), INTENT(inout), TARGET:: self CLASS(meshGeneric), INTENT(inout), TARGET:: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
INTEGER:: e INTEGER:: e
CLASS(meshVol), POINTER:: vol CLASS(meshCell), POINTER:: cell
INTEGER:: k, i, j INTEGER:: k, i, j
INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision REAL(8):: pMax !Maximum probability of collision
@ -686,21 +807,22 @@ MODULE moduleMesh
REAL(8):: vRel, rMass, eRel REAL(8):: vRel, rMass, eRel
REAL(8):: sigmaVrelTotal REAL(8):: sigmaVrelTotal
REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) 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 IF (MOD(t, everyColl) == 0) THEN
!Collisions need to be performed in this iteration !Collisions need to be performed in this iteration
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) !$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 !TODO: Simplify this, to many sublevels
!Iterate over the number of pairs !Iterate over the number of pairs
DO k = 1, nCollPairs DO k = 1, nCollPairs
!Reset tally of collisions !Reset tally of collisions
IF (collOutput) THEN IF (collOutput) THEN
vol%tallyColl(k)%tally = 0 cell%tallyColl(k)%tally = 0
END IF END IF
@ -710,8 +832,8 @@ MODULE moduleMesh
j = interactionMatrix(k)%sp_j%n j = interactionMatrix(k)%sp_j%n
!Number of particles per species in the collision pair !Number of particles per species in the collision pair
nPart_i = vol%listPart_in(i)%amount nPart_i = cell%listPart_in(i)%amount
nPart_j = vol%listPart_in(j)%amount nPart_j = cell%listPart_in(j)%amount
IF (nPart_i > 0 .AND. nPart_j > 0) THEN IF (nPart_i > 0 .AND. nPart_j > 0) THEN
!Total number of particles for the collision pair !Total number of particles for the collision pair
@ -721,15 +843,15 @@ MODULE moduleMesh
nColl = 0 nColl = 0
!Probability of collision for pair i-j !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 !Number of collisions in the cell
nColl = NINT(REAL(nPart)*pMax*0.5D0) nColl = NINT(REAL(nPart)*pMax*0.5D0)
!Converts the list of particles to an array for easy access !Converts the list of particles to an array for easy access
IF (nColl > 0) THEN IF (nColl > 0) THEN
partTemp_i = vol%listPart_in(i)%convert2Array() partTemp_i = cell%listPart_in(i)%convert2Array()
partTemp_j = vol%listPart_in(j)%convert2Array() partTemp_j = cell%listPart_in(j)%convert2Array()
END IF END IF
@ -737,10 +859,10 @@ MODULE moduleMesh
!Select random particles !Select random particles
part_i => NULL() part_i => NULL()
part_j => NULL() part_j => NULL()
rnd = random(1, nPart_i) rnd_int = random(1, nPart_i)
part_i => partTemp_i(rnd)%part part_i => partTemp_i(rnd_int)%part
rnd = random(1, nPart_j) rnd_int = random(1, nPart_j)
part_j => partTemp_j(rnd)%part part_j => partTemp_j(rnd_int)%part
!If they are the same particle, skip !If they are the same particle, skip
!TODO: Maybe try to improve this !TODO: Maybe try to improve this
IF (ASSOCIATED(part_i, part_j)) THEN IF (ASSOCIATED(part_i, part_j)) THEN
@ -764,32 +886,32 @@ MODULE moduleMesh
CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
!Update maximum sigma*v_rel !Update maximum sigma*v_rel
IF (sigmaVrelTotal > vol%sigmaVrelMax(k)) THEN IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN
vol%sigmaVrelMax(k) = sigmaVrelTotal cell%sigmaVrelMax(k) = sigmaVrelTotal
END IF END IF
ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
probabilityColl = 0.0 probabilityColl = 0.0
DO c = 1, interactionMatrix(k)%amount 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 END DO
!Selects random number between 0 and 1 !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 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 !Loop over collisions
DO c = 1, interactionMatrix(k)%amount 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) CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
!If collisions are gonna be output, count the collision !If collisions are gonna be output, count the collision
IF (collOutput) THEN 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 END IF

14
src/modules/mesh/moduleMeshBoundary.f90
View file

@ -1,4 +1,4 @@
!moduleMeshBoundary: Boundary functions !moduleMeshBoundary: Boundary functions for the mesh edges
MODULE moduleMeshBoundary MODULE moduleMeshBoundary
USE moduleMesh USE moduleMesh
@ -55,10 +55,10 @@ MODULE moduleMeshBoundary
!Scatter particle in associated volume !Scatter particle in associated volume
IF (ASSOCIATED(edge%e1)) THEN IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part) CALL edge%e1%scatter(edge%e1%nNodes, part)
ELSE ELSE
CALL edge%e2%scatter(part) CALL edge%e2%scatter(edge%e2%nNodes, part)
END IF END IF
@ -156,11 +156,11 @@ MODULE moduleMeshBoundary
newElectron%r = edge%randPos() newElectron%r = edge%randPos()
newIon%r = newElectron%r newIon%r = newElectron%r
newElectron%vol = part%vol newElectron%cell = part%cell
newIon%vol = part%vol newIon%cell = part%cell
newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r)
newIon%xi = newElectron%xi newIon%Xi = newElectron%Xi
newElectron%weight = part%weight newElectron%weight = part%weight
newIon%weight = newElectron%weight newIon%weight = newElectron%weight

9
src/modules/moduleCollisions.f90
View file

@ -111,13 +111,13 @@ MODULE moduleCollisions
IMPLICIT NONE IMPLICIT NONE
REAL(8):: n(1:3) REAL(8):: n(1:3)
REAL(8):: cosXii, sinXii, eps REAL(8):: cosXi, sinXi, eps
cosXii = random(-1.D0, 1.D0) cosXi = random(-1.D0, 1.D0)
sinXii = DSQRT(1.D0 - cosXii**2) sinXi = DSQRT(1.D0 - cosXi**2)
eps = random(0.D0, PI2) eps = random(0.D0, PI2)
n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /) n = (/ cosXi, sinXi*DCOS(eps), sinXi*DSIN(eps) /)
END FUNCTION randomDirectionVHS END FUNCTION randomDirectionVHS
@ -439,7 +439,6 @@ MODULE moduleCollisions
REAL(8), INTENT(in):: vRel REAL(8), INTENT(in):: vRel
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
TYPE(particle), POINTER:: electron => NULL(), ion => NULL() TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
REAL(8):: sigmaVrel
REAL(8), DIMENSION(1:3):: vp_i REAL(8), DIMENSION(1:3):: vp_i
TYPE(particle), POINTER:: remainingIon => NULL() TYPE(particle), POINTER:: remainingIon => NULL()

10
src/modules/moduleInject.f90
View file

@ -132,7 +132,7 @@ MODULE moduleInject
IF (doubleMesh) THEN IF (doubleMesh) THEN
nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
IF (nVolColl > 0) THEN IF (nVolColl > 0) THEN
mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj
ELSE ELSE
CALL criticalError("No connection between edge and meshColl", "initInject") CALL criticalError("No connection between edge and meshColl", "initInject")
@ -285,16 +285,16 @@ MODULE moduleInject
partInj(n)%r = randomEdge%randPos() partInj(n)%r = randomEdge%randPos()
!Volume associated to the edge: !Volume associated to the edge:
IF (ASSOCIATED(randomEdge%e1)) THEN IF (ASSOCIATED(randomEdge%e1)) THEN
partInj(n)%vol = randomEdge%e1%n partInj(n)%cell = randomEdge%e1%n
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
partInj(n)%vol = randomEdge%e2%n partInj(n)%cell = randomEdge%e2%n
ELSE ELSE
CALL criticalError("No Volume associated to edge", 'addParticles') CALL criticalError("No Volume associated to edge", 'addParticles')
END IF END IF
partInj(n)%volColl = randomEdge%eColl%n partInj(n)%cellColl = randomEdge%eColl%n
sp = self%species%n sp = self%species%n
!Assign particle type !Assign particle type
@ -305,7 +305,7 @@ MODULE moduleInject
self%v(3)%obj%randomVel() /) self%v(3)%obj%randomVel() /)
!Obtain natural coordinates of particle in cell !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 !Push new particle with the minimum time step
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle !Assign cell to new particle

3
src/modules/moduleList.f90
View file

@ -92,9 +92,12 @@ MODULE moduleList
DEALLOCATE(current) DEALLOCATE(current)
current => next current => next
END DO END DO
IF (ASSOCIATED(self%head)) NULLIFY(self%head) IF (ASSOCIATED(self%head)) NULLIFY(self%head)
IF (ASSOCIATED(self%tail)) NULLIFY(self%tail) IF (ASSOCIATED(self%tail)) NULLIFY(self%tail)
self%amount = 0 self%amount = 0
END SUBROUTINE eraseList END SUBROUTINE eraseList
SUBROUTINE setLock(self) SUBROUTINE setLock(self)

6
src/modules/moduleSpecies.f90
View file

@ -38,9 +38,9 @@ MODULE moduleSpecies
REAL(8):: r(1:3) !Position REAL(8):: r(1:3) !Position
REAL(8):: v(1:3) !Velocity REAL(8):: v(1:3) !Velocity
CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
INTEGER:: vol !Index of element in which the particle is located INTEGER:: cell !Index of element in which the particle is located
INTEGER:: volColl !Index of element in which the particle is located in the Collision Mesh 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. 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 LOGICAL:: n_in !Flag that indicates if a particle is in the domain
REAL(8):: weight=0.D0 !weight of particle REAL(8):: weight=0.D0 !weight of particle

48
src/modules/solver/electromagnetic/moduleEM.f90
View file

@ -26,12 +26,14 @@ MODULE moduleEM
CLASS(boundaryEM), INTENT(in):: self CLASS(boundaryEM), INTENT(in):: self
CLASS(meshEdge):: edge CLASS(meshEdge):: edge
INTEGER:: nNodes
INTEGER, ALLOCATABLE:: nodes(:) INTEGER, ALLOCATABLE:: nodes(:)
INTEGER:: n 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) SELECT CASE(self%typeEM)
CASE ("dirichlet") CASE ("dirichlet")
mesh%K(nodes(n), :) = 0.D0 mesh%K(nodes(n), :) = 0.D0
@ -46,40 +48,6 @@ MODULE moduleEM
END SUBROUTINE 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 !Assemble the source vector based on the charge density to solve Poisson's equation
SUBROUTINE assembleSourceVector(vectorF) SUBROUTINE assembleSourceVector(vectorF)
USE moduleMesh USE moduleMesh
@ -99,9 +67,9 @@ MODULE moduleEM
!$OMP END SINGLE !$OMP END SINGLE
!$OMP DO REDUCTION(+:vectorF) !$OMP DO REDUCTION(+:vectorF)
DO e = 1, mesh%numVols DO e = 1, mesh%numCells
nodes = mesh%vols(e)%obj%getNodes() nNodes = mesh%cells(e)%obj%nNodes
nNodes = SIZE(nodes) nodes = mesh%cells(e)%obj%getNodes(nNodes)
!Calculates charge density (rho) in element nodes !Calculates charge density (rho) in element nodes
ALLOCATE(rho(1:nNodes)) ALLOCATE(rho(1:nNodes))
rho = 0.D0 rho = 0.D0
@ -113,7 +81,7 @@ MODULE moduleEM
END DO END DO
!Calculates local F vector !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 !Assign local F to global F
DO i = 1, nNodes DO i = 1, nNodes

72
src/modules/solver/moduleSolver.f90
View file

@ -43,14 +43,14 @@ MODULE moduleSolver
END SUBROUTINE solveEM_interface END SUBROUTINE solveEM_interface
!Apply nonAnalogue scheme to a particle !Apply nonAnalogue scheme to a particle
SUBROUTINE weightingScheme_interface(part, volOld, volNew) SUBROUTINE weightingScheme_interface(part, cellOld, cellNew)
USE moduleSpecies USE moduleSpecies
USE moduleMesh USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER, INTENT(in):: volOld CLASS(meshCell), POINTER, INTENT(in):: cellOld
CLASS(meshVol), POINTER, INTENT(inout):: volNew CLASS(meshCell), POINTER, INTENT(inout):: cellNew
END SUBROUTINE weightingScheme_interface END SUBROUTINE weightingScheme_interface
@ -204,7 +204,10 @@ MODULE moduleSolver
DO n = 1, partList%amount DO n = 1, partList%amount
partNext => partCurr%next partNext => partCurr%next
partArray(nStart + n) = partCurr%part partArray(nStart + n) = partCurr%part
IF (nProbes > 0) THEN
CALL doProbes(partArray(nStart+n)) CALL doProbes(partArray(nStart+n))
END IF
DEALLOCATE(partCurr) DEALLOCATE(partCurr)
partCurr => partNext partCurr => partNext
@ -270,10 +273,13 @@ MODULE moduleSolver
IF (partInj(n)%n_in) THEN IF (partInj(n)%n_in) THEN
nn = nn + 1 nn = nn + 1
partOld(nn) = partInj(n) partOld(nn) = partInj(n)
IF (nProbes > 0) THEN
CALL doProbes(partOld(nn)) CALL doProbes(partOld(nn))
END IF END IF
END IF
END DO END DO
!$OMP SECTION !$OMP SECTION
@ -283,10 +289,13 @@ MODULE moduleSolver
IF (partTemp(n)%n_in) THEN IF (partTemp(n)%n_in) THEN
nn = nn + 1 nn = nn + 1
partOld(nn) = partTemp(n) partOld(nn) = partTemp(n)
IF (nProbes > 0) THEN
CALL doProbes(partOld(nn)) CALL doProbes(partOld(nn))
END IF END IF
END IF
END DO END DO
!$OMP SECTION !$OMP SECTION
@ -313,11 +322,12 @@ MODULE moduleSolver
!$OMP SECTION !$OMP SECTION
!Erase the list of particles inside the cell if particles have been pushed !Erase the list of particles inside the cell if particles have been pushed
IF (doMCC) THEN
DO s = 1, nSpecies DO s = 1, nSpecies
DO e = 1, mesh%numVols DO e = 1, mesh%numCells
IF (solver%pusher(s)%pushSpecies) THEN IF (solver%pusher(s)%pushSpecies) THEN
CALL mesh%vols(e)%obj%listPart_in(s)%erase() CALL mesh%cells(e)%obj%listPart_in(s)%erase()
mesh%vols(e)%obj%totalWeight(s) = 0.D0 mesh%cells(e)%obj%totalWeight(s) = 0.D0
END IF END IF
@ -325,13 +335,16 @@ MODULE moduleSolver
END DO END DO
END IF
!$OMP SECTION !$OMP SECTION
!Erase the list of particles inside the cell in coll mesh !Erase the list of particles inside the cell in coll mesh
IF (doubleMesh) THEN
DO s = 1, nSpecies DO s = 1, nSpecies
DO e = 1, meshColl%numVols DO e = 1, meshColl%numCells
IF (solver%pusher(s)%pushSpecies) THEN IF (solver%pusher(s)%pushSpecies) THEN
CALL meshColl%vols(e)%obj%listPart_in(s)%erase() CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
meshColl%vols(e)%obj%totalWeight(s) = 0.D0 meshColl%cells(e)%obj%totalWeight(s) = 0.D0
END IF END IF
@ -339,6 +352,8 @@ MODULE moduleSolver
END DO END DO
END IF
!$OMP END SECTIONS !$OMP END SECTIONS
!$OMP SINGLE !$OMP SINGLE
@ -354,11 +369,13 @@ MODULE moduleSolver
IMPLICIT NONE IMPLICIT NONE
INTEGER:: n INTEGER:: n
CLASS(meshCell), POINTER:: cell
!Loops over the particles to scatter them !Loops over the particles to scatter them
!$OMP DO !$OMP DO
DO n = 1, nPartOld 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 END DO
!$OMP END DO !$OMP END DO
@ -376,28 +393,28 @@ MODULE moduleSolver
END SUBROUTINE doEMField END SUBROUTINE doEMField
!Split particles as a function of cell volume and splits particle !Split particles as a function of cell volume and splits particle
SUBROUTINE volumeWScheme(part, volOld, volNew) SUBROUTINE volumeWScheme(part, cellOld, cellNew)
USE moduleSpecies USE moduleSpecies
USE moduleMesh USE moduleMesh
USE moduleRandom USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER, INTENT(in):: volOld CLASS(meshCell), POINTER, INTENT(in):: cellOld
CLASS(meshVol), POINTER, INTENT(inout):: volNew CLASS(meshCell), POINTER, INTENT(inout):: cellNew
REAL(8):: fractionVolume, pSplit REAL(8):: fractionVolume, pSplit
!If particle changes volume to smaller cell !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 part%weight >= part%species%weight*1.0D-1) THEN
fractionVolume = volOld%volume/volNew%volume fractionVolume = cellOld%volume/cellNew%volume
!Calculate probability of splitting particle !Calculate probability of splitting particle
pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1) pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1)
IF (random() < pSplit) THEN IF (random() < pSplit) THEN
!Split particle in two !Split particle in two
CALL splitParticle(part, 2, volNew) CALL splitParticle(part, 2, cellNew)
END IF END IF
@ -407,7 +424,7 @@ MODULE moduleSolver
!Subroutine to split the particle 'part' into a number 'nSplit' of particles. !Subroutine to split the particle 'part' into a number 'nSplit' of particles.
!'nSplit-1' particles are added to the partNAScheme list !'nSplit-1' particles are added to the partNAScheme list
SUBROUTINE splitParticle(part, nSplit, vol) SUBROUTINE splitParticle(part, nSplit, cell)
USE moduleSpecies USE moduleSpecies
USE moduleList USE moduleList
USE moduleMesh USE moduleMesh
@ -416,7 +433,7 @@ MODULE moduleSolver
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
INTEGER, INTENT(in):: nSplit INTEGER, INTENT(in):: nSplit
CLASS(meshVol), INTENT(inout):: vol CLASS(meshCell), INTENT(inout):: cell
REAL(8):: newWeight REAL(8):: newWeight
TYPE(particle), POINTER:: newPart TYPE(particle), POINTER:: newPart
INTEGER:: p INTEGER:: p
@ -438,10 +455,13 @@ MODULE moduleSolver
CALL partWScheme%add(newPart) CALL partWScheme%add(newPart)
CALL partWScheme%unsetLock() CALL partWScheme%unsetLock()
!Add particle to cell list !Add particle to cell list
CALL OMP_SET_lock(vol%lock)
sp = part%species%n sp = part%species%n
CALL vol%listPart_in(sp)%add(newPart) IF (doMCC) THEN
CALL OMP_UNSET_lock(vol%lock) CALL OMP_SET_lock(cell%lock)
CALL cell%listPart_in(sp)%add(newPart)
CALL OMP_UNSET_lock(cell%lock)
END IF
END DO END DO
@ -454,18 +474,18 @@ MODULE moduleSolver
CLASS(solverGeneric), INTENT(in):: self CLASS(solverGeneric), INTENT(in):: self
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER:: volOld, volNew CLASS(meshCell), POINTER:: cellOld, cellNew
!Assume that particle is outside the domain !Assume that particle is outside the domain
part%n_in = .FALSE. part%n_in = .FALSE.
volOld => mesh%vols(part%vol)%obj cellOld => mesh%cells(part%cell)%obj
CALL volOld%findCell(part) CALL cellOld%findCell(part)
CALL findCellColl(part) CALL findCellColl(part)
volNew => mesh%vols(part%vol)%obj
!Call the NA shcme !Call the NA shcme
IF (ASSOCIATED(self%weightingScheme)) THEN 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 END IF

26
src/modules/solver/pusher/modulePusher.f90
View file

@ -15,7 +15,7 @@ MODULE modulePusher
PURE SUBROUTINE pushCartElectrostatic(part, tauIn) PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
USE moduleSPecies USE moduleSPecies
USE moduleEM USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
@ -23,7 +23,8 @@ MODULE modulePusher
REAL(8):: qmEFt(1:3) REAL(8):: qmEFt(1:3)
!Get the electric field at particle position !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 !Update velocity
part%v = part%v + qmEFt part%v = part%v + qmEFt
@ -34,8 +35,8 @@ MODULE modulePusher
END SUBROUTINE pushCartElectrostatic END SUBROUTINE pushCartElectrostatic
PURE SUBROUTINE pushCartElectromagnetic(part, tauIn) PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
USE moduleSPecies USE moduleSpecies
USE moduleEM USE moduleMesh
USE moduleMath USE moduleMath
IMPLICIT NONE IMPLICIT NONE
@ -49,13 +50,14 @@ MODULE modulePusher
tauInHalf = tauIn *0.5D0 tauInHalf = tauIn *0.5D0
!Half of the force o f the electric field !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 !Half step for electrostatic
v_minus = part%v + qmEFt v_minus = part%v + qmEFt
!Full step rotation !Full step rotation
B = gatherMagnField(part) B = mesh%cells(part%cell)%obj%gatherMagneticField(part%Xi)
BNorm = NORM2(B) BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN IF (BNorm > 0.D0) THEN
fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
@ -112,7 +114,7 @@ MODULE modulePusher
!Push one particle. Boris pusher for 2D Cyl Electrostatic particle !Push one particle. Boris pusher for 2D Cyl Electrostatic particle
PURE SUBROUTINE push2DCylElectrostatic(part, tauIn) PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
USE moduleSpecies USE moduleSpecies
USE moduleEM USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
@ -124,7 +126,8 @@ MODULE modulePusher
part_temp = part part_temp = part
!Get electric field at particle position !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 !z
part_temp%v(1) = part%v(1) + qmEFt(1) part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn 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 !Push one particle. Boris pusher for 1D Radial Neutral particle
PURE SUBROUTINE push1DRadNeutral(part, tauIn) PURE SUBROUTINE push1DRadNeutral(part, tauIn)
USE moduleSpecies USE moduleSpecies
USE moduleEM
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
@ -188,7 +190,7 @@ MODULE modulePusher
!Push one particle. Boris pusher for 1D Radial Electrostatic particle !Push one particle. Boris pusher for 1D Radial Electrostatic particle
PURE SUBROUTINE push1DRadElectrostatic(part, tauIn) PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
USE moduleSpecies USE moduleSpecies
USE moduleEM USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
@ -200,7 +202,8 @@ MODULE modulePusher
part_temp = part part_temp = part
!Get electric field at particle position !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 !r,theta
v_p_oh_star(1) = part%v(1) + qmEFt(1) v_p_oh_star(1) = part%v(1) + qmEFt(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn x_new = part%r(1) + v_p_oh_star(1)*tauIn
@ -226,7 +229,6 @@ MODULE modulePusher
!Dummy pusher for 0D geometry !Dummy pusher for 0D geometry
PURE SUBROUTINE push0D(part, tauIn) PURE SUBROUTINE push0D(part, tauIn)
USE moduleSpecies USE moduleSpecies
USE moduleEM
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part