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

fPsi no longer allocates memory

Browse source 
I noticed that phy2logquad had a lot of overhead. Trying to reducing it
by simplifying calls to fPsi, dPsi and such.

The function for fPsi has been made so no memory is allocated and works
under the assumption that the input array has the right size (1:numNodes)
This commit is contained in:
Jorge Gonzalez 2023-01-01 12:12:06 +01:00
commit 0db76083ec
8 changed files with 408 additions and 409 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -361,14 +361,14 @@ MODULE moduleInput
!Density at centroid of cell
nodes = mesh%vols(e)%obj%getNodes()
nNodes = SIZE(nodes)
fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
ALLOCATE(fPsi(1:nNodes))
CALL mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), fPsi)
ALLOCATE(source(1:nNodes))
DO j = 1, nNodes
source(j) = density(nodes(j))
END DO
densityCen = DOT_PRODUCT(fPsi, source)
DEALLOCATE(fPsi)
!Calculate number of particles
nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
@ -380,7 +380,7 @@ MODULE moduleInput
partNew%r = mesh%vols(e)%obj%randPos()
partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r)
!Get mean velocity at particle position
fPsi = mesh%vols(e)%obj%fPsi(partNew%xi)
CALL mesh%vols(e)%obj%fPsi(partNew%xi, fPsi)
DO j = 1, nNodes
source(j) = velocity(nodes(j), 1)

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

@ -106,14 +106,13 @@ MODULE moduleMesh0D
END FUNCTION randPos0D
PURE FUNCTION fPsi0D(xi) RESULT(fPsi)
PURE SUBROUTINE fPsi0D(xi, fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(out):: fPsi(:)
ALLOCATE(fPsi(1:1))
fPsi = 1.D0
END FUNCTION fPsi0D
END SUBROUTINE fPsi0D
PURE FUNCTION gatherEF0D(self, xi) RESULT(EF)
IMPLICIT NONE

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

@ -230,13 +230,13 @@ MODULE moduleMesh1DCart
CLASS(meshVol1DCartSegm), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
xii(1) = random(-1.D0, 1.D0)
xii(2:3) = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2:3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(Xi, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%x)
END FUNCTION randPos1DCartSeg
@ -249,36 +249,34 @@ MODULE moduleMesh1DCart
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xii(1:3)
REAL(8):: Xi(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)
Xi = 0.D0
CALL self%fPsi(Xi, fPsi)
detJ = self%detJac(Xi)
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)
!Computes element functions at point Xi
PURE SUBROUTINE fPsiSegm(xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:2))
REAL(8), INTENT(out):: fPsi(:)
fPsi(1) = 1.D0 - xi(1)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiSegm
END SUBROUTINE fPsiSegm
!Computes element derivative shape function at Xii
!Computes element derivative shape function at Xi
PURE FUNCTION dPsiSegm(xi) RESULT(dPsi)
IMPLICIT NONE
@ -310,19 +308,19 @@ MODULE moduleMesh1DCart
CLASS(meshVol1DCartSegm), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: Xii(1:3)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ(1), detJ
INTEGER:: l
ALLOCATE(localK(1:2,1:2))
localK = 0.D0
Xii = 0.D0
Xi = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
dPsi = self%dPsi(Xii)
detJ = self%detJac(Xii, dPsi)
invJ = self%invJac(Xii, dPsi)
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi, dPsi)
invJ = self%invJac(Xi, dPsi)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
wSeg(l)/detJ
@ -339,17 +337,17 @@ MODULE moduleMesh1DCart
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ, f
REAL(8):: Xii(1:3)
REAL(8):: Xi(1:3)
INTEGER:: l
ALLOCATE(localF(1:2))
localF = 0.D0
Xii = 0.D0
Xi = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
detJ = self%detJac(Xii)
fPsi = self%fPsi(Xii)
Xi(1) = corSeg(l)
detJ = self%detJac(Xi)
CALL self%fPsi(Xi, fPsi)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*wSeg(l)*detJ
@ -368,7 +366,7 @@ MODULE moduleMesh1DCart
END FUNCTION insideSegm
!Gathers EF at position Xii
!Gathers EF at position Xi
PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF)
IMPLICIT NONE
@ -407,7 +405,7 @@ MODULE moduleMesh1DCart
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFSegm

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

