Jorge Gonzalez
db6b0a2c03
assigning the normal vector in 2D and 3D. 3D Cartesian geometry is working properly, although it needs testing. Still issue with ionization boundary.
527 lines
14 KiB
Fortran
527 lines
14 KiB
Fortran
!moduleMesh1D: 1D cartesian module
|
|
! x == x
|
|
! y == unused
|
|
! z == unused
|
|
MODULE moduleMesh1DCart
|
|
USE moduleMesh
|
|
USE moduleMeshBoundary
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
|
|
REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /)
|
|
|
|
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
|
|
!Element coordinates
|
|
REAL(8):: x = 0.D0
|
|
CONTAINS
|
|
PROCEDURE, PASS:: init => initNode1DCart
|
|
PROCEDURE, PASS:: getCoordinates => getCoord1DCart
|
|
|
|
END TYPE meshNode1DCart
|
|
|
|
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge1DCart
|
|
!Element coordinates
|
|
REAL(8):: x = 0.D0
|
|
!Connectivity to nodes
|
|
CLASS(meshNode), POINTER:: n1 => NULL()
|
|
CONTAINS
|
|
PROCEDURE, PASS:: init => initEdge1DCart
|
|
PROCEDURE, PASS:: getNodes => getNodes1DCart
|
|
PROCEDURE, PASS:: intersection => intersection1DCart
|
|
PROCEDURE, PASS:: randPos => randPosEdge
|
|
|
|
END TYPE meshEdge1DCart
|
|
|
|
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart
|
|
CONTAINS
|
|
PROCEDURE, PASS:: detJac => detJ1DCart
|
|
PROCEDURE, PASS:: invJac => invJ1DCart
|
|
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
|
|
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
|
|
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
|
|
|
|
END TYPE meshVol1DCart
|
|
|
|
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 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
|
|
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 => initVol1DCartSegm
|
|
PROCEDURE, PASS:: randPos => randPos1DCartSeg
|
|
PROCEDURE, PASS:: area => areaSegm
|
|
PROCEDURE, NOPASS:: fPsi => fPsiSegm
|
|
PROCEDURE, NOPASS:: dPsi => dPsiSegm
|
|
PROCEDURE, PASS:: partialDer => partialDerSegm
|
|
PROCEDURE, PASS:: elemK => elemKSegm
|
|
PROCEDURE, PASS:: elemF => elemFSegm
|
|
PROCEDURE, NOPASS:: 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
|
|
|
|
END TYPE meshVol1DCartSegm
|
|
|
|
CONTAINS
|
|
!NODE FUNCTIONS
|
|
!Init node element
|
|
SUBROUTINE initNode1DCart(self, n, r)
|
|
USE moduleSpecies
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshNode1DCart), 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 initNode1DCart
|
|
|
|
PURE FUNCTION getCoord1DCart(self) RESULT(r)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshNode1DCart), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
r = (/ self%x, 0.D0, 0.D0 /)
|
|
|
|
END FUNCTION getCoord1DCart
|
|
|
|
!EDGE FUNCTIONS
|
|
!Inits edge element
|
|
SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface)
|
|
USE moduleSpecies
|
|
USE moduleBoundary
|
|
USE moduleErrors
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge1DCart), 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
|
|
INTEGER:: s
|
|
|
|
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%boundary => boundary(bt)
|
|
ALLOCATE(self%fboundary(1:nSpecies))
|
|
!Assign functions to boundary
|
|
DO s = 1, nSpecies
|
|
CALL pointBoundaryFunction(self, s)
|
|
|
|
END DO
|
|
|
|
!Physical Surface
|
|
self%physicalSurface = physicalSurface
|
|
|
|
END SUBROUTINE initEdge1DCart
|
|
|
|
!Get nodes from edge
|
|
PURE FUNCTION getNodes1DCart(self) RESULT(n)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge1DCart), INTENT(in):: self
|
|
INTEGER, ALLOCATABLE:: n(:)
|
|
|
|
ALLOCATE(n(1))
|
|
n = (/ self%n1%n /)
|
|
|
|
END FUNCTION getNodes1DCart
|
|
|
|
PURE FUNCTION intersection1DCart(self, r0) RESULT(r)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge1DCart), INTENT(in):: self
|
|
REAL(8), DIMENSION(1:3), INTENT(in):: r0
|
|
REAL(8), DIMENSION(1:3):: r
|
|
|
|
r = (/ self%x, 0.D0, 0.D0 /)
|
|
|
|
END FUNCTION intersection1DCart
|
|
|
|
!Calculates a 'random' position in edge
|
|
FUNCTION randPosEdge(self) RESULT(r)
|
|
CLASS(meshEdge1DCart), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
r = (/ self%x, 0.D0, 0.D0 /)
|
|
|
|
END FUNCTION randPosEdge
|
|
|
|
!VOLUME FUNCTIONS
|
|
!SEGMENT FUNCTIONS
|
|
!Init segment element
|
|
SUBROUTINE initVol1DCartSegm(self, n, p)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), 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 initVol1DCartSegm
|
|
|
|
!Calculates a random position in 1D volume
|
|
FUNCTION randPos1DCartSeg(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: xii(1:3)
|
|
REAL(8), ALLOCATABLE:: fPsi(:)
|
|
|
|
xii(1) = random(-1.D0, 1.D0)
|
|
xii(2:3) = 0.D0
|
|
|
|
fPsi = self%fPsi(xii)
|
|
r(1) = DOT_PRODUCT(fPsi, self%x)
|
|
|
|
END FUNCTION randPos1DCartSeg
|
|
|
|
!Computes element area
|
|
PURE SUBROUTINE areaSegm(self)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), INTENT(inout):: self
|
|
REAL(8):: l !element length
|
|
REAL(8):: fPsi(1:2)
|
|
REAL(8):: detJ
|
|
REAL(8):: Xii(1:3)
|
|
|
|
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(meshVol1DCartSegm), 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
|
|
FUNCTION elemKSegm(self) RESULT(ke)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), 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(1), detJ
|
|
INTEGER:: l
|
|
|
|
ke = 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)
|
|
ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
|
|
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
|
|
wSeg(l)/detJ
|
|
|
|
END DO
|
|
|
|
END FUNCTION elemKSegm
|
|
|
|
PURE FUNCTION elemFSegm(self, source) RESULT(localF)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), INTENT(in):: self
|
|
REAL(8), INTENT(in):: source(1:)
|
|
REAL(8), ALLOCATABLE:: localF(:)
|
|
REAL(8):: fPsi(1:2)
|
|
REAL(8):: detJ, f
|
|
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)
|
|
f = DOT_PRODUCT(fPsi, source)
|
|
localF = localF + f*fPsi*wSeg(l)*detJ
|
|
|
|
END DO
|
|
|
|
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 moduleMath
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCartSegm), 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(meshVol1DCartSegm), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8):: dPsi(1, 1:2)
|
|
REAL(8):: phi(1:2)
|
|
REAL(8):: EF(1:3)
|
|
REAL(8):: invJ
|
|
|
|
phi = (/ self%n1%emData%phi, &
|
|
self%n2%emData%phi /)
|
|
|
|
dPsi = self%dPsi(xi)
|
|
invJ = self%invJac(xi, dPsi)
|
|
EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
|
|
EF(2) = 0.D0
|
|
EF(3) = 0.D0
|
|
|
|
END FUNCTION gatherEFSegm
|
|
|
|
!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):: 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(meshVol1DCartSegm), 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
|
|
|
|
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
|
|
!Calculates a random position in 1D volume
|
|
!Computes the element Jacobian determinant
|
|
PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCart), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
REAL(8):: 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 detJ1DCart
|
|
|
|
!Computes the invers Jacobian
|
|
PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol1DCart), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
REAL(8):: dx(1)
|
|
REAL(8):: invJ
|
|
|
|
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 invJ1DCart
|
|
|
|
|
|
END MODULE moduleMesh1DCart
|
|
|