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

Files renamed and makefile make compatible with ifort.

Browse source
This commit is contained in:
Jorge Gonzalez 2020-12-10 19:25:17 +01:00
commit af74205932
25 changed files with 5439 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

489
src/modules/mesh/1DCart/moduleMesh1D.f90 Normal file
View file

@ -0,0 +1,489 @@
!moduleMesh1D: 1D cartesian module
! x == x
! y == unused
! z == unused
MODULE moduleMesh1D
USE moduleMesh
IMPLICIT NONE
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1D
!Element coordinates
REAL(8):: x = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initNode1D
PROCEDURE, PASS:: getCoordinates => getCoord1D
END TYPE meshNode1D
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1D
!Element coordinates
REAL(8):: x = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdge1D
PROCEDURE, PASS:: getNodes => getNodes1D
PROCEDURE, PASS:: randPos => randPos1D
END TYPE meshEdge1D
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1D
CONTAINS
PROCEDURE, PASS:: detJac => detJ1D
PROCEDURE, PASS:: invJac => invJ1D
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol1D
ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_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 meshVol1D
CLASS(meshVol1D), 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(meshVol1D):: meshVol1DSegm
!Element coordinates
REAL(8):: x(1:2)
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
!Connectivity to adjacent elements
CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS
PROCEDURE, PASS:: init => initVol1DSegm
PROCEDURE, PASS:: area => areaSegm
PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: partialDer => partialDerSegm
PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: weight => weightSegm
PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: scatter => scatterSegm
PROCEDURE, PASS:: gatherEF => gatherEFSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: nextElement => nextElementSegm
PROCEDURE, PASS:: resetOutput => resetOutputSegm
END TYPE meshVol1DSegm
CONTAINS
!NODE FUNCTIONS
!Init node element
SUBROUTINE initNode1D(self, n, r)
USE moduleSpecies
USE moduleRefParam
IMPLICIT NONE
CLASS(meshNode1D), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%x = r(1)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output
ALLOCATE(self%output(1:nSpecies))
END SUBROUTINE initNode1D
PURE FUNCTION getCoord1D(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNode1D), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%x, 0.D0, 0.D0 /)
END FUNCTION getCoord1D
!EDGE FUNCTIONS
!Inits edge element
SUBROUTINE initEdge1D(self, n, p, bt, physicalSurface)
IMPLICIT NONE
CLASS(meshEdge1D), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
REAL(8), DIMENSION(1:3):: r1
self%n = n
self%n1 => mesh%nodes(p(1))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
self%x = r1(1)
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
!Boundary index
self%bt = bt
!Physical Surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdge1D
!Get nodes from edge
PURE FUNCTION getNodes1D(self) RESULT(n)
IMPLICIT NONE
CLASS(meshEdge1D), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1))
n = (/ self%n1%n /)
END FUNCTION getNodes1D
!Calculates a 'random' position in edge
FUNCTION randPos1D(self) RESULT(r)
CLASS(meshEdge1D), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%x, 0.D0, 0.D0 /)
END FUNCTION randPos1D
!VOLUME FUNCTIONS
!SEGMENT FUNCTIONS
!Init segment element
SUBROUTINE initVol1DSegm(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8), DIMENSION(1:3):: r1, r2
self%n = n
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
self%x = (/ r1(1), r2(1) /)
!Assign node volume
CALL self%area()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initVol1DSegm
!Computes element area
PURE SUBROUTINE areaSegm(self)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xii(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!1 point Gauss integral
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
l = 2.D0*detJ
self%volume = l
self%arNodes = fPsi*l
END SUBROUTINE areaSegm
!Computes element functions at point xii
PURE FUNCTION fPsiSegm(xi) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:2))
fPsi(1) = 1.D0 - xi(1)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiSegm
!Computes element derivative shape function at Xii
PURE FUNCTION dPsiSegm(xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:1, 1:2))
dPsi(1, 1) = -5.D-1
dPsi(1, 2) = 5.D-1
END FUNCTION dPsiSegm
!Computes partial derivatives of coordinates
PURE SUBROUTINE partialDerSegm(self, dPsi, dx)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1):: dx
dx(1) = DOT_PRODUCT(dPsi(1,:), self%x)
END SUBROUTINE partialDerSegm
!Computes local stiffness matrix
PURE FUNCTION elemKSegm(self) RESULT(ke)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8):: ke(1:2,1:2)
REAL(8):: Xii(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ
ke = 0.D0
Xii = (/ 0.D0, 0.D0, 0.D0 /)
dPsi = self%dPsi(Xii)
invJ = self%invJac(Xii, dPsi)
ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /)
ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /)
ke = 2.D0*ke*invJ
END FUNCTION elemKSegm
PURE FUNCTION elemFSegm(self, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xii(1:3)
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
ALLOCATE(localF(1:2))
localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
END FUNCTION elemFSegm
PURE FUNCTION weightSegm(xi) RESULT(w)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: w(1:2)
w = fPsiSegm(xi)
END FUNCTION weightSegm
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
SUBROUTINE scatterSegm(self, part)
USE moduleOutput
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
CLASS(particle), INTENT(in):: part
TYPE(outputNode), POINTER:: vertex
REAL(8):: w_p(1:2)
REAL(8):: tensorS(1:3,1:3)
w_p = self%weight(part%xi)
tensorS = outerProduct(part%v, part%v)
vertex => self%n1%output(part%sp)
vertex%den = vertex%den + part%weight*w_p(1)
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
vertex => self%n2%output(part%sp)
vertex%den = vertex%den + part%weight*w_p(2)
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
END SUBROUTINE scatterSegm
!Gathers EF at position Xii
PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol1DSegm), 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
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVol1DSegm), 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(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
xN = 0.D0
xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
END FUNCTION phy2logSegm
!Get next element for a logical position xi
SUBROUTINE nextElementSegm(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
NULLIFY(nextElement)
IF (xi(1) < -1.D0) THEN
nextElement => self%e2
ELSEIF (xi(1) > 1.D0) THEN
nextElement => self%e1
END IF
END SUBROUTINE nextElementSegm
!Reset the output of nodes in element
PURE SUBROUTINE resetOutputSegm(self)
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout):: self
INTEGER:: k
DO k = 1, nSpecies
self%n1%output(k)%den = 0.D0
self%n1%output(k)%mom = 0.D0
self%n1%output(k)%tensorS = 0.D0
self%n2%output(k)%den = 0.D0
self%n2%output(k)%mom = 0.D0
self%n2%output(k)%tensorS = 0.D0
END DO
END SUBROUTINE resetOutputSegm
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes the element Jacobian determinant
PURE FUNCTION detJ1D(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol1D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
dJ = dx(1)
END FUNCTION detJ1D
!Computes the invers Jacobian
PURE FUNCTION invJ1D(self, xi, dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol1D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1)
REAL(8):: invJ
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
invJ = 1.D0/dx(1)
END FUNCTION invJ1D
END MODULE moduleMesh1D

40
src/modules/mesh/1DCart/moduleMesh1DBoundary.f90 Normal file
View file

@ -0,0 +1,40 @@
MODULE moduleMesh1DBoundary
USE moduleMesh1D
TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DRef
CONTAINS
PROCEDURE, PASS:: fBoundary => reflection
END TYPE meshEdge1DRef
TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DAbs
CONTAINS
PROCEDURE, PASS:: fBoundary => absorption
END TYPE meshEdge1DAbs
CONTAINS
SUBROUTINE reflection(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge1DRef), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
part%v(1) = -part%v(1)
part%r(1) = 2.D0*self%x - part%r(1)
END SUBROUTINE reflection
SUBROUTINE absorption(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge1DAbs), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
part%n_in = .FALSE.
END SUBROUTINE absorption
END MODULE moduleMesh1DBoundary

268
src/modules/mesh/1DCart/moduleMesh1DRead.f90 Normal file
View file

@ -0,0 +1,268 @@
MODULE moduleMesh1DRead
USE moduleMesh
USE moduleMesh1D
USE moduleMesh1DBoundary
!TODO: make this abstract to allow different mesh formats
TYPE, EXTENDS(meshGeneric):: mesh1DGeneric
CONTAINS
PROCEDURE, PASS:: init => init1DMesh
PROCEDURE, PASS:: readMesh => readMesh1D
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init 1D mesh
SUBROUTINE init1DMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(mesh1DGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1D")
END SELECT
END SUBROUTINE init1DMesh
!Reads 1D mesh
SUBROUTINE readMesh1D(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(mesh1DGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: x
INTEGER:: p(1:2)
INTEGER:: e, et, n, eTemp, elemType, bt
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Open file mesh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes coordinates. Only relevant for x
DO e = 1, self%numNodes
READ(10, *) n, x
ALLOCATE(meshNode1D:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads the total number of elements (edges+vol)
READ(10, *) totalNumElem
self%numEdges = 0
DO e = 1, totalNumElem
READ(10, *) eTemp, elemType
IF (elemType == 15) THEN !15 is physical node in GMSH
self%numEdges = e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = totalNumelem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beginning of reading elements
DO e = 1, totalNumelem
BACKSPACE(10)
END DO
!Reads edges
DO e = 1, self%numEdges
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
SELECT CASE(boundary(bt)%obj%boundaryType)
CASE ('reflection')
ALLOCATE(meshEdge1DRef:: self%edges(e)%obj)
CASE ('absorption')
ALLOCATE(meshEdge1DAbs:: self%edges(e)%obj)
END SELECT
CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType)
END DO
!Read and initialize volumes
DO e = 1, self%numVols
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DSegm:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity betwen vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMesh1D
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVol1DSegm)
SELECT TYPE(elemB)
TYPE IS(meshVol1DSegm)
CALL connectedSegmSegm(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedSegmSegm(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectedSegmSegm
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS (meshVol1DSegm)
SELECT TYPE(elemB)
CLASS IS(meshEdge1D)
CALL connectedSegmEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedSegmEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
CLASS(meshEdge1D), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectedSegmEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8):: localK(1:2,1:2)
INTEGER:: i, j
INTEGER:: n(1:2)
SELECT TYPE(elem)
TYPE IS(meshVol1DSegm)
localK = elem%elemK()
n = (/ elem%n1%n, elem%n2%n /)
CLASS DEFAULT
n = 0
localK = 0.D0
END SELECT
DO i = 1, 2
DO j = 1, 2
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMesh1DRead

11
src/modules/mesh/2DCyl/makefile Normal file
View file

@ -0,0 +1,11 @@
all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o
moduleMeshCyl.o: moduleMeshCyl.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

1022
src/modules/mesh/2DCyl/moduleMeshCyl.f90 Normal file
View file

File diff suppressed because it is too large Load diff

80
src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 Normal file
View file

@ -0,0 +1,80 @@
!moduleMeshCylBoundary: Edge elements for Cylindrical mesh.
MODULE moduleMeshCylBoundary
USE moduleMeshCyl
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylRef
CONTAINS
PROCEDURE, PASS:: fBoundary => reflection
END TYPE meshEdgeCylRef
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAbs
CONTAINS
PROCEDURE, PASS:: fBoundary => absorption
END TYPE meshEdgeCylAbs
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAxis
CONTAINS
PROCEDURE, PASS:: fBoundary => symmetryAxis
END TYPE meshEdgeCylAxis
CONTAINS
SUBROUTINE reflection(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylRef), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
edgeNorm = DSQRT((self%r(2)-self%r(1))**2 + (self%z(2)-self%z(1))**2)
cosT = (self%z(2)-self%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - self%z(1);
rp(2) = part%r(2) - self%r(1);
rpp(1) = cosT*rp(1) - sinT*rp(2)
rpp(2) = sinT*rp(1) + cosT*rp(2)
rpp(2) = -rpp(2)
vpp(1) = cosT*part%v(1) - sinT*part%v(2)
vpp(2) = sinT*part%v(1) + cosT*part%v(2)
vpp(2) = -vpp(2)
part%r(1) = cosT*rpp(1) + sinT*rpp(2) + self%z(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + self%r(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAbs), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
!TODO: Add scatter to mesh nodes
part%n_in = .FALSE.
END SUBROUTINE absorption
SUBROUTINE symmetryAxis(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAxis), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
END SUBROUTINE symmetryAxis
END MODULE moduleMeshCylBoundary

610
src/modules/mesh/2DCyl/moduleMeshCylRead.f90 Normal file
View file

@ -0,0 +1,610 @@
MODULE moduleMeshCylRead
USE moduleMesh
USE moduleMeshCyl
USE moduleMeshCylBoundary
!TODO: make this abstract to allow different mesh formats
TYPE, EXTENDS(meshGeneric):: meshCylGeneric
CONTAINS
PROCEDURE, PASS:: init => initCylMesh
PROCEDURE, PASS:: readMesh => readMeshCylGmsh
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init mesh
SUBROUTINE initCylMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "initCylMesh")
END SELECT
END SUBROUTINE initCylMesh
!Read mesh from gmsh file
SUBROUTINE readMeshCylGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r, z
INTEGER:: p(1:4)
INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read msh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes cartesian coordinates (x=z, y=r, z=null)
DO e=1, self%numNodes
READ(10, *) n, z, r
ALLOCATE(meshNodeCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads Totalnumber of elements
READ(10, *) TotalnumElem
!counts edges and volume elements
self%numEdges = 0
DO e=1, TotalnumElem
READ(10,*) eTemp, elemType
IF (elemType==1) THEN
self%numEdges=e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read elements
DO e=1, totalNumElem
BACKSPACE(10)
END DO
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
SELECT CASE(boundary(bt)%obj%boundaryType)
CASE ('reflection')
ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj)
CASE ('absorption')
ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj)
CASE ('axis')
ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj)
END SELECT
CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
END DO
!Read and initialize volumes
DO e=1, self%numVols
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVolCylTria:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
CASE (3)
!Quadrilateral element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity between vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMeshCylGmsh
!Selects type of elements to build connection
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
!Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
!Element B is a quadrilateral
CALL connectedQuadQuad(elemA, elemB)
TYPE IS(meshVolCylTria)
!Element B is a triangle
CALL connectedQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshVolCylTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
!Element B is a quadrilateral
CALL connectedQuadTria(elemB, elemA)
TYPE IS(meshVolCylTria)
!Element B is a triangle
CALL connectedTriaTria(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemB)
CLASS IS(meshEdgeCyl)
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
!Element A is a quadrilateral
CALL connectedQuadEdge(elemA, elemB)
TYPE IS(meshVolCylTria)
!Element A is a triangle
CALL connectedTriaEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n4%n) THEN
elemA%e2 => elemB
elemB%e4 => elemA
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
END IF
END SUBROUTINE connectedQuadQuad
SUBROUTINE connectedQuadTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e4 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadTria
SUBROUTINE connectedTriaTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
END SUBROUTINE connectedTriaTria
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadEdge
SUBROUTINE connectedTriaEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedTriaEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVolCylQuad)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
TYPE IS(meshVolCylTria)
nNodes = 3
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0
ALLOCATE(n(1:1))
n = 0
END SELECT
DO i = 1, nNodes
DO j = 1, nNodes
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMeshCylRead

609
src/modules/mesh/moduleMesh.f90 Normal file
View file

@ -0,0 +1,609 @@
!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
USE moduleList
USE moduleOutput
IMPLICIT NONE
!Parent of Node element
TYPE, PUBLIC, ABSTRACT:: meshNode
!Node index
INTEGER:: n = 0
!Node volume
REAL(8):: v = 0.D0
!Output values
TYPE(outputNode), ALLOCATABLE:: output(:)
TYPE(emNode):: emData
CONTAINS
PROCEDURE(initNode_interface), DEFERRED, PASS:: init
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
END TYPE meshNode
ABSTRACT INTERFACE
!Interface of init a node (3D generic coordinates)
SUBROUTINE initNode_interface(self, n, r)
IMPORT:: meshNode
CLASS(meshNode), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
END SUBROUTINE initNode_interface
!Interface to get coordinates from node
PURE FUNCTION getCoord_interface(self) RESULT(r)
IMPORT:: meshNode
CLASS(meshNode), INTENT(in):: self
REAL(8):: r(1:3)
END FUNCTION getCoord_interface
END INTERFACE
!Containers for nodes in the mesh
TYPE:: meshNodeCont
CLASS(meshNode), ALLOCATABLE:: obj
END TYPE meshNodeCont
!Parent of Edge element
TYPE, PUBLIC, ABSTRACT:: meshEdge
!Element index
INTEGER:: n = 0
!Connectivity to vols
CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL()
!Normal vector
REAL(8):: normal(1:3)
!Physical surface in mesh
INTEGER:: physicalSurface
!id for boundary condition
INTEGER:: bt = 0
CONTAINS
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary
PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos
END TYPE meshEdge
ABSTRACT INTERFACE
SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
END SUBROUTINE initEdge_interface
SUBROUTINE boundary_interface(self, part)
USE moduleSpecies
IMPORT:: meshEdge
CLASS (meshEdge), INTENT(inout):: self
CLASS (particle), INTENT(inout):: part
END SUBROUTINE
PURE FUNCTION getNodesEdge_interface(self) RESULT(n)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION
FUNCTION randPos_interface(self) RESULT(r)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
REAL(8):: r(1:3)
END FUNCTION randPos_interface
END INTERFACE
!Containers for edges in the mesh
TYPE:: meshEdgeCont
CLASS(meshEdge), ALLOCATABLE:: obj
END TYPE meshEdgeCont
!Parent of Volume element
TYPE, PUBLIC, ABSTRACT:: meshVol
!Volume index
INTEGER:: n = 0
!Maximum collision rate
REAL(8):: sigmaVrelMax = 1.D-15
!Volume
REAL(8):: volume = 0.D0
!List of particles inside the volume
TYPE(listNode):: listPart_in
!Lock indicator for listPart_in
INTEGER(KIND=OMP_LOCK_KIND):: lock
!Number of collisions per volume
INTEGER:: nColl = 0
!Total weight of particles inside cell
REAL(8):: totalWeight = 0.D0
CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
PROCEDURE, PASS:: findCell
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside
PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement
PROCEDURE, PASS:: collision
PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput
END TYPE meshVol
ABSTRACT INTERFACE
SUBROUTINE initVol_interface(self, n, p)
IMPORT:: meshVol
CLASS(meshVol), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
END SUBROUTINE initVol_interface
SUBROUTINE scatter_interface(self, part)
USE moduleSpecies
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
CLASS(particle), INTENT(in):: part
END SUBROUTINE scatter_interface
PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: EF(1:3)
END FUNCTION gatherEF_interface
PURE FUNCTION getNodesVol_interface(self) RESULT(n)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION getNodesVol_interface
PURE FUNCTION elemF_interface(self, source) RESULT(localF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
END FUNCTION elemF_interface
SUBROUTINE nextElement_interface(self, xi, nextElement)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), 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
END FUNCTION inside_interface
SUBROUTINE collision_interface(self)
IMPORT:: meshVol
CLASS(meshVol), INTENT(inout):: self
END SUBROUTINE collision_interface
PURE SUBROUTINE resetOutput_interface(self)
IMPORT:: meshVol
CLASS(meshVol), INTENT(inout):: self
END SUBROUTINE resetOutput_interface
END INTERFACE
!Containers for volumes in the mesh
TYPE:: meshVolCont
CLASS(meshVol), ALLOCATABLE:: obj
END TYPE meshVolCont
!Abstract type of mesh
TYPE, PUBLIC, ABSTRACT:: meshGeneric
INTEGER:: numEdges, numNodes, numVols
!Array of nodes
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of boundary elements
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
!Array of volume elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:)
!Global stiffness matrix
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Permutation matrix for P L U factorization
INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
CONTAINS
PROCEDURE(initMesh_interface), DEFERRED, PASS:: init
PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh
END TYPE meshGeneric
ABSTRACT INTERFACE
!Inits the mesh
SUBROUTINE initMesh_interface(self, meshFormat)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
END SUBROUTINE initMesh_interface
!Reads the mesh from a file
SUBROUTINE readMesh_interface(self, filename)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
END SUBROUTINE readMesh_interface
!Prints Species data
SUBROUTINE printOutput_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints number of collisions
SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
!Prints EM info
SUBROUTINE printEM_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface
END INTERFACE
!Generic mesh
CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh
CONTAINS
!Find next cell for particle
RECURSIVE SUBROUTINE findCell(self, part, oldCell)
USE moduleSpecies
USE OMP_LIB
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
CLASS(particle), INTENT(inout), TARGET:: part
REAL(8):: xi(1:3)
CLASS(*), POINTER:: nextElement
xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell
IF (self%inside(xi)) THEN
part%vol = self%n
part%xi = xi
part%n_in = .TRUE.
!Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock)
CALL self%listPart_in%add(part)
self%totalWeight = self%totalWeight + part%weight
CALL OMP_UNSET_LOCK(self%lock)
ELSE
!If not, searches for a neighbour and repeats the process.
CALL self%nextElement(xi, nextElement)
!Defines the next step
SELECT TYPE(nextElement)
CLASS IS(meshVol)
!Particle moved to new cell, repeat find procedure
CALL nextElement%findCell(part, self)
CLASS IS (meshEdge)
!Particle encountered an edge, execute boundary
CALL nextElement%fBoundary(part)
!If particle is still inside the domain, call findCell
IF (part%n_in) THEN
IF(PRESENT(oldCell)) THEN
CALL self%findCell(part, oldCell)
ELSE
CALL self%findCell(part)
END IF
END IF
CLASS DEFAULT
WRITE(*,*) "ERROR, CHECK findCell"
END SELECT
END IF
END SUBROUTINE findCell
!Computes collisions in element
SUBROUTINE collision(self)
USE moduleCollisions
USE moduleSpecies
USE moduleList
use moduleRefParam
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
INTEGER:: nPart !Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision
INTEGER:: rnd !random index
REAL(8):: rndReal
TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision
INTEGER:: ij, k
REAL(8):: sigmaVrelMaxNew
TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
self%nColl = 0
nPart = self%listPart_in%amount
IF (nPart > 1) THEN
pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
self%nColl = INT(REAL(nPart)*pMax*0.5D0)
!Converts the list of particles to an array for easy access
IF (self%nColl > 0) THEN
partTemp = self%listPart_in%convert2Array()
END IF
DO n = 1, self%nColl
!Select random numbers
CALL RANDOM_NUMBER(rndReal)
rnd = 1 + FLOOR(nPart*rndReal)
part_i => partTemp(rnd)%part
CALL RANDOM_NUMBER(rndReal)
rnd = 1 + FLOOR(nPart*rndReal)
part_j => partTemp(rnd)%part
ij = interactionIndex(part_i%sp, part_j%sp)
sigmaVrelMaxNew = 0.D0
DO k = 1, interactionMatrix(ij)%amount
CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
END DO
!Update maximum cross section*v_rel per each collision
IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
self%sigmaVrelMax = sigmaVrelMaxNew
END IF
END DO
END IF
self%totalWeight = 0.D0
!Reset output in nodes
CALL self%resetOutput()
!Erase the list of particles inside the cell
CALL self%listPart_in%erase()
END SUBROUTINE collision
SUBROUTINE printOutputGmsh(self, t)
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
time = DBLE(t)*tauMin*ti_ref
DO i = 1, nSpecies
WRITE(tstring, '(I6.6)') t
fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Density (m^-3)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Velocity (m/s)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 3
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Pressure (Pa)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Temperature (K)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
WRITE(60, "(A)") '$EndNodeData'
CLOSE (60)
END DO
END SUBROUTINE printOutputGmsh
SUBROUTINE printCollGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, '(I6.6)') t
fileName='OUTPUT_' // tstring// '_Collisions.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
END DO
WRITE(60, "(A)") '$EndElementData'
CLOSE(60)
END IF
END SUBROUTINE printCollGmsh
SUBROUTINE printEMGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
REAL(8):: xi(1:3)
xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, '(I6.6)') t
fileName='OUTPUT_' // tstring// '_EMField.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (20, file = path // folder // '/' // fileName)
WRITE(20, "(A)") '$MeshFormat'
WRITE(20, "(A)") '2.2 0 8'
WRITE(20, "(A)") '$EndMeshFormat'
WRITE(20, "(A)") '$NodeData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Potential (V)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 1
WRITE(20, *) self%numNodes
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO
WRITE(20, "(A)") '$EndNodeData'
WRITE(20, "(A)") '$ElementData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Electric Field (V/m)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 3
WRITE(20, *) self%numVols
DO e=1, self%numVols
WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
END DO
WRITE(20, "(A)") '$EndElementData'
CLOSE(20)
END IF
END SUBROUTINE printEMGmsh
END MODULE moduleMesh