@ -232,13 +232,13 @@ MODULE moduleMesh1DRad
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
xii(1) = random(-1.D0, 1.D0)
xii(2:3) = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2:3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(Xi, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%r)
END FUNCTION randPos1DRadSeg
@ -249,47 +249,47 @@ MODULE moduleMesh1DRad
CLASS(meshVol1DRadSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: fPsi(1:2), fPsi_node(1:2)
REAL(8):: r
REAL(8):: detJ
REAL(8):: Xii(1:3)
REAL(8):: Xi(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)
Xi = 0.D0
CALL self%fPsi(Xi, fPsi)
detJ = self%detJac(Xi)
!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)
Xi = (/-5.D-1, 0.D0, 0.D0/)
CALL self%fPsi(Xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,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)
Xi = (/ 5.D-1, 0.D0, 0.D0/)
CALL self%fPsi(Xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,self%r)
self%arNodes(2) = fPsi(2)*r*l
END SUBROUTINE areaRad
!Computes element functions at point xii
PURE FUNCTION fPsiRad(xi) RESULT(fPsi)
!Computes element functions at point Xi
PURE SUBROUTINE fPsiRad(xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:2))
REAL(8), INTENT(out):: fPsi(:)
fPsi(1) = 1.D0 - xi(1)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiRad
END SUBROUTINE fPsiRad
!Computes element derivative shape function at Xii
!Computes element derivative shape function at Xi
PURE FUNCTION dPsiRad(xi) RESULT(dPsi)
IMPLICIT NONE
@ -322,7 +322,7 @@ MODULE moduleMesh1DRad
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: Xii(1:3)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ(1), detJ
REAL(8):: r, fPsi(1:2)
@ -330,13 +330,13 @@ MODULE moduleMesh1DRad
ALLOCATE(localK(1:2, 1:2))
localK = 0.D0
Xii = 0.D0
Xi = 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)
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi, dPsi)
invJ = self%invJac(Xi, dPsi)
CALL self%fPsi(Xi, fPsi)
r = DOT_PRODUCT(fPsi, self%r)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
@ -357,17 +357,17 @@ MODULE moduleMesh1DRad
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ, f, r
REAL(8):: Xii(1:3)
REAL(8):: Xi(1:3)
INTEGER:: l
ALLOCATE(localF(1:2))
localF = 0.D0
Xii = 0.D0
Xi = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
detJ = self%detJac(Xii)
fPsi = self%fPsi(Xii)
Xi(1) = corSeg(l)
detJ = self%detJac(Xi)
CALL self%fPsi(Xi, fPsi)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*r*wSeg(l)*detJ
@ -387,7 +387,7 @@ MODULE moduleMesh1DRad
END FUNCTION insideRad
!Gathers EF at position Xii
!Gathers EF at position Xi
PURE FUNCTION gatherEFRad(self, xi) RESULT(EF)
IMPLICIT NONE
@ -426,7 +426,7 @@ MODULE moduleMesh1DRad
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFRad

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

@ -11,8 +11,8 @@ MODULE moduleMesh2DCart
REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /)
REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /)
REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /)
REAL(8), PARAMETER:: Xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /)
REAL(8), PARAMETER:: Xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /)
REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCart
@ -47,8 +47,8 @@ MODULE moduleMesh2DCart
END TYPE meshVol2DCart
ABSTRACT INTERFACE
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3)
PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
END FUNCTION dPsi_interface
@ -210,14 +210,14 @@ MODULE moduleMesh2DCart
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
xii(1) = random(-1.D0, 1.D0)
xii(2) = random(-1.D0, 1.D0)
xii(3) = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(Xi, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
@ -319,78 +319,77 @@ MODULE moduleMesh2DCart
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(inout):: self
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:4)
self%volume = 0.D0
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
xi = 0.D0
detJ = self%detJac(xi)*4.D0 !4
fPsi = self%fPsi(xi)
Xi = 0.D0
detJ = self%detJac(Xi)*4.D0 !4
CALL self%fPsi(Xi, fPsi)
self%volume = detJ
self%arNodes = fPsi*detJ
END SUBROUTINE areaQuad
!Computes element functions in point xi
PURE FUNCTION fPsiQuad(xi) RESULT(fPsi)
!Computes element functions in point Xi
PURE SUBROUTINE fPsiQuad(Xi, fPsi)
IMPLICIT NONE
REAL(8),INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(out):: fPsi(:)
ALLOCATE(fPsi(1:4))
fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2))
fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2))
fPsi(3) = (1.D0+Xi(1)) * (1.D0+Xi(2))
fPsi(4) = (1.D0-Xi(1)) * (1.D0+Xi(2))
fPsi(1) = (1.D0-xi(1))*(1.D0-xi(2))
fPsi(2) = (1.D0+xi(1))*(1.D0-xi(2))
fPsi(3) = (1.D0+xi(1))*(1.D0+xi(2))
fPsi(4) = (1.D0-xi(1))*(1.D0+xi(2))
fPsi = fPsi*0.25D0
END FUNCTION fPsiQuad
END SUBROUTINE fPsiQuad
!Derivative element function at coordinates xi
PURE FUNCTION dPsiQuad(xi) RESULT(dPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiQuad(Xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:2,1:4))
dPsi(1,:) = dPsiQuadXi1(xi(2))
dPsi(2,:) = dPsiQuadXi2(xi(1))
dPsi(1,:) = dPsiQuadXi1(Xi(2))
dPsi(2,:) = dPsiQuadXi2(Xi(1))
END FUNCTION dPsiQuad
!Derivative element function (xi1)
PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1)
!Derivative element function (Xi1)
PURE FUNCTION dPsiQuadXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8),INTENT(in):: xi2
REAL(8),INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:4)
dPsiXi1(1) = -(1.D0-xi2)
dPsiXi1(2) = (1.D0-xi2)
dPsiXi1(3) = (1.D0+xi2)
dPsiXi1(4) = -(1.D0+xi2)
dPsiXi1(1) = -(1.D0-Xi2)
dPsiXi1(2) = (1.D0-Xi2)
dPsiXi1(3) = (1.D0+Xi2)
dPsiXi1(4) = -(1.D0+Xi2)
dPsiXi1 = dPsiXi1*0.25D0
END FUNCTION dPsiQuadXi1
!Derivative element function (xi2)
PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2)
!Derivative element function (Xi2)
PURE FUNCTION dPsiQuadXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8),INTENT(in):: xi1
REAL(8),INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:4)
dPsiXi2(1) = -(1.D0-xi1)
dPsiXi2(2) = -(1.D0+xi1)
dPsiXi2(3) = (1.D0+xi1)
dPsiXi2(4) = (1.D0-xi1)
dPsiXi2(1) = -(1.D0-Xi1)
dPsiXi2(2) = -(1.D0+Xi1)
dPsiXi2(3) = (1.D0+Xi1)
dPsiXi2(4) = (1.D0-Xi1)
dPsiXi2 = dPsiXi2*0.25D0
END FUNCTION dPsiQuadXi2
@ -415,24 +414,24 @@ MODULE moduleMesh2DCart
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:2,1:4)
REAL(8):: invJ(1:2,1:2), detJ
INTEGER:: l, m
ALLOCATE(localK(1:4, 1:4))
localK=0.D0
xi=0.D0
Xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 3
xi(2) = corQuad(l)
dPsi(1,:) = self%dPsiXi1(xi(2))
Xi(2) = corQuad(l)
dPsi(1,:) = self%dPsiXi1(Xi(2))
DO m = 1, 3
xi(1) = corQuad(m)
dPsi(2,:) = self%dPsiXi2(xi(1))
fPsi = self%fPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
Xi(1) = corQuad(m)
dPsi(2,:) = self%dPsiXi2(Xi(1))
CALL self%fPsi(Xi, fPsi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ
END DO
@ -447,20 +446,20 @@ MODULE moduleMesh2DCart
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: detJ, f
INTEGER:: l, m
ALLOCATE(localF(1:4))
localF = 0.D0
xi = 0.D0
Xi = 0.D0
DO l=1, 3
xi(1) = corQuad(l)
Xi(1) = corQuad(l)
DO m = 1, 3
xi(2) = corQuad(m)
detJ = self%detJac(xi)
fPsi = self%fPsi(xi)
Xi(2) = corQuad(m)
detJ = self%detJac(Xi)
CALL self%fPsi(Xi, fPsi)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ
@ -470,23 +469,23 @@ MODULE moduleMesh2DCart
END FUNCTION elemFQuad
!Checks if a particle is inside a quad element
PURE FUNCTION insideQuad(xi) RESULT(ins)
PURE FUNCTION insideQuad(Xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins
ins = (xi(1) >= -1.D0 .AND. xi(1) <= 1.D0) .AND. &
(xi(2) >= -1.D0 .AND. xi(2) <= 1.D0)
ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. &
(Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0)
END FUNCTION insideQuad
!Gathers the electric field at position xi
PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)
!Gathers the electric field at position Xi
PURE FUNCTION gatherEFQuad(self,Xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: dPsi(1:2,1:4)
REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates
REAL(8):: invJ(1:2,1:2), detJ
@ -498,9 +497,9 @@ MODULE moduleMesh2DCart
self%n3%emData%phi, &
self%n4%emData%phi /)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
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)
@ -508,11 +507,11 @@ MODULE moduleMesh2DCart
END FUNCTION gatherEFQuad
PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF)
PURE FUNCTION gatherMFQuad(self,Xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
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)
@ -530,7 +529,7 @@ MODULE moduleMesh2DCart
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(Xi, fPsi)
MF = MATMUL(fPsi(:), MF_Nodes)
END FUNCTION gatherMFQuad
@ -548,47 +547,47 @@ MODULE moduleMesh2DCart
END FUNCTION getNodesQuad
!Transforms physical coordinates to element coordinates
PURE FUNCTION phy2logQuad(self,r) RESULT(xN)
PURE FUNCTION phy2logQuad(self,r) RESULT(XiN)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2)
REAL(8):: XiN(1:3)
REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2)
REAL(8):: dPsi(1:2,1:4), fPsi(1:4)
REAL(8):: conv
!Iterative newton method to transform coordinates
conv=1.D0
xO=0.D0
XiO=0.D0
DO WHILE(conv>1.D-4)
dPsi = self%dPsi(xO)
invJ = self%invJac(xO, dPsi)
fPsi = self%fPsi(xO)
dPsi = self%dPsi(XiO)
invJ = self%invJac(XiO, dPsi)
CALL self%fPsi(XiO, fPsi)
f(1) = DOT_PRODUCT(fPsi,self%x)-r(1)
f(2) = DOT_PRODUCT(fPsi,self%y)-r(2)
detJ = self%detJac(xO,dPsi)
xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(xN-xO),1)
xO=xN
detJ = self%detJac(XiO,dPsi)
XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(XiN-XiO),1)
XiO=XiN
END DO
END FUNCTION phy2logQuad
!Gets the next element for a logical position xi
SUBROUTINE nextElementQuad(self, xi, nextElement)
!Gets the next element for a logical position Xi
SUBROUTINE nextElementQuad(self, Xi, nextElement)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:4)
REAL(8):: XiArray(1:4)
INTEGER:: nextInt
xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /)
nextInt = MAXLOC(xiArray,1)
XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
nextInt = MAXLOC(XiArray,1)
!Selects the higher value of directions and searches in that direction
NULLIFY(nextElement)
SELECT CASE (nextInt)
@ -649,14 +648,14 @@ MODULE moduleMesh2DCart
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3)
xii(1) = random( 0.D0, 1.D0)
xii(2) = random( 0.D0, 1.D0 - xii(1))
xii(3) = 0.D0
Xi(1) = random( 0.D0, 1.D0)
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
Xi(3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(Xi, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
@ -669,55 +668,53 @@ MODULE moduleMesh2DCart
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(inout):: self
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
detJ = self%detJac(xi)/2.D0
fPsi = self%fPsi(xi)
Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
detJ = self%detJac(Xi)/2.D0
CALL self%fPsi(Xi, fPsi)
self%volume = detJ
self%arNodes = fPsi*detJ
END SUBROUTINE areaTria
!Shape functions for triangular element
PURE FUNCTION fPsiTria(xi) RESULT(fPsi)
PURE SUBROUTINE fPsiTria(Xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(out):: fPsi(:)
ALLOCATE(fPsi(1:3))
fPsi(1) = 1.D0 - Xi(1) - Xi(2)
fPsi(2) = Xi(1)
fPsi(3) = Xi(2)
fPsi(1) = 1.D0 - xi(1) - xi(2)
fPsi(2) = xi(1)
fPsi(3) = xi(2)
END SUBROUTINE fPsiTria
END FUNCTION fPsiTria
!Derivative element function at coordinates xi
PURE FUNCTION dPsiTria(xi) RESULT(dPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiTria(Xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:2,1:3))
dPsi(1,:) = dPsiTriaXi1(xi(2))
dPsi(2,:) = dPsiTriaXi2(xi(1))
dPsi(1,:) = dPsiTriaXi1(Xi(2))
dPsi(2,:) = dPsiTriaXi2(Xi(1))
END FUNCTION dPsiTria
!Derivative element function (xi1)
PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1)
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8), INTENT(in):: xi2
REAL(8), INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:3)
dPsiXi1(1) = -1.D0
@ -726,11 +723,11 @@ MODULE moduleMesh2DCart
END FUNCTION dPsiTriaXi1
!Derivative element function (xi1)
PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2)
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8), INTENT(in):: xi1
REAL(8), INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:3)
dPsiXi2(1) = -1.D0
@ -759,22 +756,22 @@ MODULE moduleMesh2DCart
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
REAL(8):: invJ(1:2,1:2), detJ
INTEGER:: l
ALLOCATE(localK(1:4, 1:4))
localK=0.D0
xi=0.D0
Xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 4
xi(1) = xi1Tria(l)
xi(2) = xi2Tria(l)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
fPsi = self%fPsi(xi)
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
CALL self%fPsi(Xi, fPsi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
END DO
@ -789,19 +786,19 @@ MODULE moduleMesh2DCart
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:3)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
INTEGER:: l
ALLOCATE(localF(1:3))
localF = 0.D0
xi = 0.D0
Xi = 0.D0
!Start 2D Gauss Quad Integral
DO l=1, 4
xi(1) = xi1Tria(l)
xi(2) = xi2Tria(l)
detJ = self%detJac(xi)
fPsi = self%fPsi(xi)
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
detJ = self%detJac(Xi)
CALL self%fPsi(Xi, fPsi)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wTria(l)*detJ
@ -809,24 +806,24 @@ MODULE moduleMesh2DCart
END FUNCTION elemFTria
PURE FUNCTION insideTria(xi) RESULT(ins)
PURE FUNCTION insideTria(Xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins
ins = xi(1) >= 0.D0 .AND. &
xi(2) >= 0.D0 .AND. &
1.D0 - xi(1) - xi(2) >= 0.D0
ins = Xi(1) >= 0.D0 .AND. &
Xi(2) >= 0.D0 .AND. &
1.D0 - Xi(1) - Xi(2) >= 0.D0
END FUNCTION insideTria
!Gathers the electric field at position xi
PURE FUNCTION gatherEFTria(self,xi) RESULT(EF)
!Gathers the electric field at position Xi
PURE FUNCTION gatherEFTria(self,Xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: dPsi(1:2,1:3)
REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates
REAL(8):: invJ(1:2,1:2), detJ
@ -837,9 +834,9 @@ MODULE moduleMesh2DCart
self%n2%emData%phi, &
self%n3%emData%phi /)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
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)
@ -847,11 +844,11 @@ MODULE moduleMesh2DCart
END FUNCTION gatherEFTria
PURE FUNCTION gatherMFTria(self,xi) RESULT(MF)
PURE FUNCTION gatherMFTria(self,Xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: fPsi(1:3)
REAL(8):: MF_Nodes(1:3,1:3)
REAL(8):: MF(1:3)
@ -866,7 +863,7 @@ MODULE moduleMesh2DCart
self%n2%emData%B(3), &
self%n3%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(Xi, fPsi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTria
@ -884,37 +881,37 @@ MODULE moduleMesh2DCart
END FUNCTION getNodesTria
!Transforms physical coordinates to element coordinates
PURE FUNCTION phy2logTria(self,r) RESULT(xi)
PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: deltaR(1:2)
REAL(8):: dPsi(1:2,1:3)
!Direct method to convert coordinates
xi = 0.D0 !Irrelevant, required for input
Xi = 0.D0 !Irrelevant, required for input
deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)
detJ = self%detJac(xi, dPsi)
xi(1:2) = MATMUL(invJ,deltaR)/detJ
dPsi = self%dPsi(Xi)
invJ = self%invJac(Xi, dPsi)
detJ = self%detJac(Xi, dPsi)
Xi(1:2) = MATMUL(invJ,deltaR)/detJ
END FUNCTION phy2logTria
SUBROUTINE nextElementTria(self, xi, nextElement)
SUBROUTINE nextElementTria(self, Xi, nextElement)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:3)
REAL(8):: XiArray(1:3)
INTEGER:: nextInt
xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /)
nextInt = MINLOC(xiArray,1)
XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /)
nextInt = MINLOC(XiArray,1)
NULLIFY(nextElement)
SELECT CASE (nextInt)
CASE (1)
@ -929,11 +926,11 @@ MODULE moduleMesh2DCart
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D
!Computes element Jacobian determinant
PURE FUNCTION detJ2DCart(self, xi, dPsi_in) RESULT(dJ)
PURE FUNCTION detJ2DCart(self, Xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol2DCart), INTENT(in):: self
REAL(8), INTENT(in):: xi(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
@ -943,7 +940,7 @@ MODULE moduleMesh2DCart
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
dPsi = self%dPsi(Xi)
END IF
@ -953,11 +950,11 @@ MODULE moduleMesh2DCart
END FUNCTION detJ2DCart
!Computes element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ2DCart(self,xi,dPsi_in) RESULT(invJ)
PURE FUNCTION invJ2DCart(self,Xi,dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol2DCart), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1:2), dy(1:2)
@ -967,7 +964,7 @@ MODULE moduleMesh2DCart
dPsi=dPsi_in
ELSE
dPsi = self%dPsi(xi)
dPsi = self%dPsi(Xi)
END IF

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

@ -310,49 +310,52 @@ MODULE moduleMesh2DCyl
CLASS(meshVol2DCylQuad), INTENT(inout):: self
REAL(8):: r, xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:4)
REAL(8):: fPsi(1:4), fPsi_node(1:4)
self%volume = 0.D0
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
xi = 0.D0
detJ = self%detJac(xi)*PI8 !4*2*pi
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
!Computes volume per node
xi = (/-5.D-1, -5.D-1, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
CALL self%fPsi(xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,self%r)
self%arNodes(1) = fPsi(1)*r*detJ
xi = (/ 5.D-1, -5.D-1, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
CALL self%fPsi(xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,self%r)
self%arNodes(2) = fPsi(2)*r*detJ
xi = (/ 5.D-1, 5.D-1, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
CALL self%fPsi(xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,self%r)
self%arNodes(3) = fPsi(3)*r*detJ
xi = (/-5.D-1, 5.D-1, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
CALL self%fPsi(xi, fPsi_node)
r = DOT_PRODUCT(fPsi_node,self%r)
self%arNodes(4) = fPsi(4)*r*detJ
END SUBROUTINE areaQuad
!Computes element functions in point xi
PURE FUNCTION fPsiQuad(xi) RESULT(fPsi)
PURE SUBROUTINE fPsiQuad(xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:4))
REAL(8), INTENT(out):: fPsi(:)
fPsi(1) = (1.D0-xi(1)) * (1.D0-xi(2))
fPsi(2) = (1.D0+xi(1)) * (1.D0-xi(2))
fPsi(3) = (1.D0+xi(1)) * (1.D0+xi(2))
fPsi(4) = (1.D0-xi(1)) * (1.D0+xi(2))
fPsi = fPsi*0.25D0
END FUNCTION fPsiQuad
END SUBROUTINE fPsiQuad
!Derivative element function at coordinates xi
PURE FUNCTION dPsiQuad(xi) RESULT(dPsi)
@ -379,6 +382,7 @@ MODULE moduleMesh2DCyl
dPsiXi1(2) = (1.D0 - xi2)
dPsiXi1(3) = (1.D0 + xi2)
dPsiXi1(4) = -(1.D0 + xi2)
dPsiXi1 = dPsiXi1*0.25D0
END FUNCTION dPsiQuadXi1
@ -394,6 +398,7 @@ MODULE moduleMesh2DCyl
dPsiXi2(2) = -(1.D0 + xi1)
dPsiXi2(3) = (1.D0 + xi1)
dPsiXi2(4) = (1.D0 - xi1)
dPsiXi2 = dPsiXi2 * 0.25D0
END FUNCTION dPsiQuadXi2
@ -427,7 +432,7 @@ MODULE moduleMesh2DCyl
xii(2) = random(-1.D0, 1.D0)
xii(3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(xii, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%z)
r(2) = DOT_PRODUCT(fPsi, self%r)
@ -457,7 +462,7 @@ MODULE moduleMesh2DCyl
DO m = 1, 3
xi(1) = corQuad(m)
dPsi(2,:) = self%dPsiXi2(xi(1))
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
r = DOT_PRODUCT(fPsi,self%r)
@ -492,7 +497,7 @@ MODULE moduleMesh2DCyl
DO m = 1, 3
xi(2) = corQuad(m)
detJ = self%detJac(xi)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
r = DOT_PRODUCT(fPsi,self%r)
f = DOT_PRODUCT(fPsi,source)
localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ
@ -564,7 +569,7 @@ MODULE moduleMesh2DCyl
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
MF = MATMUL(fPsi(:), MF_Nodes)
END FUNCTION gatherMFQuad
@ -582,30 +587,30 @@ MODULE moduleMesh2DCyl
END FUNCTION getNodesQuad
!Transforms physical coordinates to element coordinates
PURE FUNCTION phy2logQuad(self,r) RESULT(xN)
PURE FUNCTION phy2logQuad(self,r) RESULT(XiN)
IMPLICIT NONE
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2)
REAL(8):: XiN(1:3)
REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2)
REAL(8):: dPsi(1:2,1:4), fPsi(1:4)
REAL(8):: conv
!Iterative newton method to transform coordinates
conv=1.D0
xO=0.D0
XiO=0.D0
DO WHILE(conv>1.D-4)
dPsi = self%dPsi(xO)
invJ = self%invJac(xO, dPsi)
fPsi = self%fPsi(xO)
f(1) = DOT_PRODUCT(fPsi,self%z)-r(1)
f(2) = DOT_PRODUCT(fPsi,self%r)-r(2)
detJ = self%detJac(xO,dPsi)
xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(xN-xO),1)
xO=xN
DO WHILE(conv>1.D-3)
CALL self%fPsi(XiO, fPsi)
f = (/ DOT_PRODUCT(fPsi,self%z)-r(1), &
DOT_PRODUCT(fPsi,self%r)-r(2) /)
dPsi = self%dPsi(XiO)
invJ = self%invJac(XiO, dPsi)
detJ = self%detJac(XiO,dPsi)
XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(XiN-XiO),1)
XiO=XiN
END DO
@ -690,7 +695,7 @@ MODULE moduleMesh2DCyl
xii(2) = random( 0.D0, 1.D0 - xii(1))
xii(3) = 0.D0
fPsi = self%fPsi(xii)
CALL self%fPsi(xii, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%z)
r(2) = DOT_PRODUCT(fPsi, self%r)
@ -713,7 +718,7 @@ MODULE moduleMesh2DCyl
!2D 1 point Gauss Quad Integral
xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
detJ = self%detJac(xi)*PI !2PI*1/2
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
@ -723,19 +728,17 @@ MODULE moduleMesh2DCyl
END SUBROUTINE areaTria
!Shape functions for triangular element
PURE FUNCTION fPsiTria(xi) RESULT(fPsi)
PURE SUBROUTINE fPsiTria(xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:3))
REAL(8), INTENT(out):: fPsi(:)
fPsi(1) = 1.D0 - xi(1) - xi(2)
fPsi(2) = xi(1)
fPsi(3) = xi(2)
END FUNCTION fPsiTria
END SUBROUTINE fPsiTria
!Derivative element function at coordinates xi
PURE FUNCTION dPsiTria(xi) RESULT(dPsi)
@ -813,7 +816,7 @@ MODULE moduleMesh2DCyl
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
r = DOT_PRODUCT(fPsi,self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
@ -843,7 +846,7 @@ MODULE moduleMesh2DCyl
xi(1) = xi1Tria(l)
xi(2) = xi2Tria(l)
detJ = self%detJac(xi)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
r = DOT_PRODUCT(fPsi,self%r)
f = DOT_PRODUCT(fPsi,source)
localF = localF + r*f*fPsi*wTria(l)*detJ
@ -910,7 +913,7 @@ MODULE moduleMesh2DCyl
self%n2%emData%B(3), &
self%n3%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(xi, fPsi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTria

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

@ -41,8 +41,8 @@ MODULE moduleMesh3DCart
END TYPE meshVol3DCart
ABSTRACT INTERFACE
PURE FUNCTION dPsi_interface(xii) RESULT(dPsi)
REAL(8), INTENT(in):: xii(1:3)
PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
END FUNCTION dPsi_interface
@ -71,8 +71,8 @@ MODULE moduleMesh3DCart
PROCEDURE, PASS:: calcVol => volumeTetra
PROCEDURE, NOPASS:: fPsi => fPsiTetra
PROCEDURE, NOPASS:: dPsi => dPsiTetra
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXii1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXii2
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXi1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXi2
PROCEDURE, PASS:: partialDer => partialDerTetra
PROCEDURE, PASS:: elemK => elemKTetra
PROCEDURE, PASS:: elemF => elemFTetra
@ -213,14 +213,14 @@ MODULE moduleMesh3DCart
CLASS(meshEdge3DCartTria), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3)
xii(1) = random( 0.D0, 1.D0)
xii(2) = random( 0.D0, 1.D0 - xii(1))
xii(3) = 0.D0
Xi(1) = random( 0.D0, 1.D0)
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
Xi(3) = 0.D0
fPsi = self%fPsi(xii)
fPsi = self%fPsi(Xi)
r = (/DOT_PRODUCT(fPsi, self%x), &
DOT_PRODUCT(fPsi, self%y), &
DOT_PRODUCT(fPsi, self%z)/)
@ -228,17 +228,17 @@ MODULE moduleMesh3DCart
END FUNCTION randPosEdgeTria
!Shape functions for triangular surface
PURE FUNCTION fPsiEdgeTria(xii) RESULT(fPsi)
PURE FUNCTION fPsiEdgeTria(Xi) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:3))
fPsi(1) = 1.D0 - xii(1) - xii(2)
fPsi(2) = xii(1)
fPsi(3) = xii(2)
fPsi(1) = 1.D0 - Xi(1) - Xi(2)
fPsi(2) = Xi(1)
fPsi(3) = Xi(2)
END FUNCTION fPsiEdgeTria
@ -254,6 +254,7 @@ MODULE moduleMesh3DCart
INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node
REAL(8):: Xi(1:3), fPsi(1:4)
REAL(8):: volNodes(1:4) !Volume of each node
self%n = n
@ -274,8 +275,9 @@ MODULE moduleMesh3DCart
CALL self%calcVol()
!Assign proportional volume to each node
!TODO: Review this to apply to all elements in the future
volNodes = self%fPsi((/0.25D0, 0.25D0, 0.25D0/))*self%volume
Xi = (/0.25D0, 0.25D0, 0.25D0/)
CALL self%fPsi(Xi, fPsi)
volNodes = fPsi*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)
@ -295,19 +297,18 @@ MODULE moduleMesh3DCart
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
xii(1) = random( 0.D0, 1.D0)
xii(2) = random( 0.D0, 1.D0 - xii(1))
xii(3) = random( 0.D0, 1.D0 - xii(1) - xii(2))
Xi(1) = random( 0.D0, 1.D0)
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2))
ALLOCATE(fPsi(1:4))
fPsi = self%fPsi(xii)
CALL self%fPsi(Xi, fPsi)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = DOT_PRODUCT(fPsi, self%z)
r = (/ DOT_PRODUCT(fPsi, self%x), &
DOT_PRODUCT(fPsi, self%y), &
DOT_PRODUCT(fPsi, self%z) /)
END FUNCTION randPosVolTetra
@ -316,83 +317,81 @@ MODULE moduleMesh3DCart
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout):: self
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
self%volume = 0.D0
xii = (/0.25D0, 0.25D0, 0.25D0/)
self%volume = self%detJac(xii)
Xi = (/0.25D0, 0.25D0, 0.25D0/)
self%volume = self%detJac(Xi)
END SUBROUTINE volumeTetra
!Computes element functions in point xii
PURE FUNCTION fPsiTetra(xi) RESULT(fPsi)
!Computes element functions in point Xi
PURE SUBROUTINE fPsiTetra(Xi, fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(out):: fPsi(:)
ALLOCATE(fPsi(1:4))
fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3)
fPsi(2) = Xi(1)
fPsi(3) = Xi(2)
fPsi(4) = Xi(3)
fPsi(1) = 1.D0 - xi(1) - xi(2) - xi(3)
fPsi(2) = xi(1)
fPsi(3) = xi(2)
fPsi(4) = xi(3)
END SUBROUTINE fPsiTetra
END FUNCTION fPsiTetra
!Derivative element function at coordinates xii
PURE FUNCTION dPsiTetra(xii) RESULT(dPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiTetra(Xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:3,1:4))
dPsi(1,:) = dPsiTetraXii1(xii(2), xii(3))
dPsi(2,:) = dPsiTetraXii2(xii(1), xii(3))
dPsi(3,:) = dPsiTetraXii3(xii(1), xii(2))
dPsi(1,:) = dPsiTetraXi1(Xi(2), Xi(3))
dPsi(2,:) = dPsiTetraXi2(Xi(1), Xi(3))
dPsi(3,:) = dPsiTetraXi3(Xi(1), Xi(2))
END FUNCTION dPsiTetra
!Derivative element function respect to xii1
PURE FUNCTION dPsiTetraXii1(xii2, xii3) RESULT(dPsiXii1)
!Derivative element function respect to Xi1
PURE FUNCTION dPsiTetraXi1(Xi2, Xi3) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8), INTENT(in):: xii2, xii3
REAL(8):: dPsiXii1(1:4)
REAL(8), INTENT(in):: Xi2, Xi3
REAL(8):: dPsiXi1(1:4)
dPsiXii1(1) = -1.D0
dPsiXii1(2) = 1.D0
dPsiXii1(3) = 0.D0
dPsiXii1(4) = 0.D0
dPsiXi1(1) = -1.D0
dPsiXi1(2) = 1.D0
dPsiXi1(3) = 0.D0
dPsiXi1(4) = 0.D0
END FUNCTION dPsiTetraXii1
END FUNCTION dPsiTetraXi1
!Derivative element function respect to xii2
PURE FUNCTION dPsiTetraXii2(xii1, xii3) RESULT(dPsiXii2)
!Derivative element function respect to Xi2
PURE FUNCTION dPsiTetraXi2(Xi1, Xi3) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8), INTENT(in):: xii1, xii3
REAL(8):: dPsiXii2(1:4)
REAL(8), INTENT(in):: Xi1, Xi3
REAL(8):: dPsiXi2(1:4)
dPsiXii2(1) = -1.D0
dPsiXii2(2) = 0.D0
dPsiXii2(3) = 1.D0
dPsiXii2(4) = 0.D0
dPsiXi2(1) = -1.D0
dPsiXi2(2) = 0.D0
dPsiXi2(3) = 1.D0
dPsiXi2(4) = 0.D0
END FUNCTION dPsiTetraXii2
END FUNCTION dPsiTetraXi2
!Derivative element function respect to xii3
PURE FUNCTION dPsiTetraXii3(xii1, xii2) RESULT(dPsiXii3)
!Derivative element function respect to Xi3
PURE FUNCTION dPsiTetraXi3(Xi1, Xi2) RESULT(dPsiXi3)
IMPLICIT NONE
REAL(8), INTENT(in):: xii1, xii2
REAL(8):: dPsiXii3(1:4)
REAL(8), INTENT(in):: Xi1, Xi2
REAL(8):: dPsiXi3(1:4)
dPsiXii3(1) = -1.D0
dPsiXii3(2) = 0.D0
dPsiXii3(3) = 0.D0
dPsiXii3(4) = 1.D0
dPsiXi3(1) = -1.D0
dPsiXi3(2) = 0.D0
dPsiXi3(3) = 0.D0
dPsiXi3(4) = 1.D0
END FUNCTION dPsiTetraXii3
END FUNCTION dPsiTetraXi3
!Computes the derivatives in global coordinates
PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz)
@ -421,19 +420,19 @@ MODULE moduleMesh3DCart
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: invJ(1:3,1:3), detJ
ALLOCATE(localK(1:4,1:4))
localK = 0.D0
xii = 0.D0
Xi = 0.D0
!TODO: One point Gauss integral. Upgrade when possible
xii = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(xii)
detJ = self%detJac(xii, dPsi)
invJ = self%invJac(xii, dPsi)
fPsi = self%fPsi(xii)
Xi = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi, dPsi)
invJ = self%invJac(Xi, dPsi)
CALL self%fPsi(Xi, fPsi)
localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
END FUNCTION elemKTetra
@ -445,40 +444,40 @@ MODULE moduleMesh3DCart
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: xii(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
ALLOCATE(localF(1:4))
localF = 0.D0
xii = 0.D0
!TODO: One point Gauss integral. Upgrade when possible
xii = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(xii)
detJ = self%detJac(xii, dPsi)
fPsi = self%fPsi(xii)
Xi = 0.D0
Xi = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi, dPsi)
CALL self%fPsi(Xi, fPsi)
f = DOT_PRODUCT(fPsi, source)
localF = f*fPsi*1.D0*detJ
END FUNCTION elemFTetra
PURE FUNCTION insideTetra(xi) RESULT(ins)
PURE FUNCTION insideTetra(Xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
LOGICAL:: ins
ins = xi(1) >= 0.D0 .AND. &
xi(2) >= 0.D0 .AND. &
xi(3) >= 0.D0 .AND. &
1.D0 - xi(1) - xi(2) - xi(3) >= 0.D0
ins = Xi(1) >= 0.D0 .AND. &
Xi(2) >= 0.D0 .AND. &
Xi(3) >= 0.D0 .AND. &
1.D0 - Xi(1) - Xi(2) - Xi(3) >= 0.D0
END FUNCTION insideTetra
PURE FUNCTION gatherEFTetra(self, xi) RESULT(EF)
PURE FUNCTION gatherEFTetra(self, Xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
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
@ -490,9 +489,9 @@ MODULE moduleMesh3DCart
self%n3%emData%phi, &
self%n4%emData%phi /)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi, dPsi)
invJ = self%invJac(xi, dPsi)
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)
@ -500,11 +499,11 @@ MODULE moduleMesh3DCart
END FUNCTION gatherEFTetra
PURE FUNCTION gatherMFTetra(self, xi) RESULT(MF)
PURE FUNCTION gatherMFTetra(self, Xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
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)
@ -522,7 +521,7 @@ MODULE moduleMesh3DCart
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
fPsi = self%fPsi(xi)
CALL self%fPsi(Xi, fPsi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFTetra
@ -538,37 +537,37 @@ MODULE moduleMesh3DCart
END FUNCTION getNodesTetra
PURE FUNCTION phy2logTetra(self,r) RESULT(xi)
PURE FUNCTION phy2logTetra(self,r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: invJ(1:3, 1:3), detJ
REAL(8):: deltaR(1:3)
REAL(8):: dPsi(1:3, 1:4)
xi = 0.D0
Xi = 0.D0
deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)
detJ = self%detJac(xi, dPsi)
xi = MATMUL(invJ, deltaR)/detJ
dPsi = self%dPsi(Xi)
invJ = self%invJac(Xi, dPsi)
detJ = self%detJac(Xi, dPsi)
Xi = MATMUL(invJ, deltaR)/detJ
END FUNCTION phy2logTetra
SUBROUTINE nextElementTetra(self, xi, nextElement)
SUBROUTINE nextElementTetra(self, Xi, nextElement)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:4)
REAL(8):: XiArray(1:4)
INTEGER:: nextInt
!TODO: Review when connectivity
xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /)
nextInt = MINLOC(xiArray, 1)
XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), Xi(2), Xi(1) /)
nextInt = MINLOC(XiArray, 1)
NULLIFY(nextElement)
SELECT CASE(nextInt)
CASE (1)
@ -585,11 +584,11 @@ MODULE moduleMesh3DCart
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D
!Computes element Jacobian determinant
PURE FUNCTION detJ3DCart(self, xi, dPsi_in) RESULT(dJ)
PURE FUNCTION detJ3DCart(self, Xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol3DCart), INTENT(in)::self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:, 1:)
REAL(8):: dJ
REAL(8), ALLOCATABLE:: dPsi(:,:)
@ -599,7 +598,7 @@ MODULE moduleMesh3DCart
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
dPsi = self%dPsi(Xi)
END IF
@ -610,11 +609,11 @@ MODULE moduleMesh3DCart
END FUNCTION detJ3DCart
PURE FUNCTION invJ3DCart(self,xi,dPsi_in) RESULT(invJ)
PURE FUNCTION invJ3DCart(self,Xi,dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol3DCart), INTENT(in):: self
REAL(8), INTENT(in):: xi(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
@ -624,7 +623,7 @@ MODULE moduleMesh3DCart
dPsi=dPsi_in
ELSE
dPsi = self%dPsi(xi)
dPsi = self%dPsi(Xi)
END IF

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

@ -211,11 +211,11 @@ MODULE moduleMesh
END FUNCTION getNodesVol_interface
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
PURE SUBROUTINE fPsi_interface(xi, fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8), INTENT(out):: fPsi(:)
END FUNCTION fPsi_interface
END SUBROUTINE fPsi_interface
PURE FUNCTION elemK_interface(self) RESULT(localK)
IMPORT:: meshVol
@ -496,11 +496,14 @@ MODULE moduleMesh
INTEGER:: i, nNodes
CLASS(meshNode), POINTER:: node
fPsi = self%fPsi(part%xi)
tensorS = outerProduct(part%v, part%v)
sp = part%species%n
volNodes = self%getNodes()
nNodes = SIZE(volNodes)
ALLOCATE(fPsi(1:nNodes))
CALL self%fPsi(part%xi, fPsi)
tensorS = outerProduct(part%v, part%v)
sp = part%species%n
DO i = 1, nNodes
node => mesh%nodes(volNodes(i))%obj