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

Small improvement

Browse source 
Very small improvement in performance.

Still, partialDer takes too long to compute.
Trying to find ways to improve it.
This commit is contained in:
Jorge Gonzalez 2023-01-06 15:18:04 +01:00
commit 7b7a5c45ca
6 changed files with 788 additions and 837 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -6,6 +6,7 @@ MODULE moduleMesh0D
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D
INTEGER:: n1
CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode0D
PROCEDURE, PASS:: getCoordinates => getCoord0D
@ -17,17 +18,18 @@ MODULE moduleMesh0D
PROCEDURE, PASS:: init => initCell0D
PROCEDURE, PASS:: getNodes => getNodes0D
PROCEDURE, PASS:: randPos => randPos0D
PROCEDURE, PASS:: fPsi => fPsi0D
PROCEDURE, PASS:: dPsi => dPsi0D
PROCEDURE, PASS:: detJac => detJ0D
PROCEDURE, PASS:: invJac => invJ0D
PROCEDURE, PASS:: elemK => elemK0D
PROCEDURE, PASS:: elemF => elemF0D
PROCEDURE, NOPASS:: fPsi => fPsi0D
PROCEDURE, NOPASS:: dPsi => dPsi0D
PROCEDURE, PASS:: partialDer => partialDer0D
PROCEDURE, NOPASS:: detJac => detJ0D
PROCEDURE, NOPASS:: invJac => invJ0D
PROCEDURE, PASS:: gatherElectricField => gatherEF0D
PROCEDURE, PASS:: gatherMagneticField => gatherMF0D
PROCEDURE, PASS:: elemK => elemK0D
PROCEDURE, PASS:: elemF => elemF0D
PROCEDURE, PASS:: phy2log => phy2log0D
PROCEDURE, NOPASS:: inside => inside0D
PROCEDURE, PASS:: nextElement => nextElement0D
PROCEDURE, PASS:: neighbourElement => neighbourElement0D
END TYPE meshCell0D
@ -89,6 +91,7 @@ MODULE moduleMesh0D
END SUBROUTINE initCell0D
!Get the nodes of the volume
PURE FUNCTION getNodes0D(self, nNodes) RESULT(n)
IMPLICIT NONE
@ -100,6 +103,7 @@ MODULE moduleMesh0D
END FUNCTION getNodes0D
!Calculate random position inside the volume
FUNCTION randPos0D(self) RESULT(r)
IMPLICIT NONE
@ -110,10 +114,9 @@ MODULE moduleMesh0D
END FUNCTION randPos0D
PURE FUNCTION fPsi0D(self, Xi, nNodes) RESULT(fPsi)
PURE FUNCTION fPsi0D(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
@ -122,10 +125,9 @@ MODULE moduleMesh0D
END FUNCTION fPsi0D
PURE FUNCTION dPsi0D(self, Xi, nNodes) RESULT(dPsi)
PURE FUNCTION dPsi0D(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
@ -134,31 +136,17 @@ MODULE moduleMesh0D
END FUNCTION dPsi0D
PURE FUNCTION detJ0D(self, Xi, nNodes, dPsi_in) RESULT(dJ)
PURE FUNCTION partialDer0D(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE
CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8):: dJ
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8):: pDer(1:3, 1:3)
dJ = 0.D0
pDer = 0.D0
END FUNCTION detJ0D
PURE FUNCTION invJ0D(self, Xi, nNodes, dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8):: invJ(1:3,1:3)
invJ = 0.D0
END FUNCTION invJ0D
END FUNCTION partialDer0D
PURE FUNCTION elemK0D(self, nNodes) RESULT(localK)
IMPLICIT NONE
@ -234,15 +222,36 @@ MODULE moduleMesh0D
END FUNCTION inside0D
SUBROUTINE nextElement0D(self, Xi, nextElement)
SUBROUTINE neighbourElement0D(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell0D), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
nextElement => NULL()
neighbourElement => NULL()
END SUBROUTINE nextElement0D
END SUBROUTINE neighbourElement0D
!COMMON FUNCTIONS
PURE FUNCTION detJ0D(pDer) RESULT(dJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: dJ
dJ = 0.D0
END FUNCTION detJ0D
PURE FUNCTION invJ0D(pDer) RESULT(invJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3)
invJ = 0.D0
END FUNCTION invJ0D
END MODULE moduleMesh0D

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

@ -14,6 +14,7 @@ MODULE moduleMesh1DCart
!Element coordinates
REAL(8):: x = 0.D0
CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode1DCart
PROCEDURE, PASS:: getCoordinates => getCoord1DCart
@ -25,6 +26,7 @@ MODULE moduleMesh1DCart
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge1DCart
PROCEDURE, PASS:: getNodes => getNodes1DCart
PROCEDURE, PASS:: intersection => intersection1DCart
@ -32,27 +34,7 @@ MODULE moduleMesh1DCart
END TYPE meshEdge1DCart
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DCart
CONTAINS
PROCEDURE, PASS:: detJac => detJ1DCart
PROCEDURE, PASS:: invJac => invJ1DCart
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshCell1DCart
ABSTRACT INTERFACE
PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx)
IMPORT meshCell1DCart
CLASS(meshCell1DCart), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1):: dx
END SUBROUTINE partialDer_interface
END INTERFACE
TYPE, PUBLIC, EXTENDS(meshCell1DCart):: meshCell1DCartSegm
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DCartSegm
!Element coordinates
REAL(8):: x(1:2)
!Connectivity to nodes
@ -61,20 +43,24 @@ MODULE moduleMesh1DCart
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCell1DCartSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: randPos => randPos1DCartSegm
PROCEDURE, PASS:: area => areaSegm
PROCEDURE, PASS:: fPsi => fPsiSegm
PROCEDURE, PASS:: dPsi => dPsiSegm
PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: partialDer => partialDerSegm
PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: detJac => detJ1DCart
PROCEDURE, NOPASS:: invJac => invJ1DCart
PROCEDURE, PASS:: gatherElectricField => gatherEFSegm
PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm
PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: nextElement => nextElementSegm
PROCEDURE, PASS:: neighbourElement => neighbourElementSegm
!PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: area => areaSegm
END TYPE meshCell1DCartSegm
@ -219,7 +205,19 @@ MODULE moduleMesh1DCart
END SUBROUTINE initCell1DCartSegm
!Calculates a random position in 1D volume
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesSegm
!Random position in 1D volume
FUNCTION randPos1DCartSegm(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
@ -227,135 +225,63 @@ MODULE moduleMesh1DCart
CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
REAL(8):: fPsi(1:2)
Xi = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2:3) = 0.D0
fPsi = self%fPsi(Xi, 2)
r = 0.D0
r(1) = DOT_PRODUCT(fPsi, self%x)
END FUNCTION randPos1DCartSegm
!Computes element area
PURE SUBROUTINE areaSegm(self)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xi(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!1 point Gauss integral
Xi = 0.D0
fPsi = self%fPsi(Xi, 2)
detJ = self%detJac(Xi, 2)
l = 2.D0*detJ
self%volume = l
self%arNodes = fPsi*l
END SUBROUTINE areaSegm
!Computes element functions at point Xi
PURE FUNCTION fPsiSegm(self, Xi, nNodes) RESULT(fPsi)
PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = 1.D0 - Xi(1)
fPsi(2) = 1.D0 + Xi(1)
fPsi = (/ 1.D0 - Xi(1), &
1.D0 + Xi(1) /)
fPsi = fPsi * 5.D-1
fPsi = fPsi * 0.50D0
END FUNCTION fPsiSegm
!Computes element derivative shape function at Xi
PURE FUNCTION dPsiSegm(self, Xi, nNodes) RESULT(dPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
dPsi(1, 1) = -5.D-1
dPsi(1, 2) = 5.D-1
dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
END FUNCTION dPsiSegm
!Computes partial derivatives of coordinates
PURE SUBROUTINE partialDerSegm(self, nNodes, dPsi, dx)
!Partial derivative in global coordinates
PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1):: dx
REAL(8):: pDer(1:3, 1:3)
dx(1) = DOT_PRODUCT(dPsi(1,:), self%x)
pDer = 0.D0
END SUBROUTINE partialDerSegm
pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%x(1:2))
pDer(2,2) = 1.D0
pDer(3,3) = 1.D0
!Computes local stiffness matrix
PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
localK = 0.D0
Xi = 0.D0
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
detJ = self%detJac(Xi, 2, dPsi)
invJ = self%invJac(Xi, 2, dPsi)
localK = localK + 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, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:2)
REAL(8):: detJ, f
REAL(8):: Xi(1:3)
INTEGER:: l
localF = 0.D0
Xi = 0.D0
DO l = 1, 3
Xi(1) = corSeg(l)
detJ = self%detJac(Xi, 2)
fPsi = self%fPsi(Xi, 2)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*wSeg(l)*detJ
END DO
END FUNCTION elemFSegm
END FUNCTION partialDerSegm
PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array)
IMPLICIT NONE
@ -391,6 +317,68 @@ MODULE moduleMesh1DCart
END FUNCTION gatherMFSegm
!Computes element local stiffness matrix
PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
localK = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
wSeg(l)/detJ
END DO
END FUNCTION elemKSegm
!Computes the local source vector for a force f
PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:2)
REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
INTEGER:: l
localF = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 2)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*wSeg(l)*detJ
END DO
END FUNCTION elemFSegm
PURE FUNCTION insideSegm(Xi) RESULT(ins)
IMPLICIT NONE
@ -402,101 +390,87 @@ MODULE moduleMesh1DCart
END FUNCTION insideSegm
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesSegm
PURE FUNCTION phy2logSegm(self, r) RESULT(xN)
PURE FUNCTION phy2logSegm(self, r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
REAL(8):: Xi(1:3)
xN = 0.D0
xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
Xi = 0.D0
Xi(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
END FUNCTION phy2logSegm
!Get next element for a logical position Xi
SUBROUTINE nextElementSegm(self, Xi, nextElement)
!Get the next element for a logical position Xi
SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell1DCartSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
NULLIFY(nextElement)
NULLIFY(neighbourElement)
IF (Xi(1) < -1.D0) THEN
nextElement => self%e2
neighbourElement => self%e2
ELSEIF (Xi(1) > 1.D0) THEN
nextElement => self%e1
neighbourElement => self%e1
END IF
END SUBROUTINE nextElementSegm
END SUBROUTINE neighbourElementSegm
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Calculates a random position in 1D volume
!Computes the element Jacobian determinant
PURE FUNCTION detJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ)
!Computes element area
PURE SUBROUTINE areaSegm(self)
IMPLICIT NONE
CLASS(meshCell1DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dJ
REAL(8):: dx(1)
CLASS(meshCell1DCartSegm), INTENT(inout):: self
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:2)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
self%volume = 0.D0
self%arNodes = 0.D0
!1D 1 point Gauss Quad Integral
Xi = 0.D0
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 2)
!Computes total volume of the cell
self%volume = detJ*2.D0
!Computes volume per node
self%arNodes = fPsi*self%volume
END IF
END SUBROUTINE areaSegm
CALL self%partialDer(2, dPsi, dx)
dJ = dx(1)
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes element Jacobian determinant
PURE FUNCTION detJ1DCart(pDer) RESULT(dJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: dJ
dJ = pDer(1, 1)
END FUNCTION detJ1DCart
!Computes the invers Jacobian
PURE FUNCTION invJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ)
!Computes element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ1DCart(pDer) RESULT(invJ)
IMPLICIT NONE
CLASS(meshCell1DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3)
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(Xi, 2)
END IF
invJ = 0.D0
CALL self%partialDer(2, dPsi, dx)
invJ(1,1) = 1.D0/dx(1)
invJ(1, 1) = 1.D0/pDer(1, 1)
invJ(2, 2) = 1.D0
invJ(3, 3) = 1.D0
END FUNCTION invJ1DCart

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

@ -14,6 +14,7 @@ MODULE moduleMesh1DRad
!Element coordinates
REAL(8):: r = 0.D0
CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode1DRad
PROCEDURE, PASS:: getCoordinates => getCoord1DRad
@ -25,6 +26,7 @@ MODULE moduleMesh1DRad
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge1DRad
PROCEDURE, PASS:: getNodes => getNodes1DRad
PROCEDURE, PASS:: intersection => intersection1DRad
@ -32,28 +34,7 @@ MODULE moduleMesh1DRad
END TYPE meshEdge1DRad
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DRad
CONTAINS
PROCEDURE, PASS:: detJac => detJ1DRad
PROCEDURE, PASS:: invJac => invJ1DRad
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshCell1DRad
ABSTRACT INTERFACE
PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx)
IMPORT meshCell1DRad
CLASS(meshCell1DRad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1):: dx
END SUBROUTINE partialDer_interface
END INTERFACE
TYPE, PUBLIC, EXTENDS(meshCell1DRad):: meshCell1DRadSegm
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DRadSegm
!Element coordinates
REAL(8):: r(1:2)
!Connectivity to nodes
@ -62,20 +43,24 @@ MODULE moduleMesh1DRad
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCell1DRadSegm
PROCEDURE, PASS:: randPos => randPos1DRadSeg
PROCEDURE, PASS:: area => areaRad
PROCEDURE, PASS:: fPsi => fPsiRad
PROCEDURE, PASS:: dPsi => dPsiRad
PROCEDURE, PASS:: partialDer => partialDerRad
PROCEDURE, PASS:: elemK => elemKRad
PROCEDURE, PASS:: elemF => elemFRad
PROCEDURE, PASS:: gatherElectricField => gatherEFRad
PROCEDURE, PASS:: gatherMagneticField => gatherMFRad
PROCEDURE, NOPASS:: inside => insideRad
PROCEDURE, PASS:: getNodes => getNodesRad
PROCEDURE, PASS:: phy2log => phy2logRad
PROCEDURE, PASS:: nextElement => nextElementRad
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: randPos => randPos1DRadSegm
PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: partialDer => partialDerSegm
PROCEDURE, NOPASS:: detJac => detJ1DRad
PROCEDURE, NOPASS:: invJac => invJ1DRad
PROCEDURE, PASS:: gatherElectricField => gatherEFSegm
PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm
PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: neighbourElement => neighbourElementSegm
!PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: area => areaSegm
END TYPE meshCell1DRadSegm
@ -139,7 +124,6 @@ MODULE moduleMesh1DRad
self%r = r1(1)
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
self%normal = self%normal/NORM2(self%normal)
!Boundary index
self%boundary => boundary(bt)
@ -221,8 +205,20 @@ MODULE moduleMesh1DRad
END SUBROUTINE initCell1DRadSegm
!Calculates a random position in 1D volume
FUNCTION randPos1DRadSeg(self) RESULT(r)
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesSegm
!Random position in 1D volume
FUNCTION randPos1DRadSegm(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
@ -231,152 +227,63 @@ MODULE moduleMesh1DRad
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:2)
Xi = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2:3) = 0.D0
fPsi = self%fPsi(Xi, 2)
r = 0.D0
r(1) = DOT_PRODUCT(fPsi, self%r)
END FUNCTION randPos1DRadSeg
!Computes element area
PURE SUBROUTINE areaRad(self)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2), fPsi_node(1:2)
REAL(8):: r
REAL(8):: detJ
REAL(8):: Xi(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!1 point Gauss integral
Xi = 0.D0
fPsi = self%fPsi(Xi, 2)
detJ = self%detJac(Xi, 2)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi, self%r)
l = 2.D0*detJ
self%volume = r*l
!Computes volume per node
Xi = (/-5.D-1, 0.D0, 0.D0/)
r = self%gatherF(Xi, 2, self%r)
self%arNodes(1) = fPsi(1)*r*l
Xi = (/ 5.D-1, 0.D0, 0.D0/)
r = self%gatherF(Xi, 2, self%r)
self%arNodes(2) = fPsi(2)*r*l
END SUBROUTINE areaRad
END FUNCTION randPos1DRadSegm
!Computes element functions at point Xi
PURE FUNCTION fPsiRad(self, Xi, nNodes) RESULT(fPsi)
PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = 1.D0 - Xi(1)
fPsi(2) = 1.D0 + Xi(1)
fPsi = (/ 1.D0 - Xi(1), &
1.D0 + Xi(1) /)
fPsi = fPsi * 5.D-1
fPsi = fPsi * 0.50D0
END FUNCTION fPsiRad
END FUNCTION fPsiSegm
!Computes element derivative shape function at Xi
PURE FUNCTION dPsiRad(self, Xi, nNodes) RESULT(dPsi)
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
dPsi(1, 1) = -5.D-1
dPsi(1, 2) = 5.D-1
dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
END FUNCTION dPsiRad
END FUNCTION dPsiSegm
!Computes partial derivatives of coordinates
PURE SUBROUTINE partialDerRad(self, nNodes, dPsi, dx)
!Partial derivative in global coordinates
PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1):: dx
REAL(8):: pDer(1:3, 1:3)
dx(1) = DOT_PRODUCT(dPsi(1,:), self%r)
pDer = 0.D0
END SUBROUTINE partialDerRad
pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%r(1:2))
pDer(2,2) = 1.D0
pDer(3,3) = 1.D0
!Computes local stiffness matrix
PURE FUNCTION elemKRad(self, nNodes) RESULT(localK)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
END FUNCTION partialDerSegm
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2)
REAL(8):: invJ(1:3,1:3), detJ
REAL(8):: r, fPsi(1:2)
INTEGER:: l
localK = 0.D0
Xi = 0.D0
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
detJ = self%detJac(Xi, 2, dPsi)
invJ = self%invJac(Xi, 2, dPsi)
fPsi = self%fPsi(Xi, 2)
r = DOT_PRODUCT(fPsi, self%r)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
r*wSeg(l)/detJ
END DO
localK = localK*PI2
END FUNCTION elemKRad
PURE FUNCTION elemFRad(self, nNodes, source) RESULT(localF)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:2)
REAL(8):: detJ, f, r
REAL(8):: Xi(1:3)
INTEGER:: l
localF = 0.D0
Xi = 0.D0
DO l = 1, 3
Xi(1) = corSeg(l)
detJ = self%detJac(Xi, 2)
fPsi = self%fPsi(Xi, 2)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*r*wSeg(l)*detJ
END DO
END FUNCTION elemFRad
PURE FUNCTION gatherEFRad(self, Xi) RESULT(array)
PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
@ -388,9 +295,9 @@ MODULE moduleMesh1DRad
array = -self%gatherDF(Xi, 2, phi)
END FUNCTION gatherEFRad
END FUNCTION gatherEFSegm
PURE FUNCTION gatherMFRad(self, Xi) RESULT(array)
PURE FUNCTION gatherMFSegm(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
@ -408,9 +315,79 @@ MODULE moduleMesh1DRad
array = self%gatherF(Xi, 2, B)
END FUNCTION gatherMFRad
END FUNCTION gatherMFSegm
PURE FUNCTION insideRad(Xi) RESULT(ins)
!Computes element local stiffness matrix
PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: r
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
localK = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
r = self%gatherF(Xi, 4, self%r)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
r*wSeg(l)/detJ
END DO
localK = localK*PI2
END FUNCTION elemKSegm
!Computes the local source vector for a force f
PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:2)
REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
REAL(8):: r
INTEGER:: l
localF = 0.D0
Xi = 0.D0
!Start 1D Gauss Quad Integral
DO l = 1, 3
Xi(1) = corSeg(l)
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 2)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi, source)
localF = localF + r*f*fPsi*wSeg(l)*detJ
END DO
localF = localF*PI2
END FUNCTION elemFSegm
PURE FUNCTION insideSegm(Xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi(1:3)
@ -419,102 +396,97 @@ MODULE moduleMesh1DRad
ins = Xi(1) >=-1.D0 .AND. &
Xi(1) <= 1.D0
END FUNCTION insideRad
END FUNCTION insideSegm
!Get nodes from 1D volume
PURE FUNCTION getNodesRad(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesRad
PURE FUNCTION phy2logRad(self, r) RESULT(xN)
PURE FUNCTION phy2logSegm(self, r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
REAL(8):: Xi(1:3)
xN = 0.D0
xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0
Xi = 0.D0
END FUNCTION phy2logRad
Xi(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0
!Get next element for a logical position Xi
SUBROUTINE nextElementRad(self, Xi, nextElement)
END FUNCTION phy2logSegm
!Get the next element for a logical position Xi
SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
NULLIFY(nextElement)
NULLIFY(neighbourElement)
IF (Xi(1) < -1.D0) THEN
nextElement => self%e2
neighbourElement => self%e2
ELSEIF (Xi(1) > 1.D0) THEN
nextElement => self%e1
neighbourElement => self%e1
END IF
END SUBROUTINE nextElementRad
END SUBROUTINE neighbourElementSegm
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes the element Jacobian determinant
PURE FUNCTION detJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(dJ)
!Computes element area
PURE SUBROUTINE areaSegm(self)
USE moduleConstParam, ONLY: PI
IMPLICIT NONE
CLASS(meshCell1DRad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dJ
REAL(8):: dx(1)
CLASS(meshCell1DRadSegm), INTENT(inout):: self
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:2)
REAL(8):: r
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
self%volume = 0.D0
self%arNodes = 0.D0
!1D 1 point Gauss Quad Integral
Xi = 0.D0
dPsi = self%dPsi(Xi, 2)
pDer = self%partialDer(2, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 2)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi, self%r)
self%volume = r*detJ*2.D0*PI !2PI
!Computes volume per node
Xi = (/-5.D-1, 0.D0, 0.D0/)
r = self%gatherF(Xi, 2, self%r)
self%arNodes(1) = fPsi(1)*self%volume
Xi = (/ 5.D-1, 0.D0, 0.D0/)
r = self%gatherF(Xi, 2, self%r)
self%arNodes(2) = fPsi(2)*self%volume
END IF
END SUBROUTINE areaSegm
CALL self%partialDer(nNodes, dPsi, dx)
dJ = dx(1)
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes element Jacobian determinant
PURE FUNCTION detJ1DRad(pDer) RESULT(dJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: dJ
dJ = pDer(1, 1)
END FUNCTION detJ1DRad
!Computes the invers Jacobian
PURE FUNCTION invJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(invJ)
!Computes element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ1DRad(pDer) RESULT(invJ)
IMPLICIT NONE
CLASS(meshCell1DRad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dx(1)
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(Xi, 2)
END IF
invJ = 0.D0
CALL self%partialDer(nNodes, dPsi, dx)
invJ(1,1) = 1.D0/dx(1)
invJ(1, 1) = 1.D0/pDer(1, 1)
invJ(2, 2) = 1.D0
invJ(3, 3) = 1.D0
END FUNCTION invJ1DRad

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

@ -19,6 +19,7 @@ MODULE moduleMesh2DCart
!Element coordinates
REAL(8):: x = 0.D0, y = 0.D0
CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode2DCart
PROCEDURE, PASS:: getCoordinates => getCoord2DCart
@ -30,6 +31,7 @@ MODULE moduleMesh2DCart
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
CONTAINS
!meshEdge DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initEdge2DCart
PROCEDURE, PASS:: getNodes => getNodes2DCart
PROCEDURE, PASS:: intersection => intersection2DCartEdge
@ -37,28 +39,8 @@ MODULE moduleMesh2DCart
END TYPE meshEdge2DCart
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCart
CONTAINS
PROCEDURE, PASS:: detJac => detJ2DCart
PROCEDURE, PASS:: invJac => invJ2DCart
PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer
END TYPE meshCell2DCart
ABSTRACT INTERFACE
PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy)
IMPORT meshCell2DCart
CLASS(meshCell2DCart), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy
END SUBROUTINE partialDer_interface
END INTERFACE
!Quadrilateral volume element
TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartQuad
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartQuad
!Element coordinates
REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0
!Connectivity to nodes
@ -68,25 +50,29 @@ MODULE moduleMesh2DCart
REAL(8):: arNodes(1:4) = 0.D0
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCellQuad2DCart
PROCEDURE, PASS:: getNodes => getNodesQuad
PROCEDURE, PASS:: randPos => randPosCellQuad
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, PASS:: fPsi => fPsiQuad
PROCEDURE, PASS:: dPsi => dPsiQuad
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad
PROCEDURE, PASS:: elemK => elemKQuad
PROCEDURE, PASS:: elemF => elemFQuad
PROCEDURE, NOPASS:: fPsi => fPsiQuad
PROCEDURE, NOPASS:: dPsi => dPsiQuad
PROCEDURE, PASS:: partialDer => partialDerQuad
PROCEDURE, NOPASS:: detJac => detJ2DCart
PROCEDURE, NOPASS:: invJac => invJ2DCart
PROCEDURE, PASS:: gatherElectricField => gatherEFQuad
PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad
PROCEDURE, PASS:: elemK => elemKQuad
PROCEDURE, PASS:: elemF => elemFQuad
PROCEDURE, NOPASS:: inside => insideQuad
PROCEDURE, PASS:: getNodes => getNodesQuad
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: nextElement => nextElementQuad
PROCEDURE, PASS:: neighbourElement => neighbourElementQuad
!PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: area => areaQuad
END TYPE meshCell2DCartQuad
!Triangular volume element
TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartTria
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartTria
!Element coordinates
REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0
!Connectivity to nodes
@ -96,20 +82,24 @@ MODULE moduleMesh2DCart
REAL(8):: arNodes(1:3) = 0.D0
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCellTria2DCart
PROCEDURE, PASS:: getNodes => getNodesTria
PROCEDURE, PASS:: randPos => randPosCellTria
PROCEDURE, PASS:: area => areaTria
PROCEDURE, PASS:: fPsi => fPsiTria
PROCEDURE, PASS:: dPsi => dPsiTria
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerTria
PROCEDURE, PASS:: elemK => elemKTria
PROCEDURE, PASS:: elemF => elemFTria
PROCEDURE, NOPASS:: fPsi => fPsiTria
PROCEDURE, NOPASS:: dPsi => dPsiTria
PROCEDURE, PASS:: partialDer => partialDerTria
PROCEDURE, NOPASS:: detJac => detJ2DCart
PROCEDURE, NOPASS:: invJac => invJ2DCart
PROCEDURE, PASS:: gatherElectricField => gatherEFTria
PROCEDURE, PASS:: gatherMagneticField => gatherMFTria
PROCEDURE, PASS:: elemK => elemKTria
PROCEDURE, PASS:: elemF => elemFTria
PROCEDURE, NOPASS:: inside => insideTria
PROCEDURE, PASS:: getNodes => getNodesTria
PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: nextElement => nextElementTria
PROCEDURE, PASS:: neighbourElement => neighbourElementTria
!PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: area => areaTria
END TYPE meshCell2DCartTria
@ -175,6 +165,7 @@ MODULE moduleMesh2DCart
r2 = self%n2%getCoordinates()
self%x = (/r1(1), r2(1)/)
self%y = (/r1(2), r2(2)/)
self%weight = 1.D0
!Normal vector
self%normal = (/ -(self%y(2)-self%y(1)), &
self%x(2)-self%x(1) , &
@ -290,84 +281,17 @@ MODULE moduleMesh2DCart
END SUBROUTINE initCellQuad2DCart
!Computes element area
PURE SUBROUTINE areaQuad(self)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout):: self
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)*4.D0 !4
fPsi = self%fPsi(Xi, 4)
self%volume = detJ
self%arNodes = fPsi*detJ
END SUBROUTINE areaQuad
!Computes element functions in point Xi
PURE FUNCTION fPsiQuad(self, Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2))
fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2))
fPsi(3) = (1.D0+Xi(1)) * (1.D0+Xi(2))
fPsi(4) = (1.D0-Xi(1)) * (1.D0+Xi(2))
fPsi = fPsi*0.25D0
END FUNCTION fPsiQuad
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiQuad(self, Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
dPsi(1,:) = (/ -(1.D0 - Xi(2)), &
(1.D0 - Xi(2)), &
(1.D0 + Xi(2)), &
-(1.D0 + Xi(2)) /)
dPsi(2,:) = (/ -(1.D0 - Xi(1)), &
-(1.D0 + Xi(1)), &
(1.D0 + Xi(1)), &
(1.D0 - Xi(1)) /)
dPsi = dPsi * 0.25D0
END FUNCTION dPsiQuad
!Partial derivative in global coordinates
PURE SUBROUTINE partialDerQuad(self, nNodes, dPsi, dx, dy)
!Get nodes from quadrilateral element
PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy
INTEGER:: n(1:nNodes)
dx = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /)
dy = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /)
n = (/ self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
END SUBROUTINE partialDerQuad
END FUNCTION getNodesQuad
!Random position in quadrilateral volume
FUNCTION randPosCellQuad(self) RESULT(r)
@ -379,78 +303,77 @@ MODULE moduleMesh2DCart
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
Xi = 0.D0
Xi(1) = random(-1.D0, 1.D0)
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(Xi, 4)
r = 0.D0
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randPosCellQuad
!Computes element local stiffness matrix
PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK)
!Computes element functions in point Xi
PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi = (/ (1.D0-Xi(1)) * (1.D0-Xi(2)), &
(1.D0+Xi(1)) * (1.D0-Xi(2)), &
(1.D0+Xi(1)) * (1.D0+Xi(2)), &
(1.D0-Xi(1)) * (1.D0+Xi(2)) /)
fPsi = fPsi * 0.25D0
END FUNCTION fPsiQuad
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), &
(1.D0 - Xi(2)), &
(1.D0 + Xi(2)), &
-(1.D0 + Xi(2)) /)
dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), &
-(1.D0 + Xi(1)), &
(1.D0 + Xi(1)), &
(1.D0 - Xi(1)) /)
dPsi = dPsi * 0.25D0
END FUNCTION dPsiQuad
!Partial derivative in global coordinates
PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3,1:4)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l, m
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8):: pDer(1:3, 1:3)
localK=0.D0
Xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 3
Xi(2) = corQuad(l)
DO m = 1, 3
Xi(1) = corQuad(m)
fPsi = self%fPsi(Xi, 4)
dPsi = self%dPsi(Xi, 4)
detJ = self%detJac(Xi, 4, dPsi)
invJ = self%invJac(Xi, 4, dPsi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
wQuad(l)*wQuad(m)/detJ
pDer = 0.D0
END DO
END DO
pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /)
pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /)
pDer(3,3) = 1.D0
END FUNCTION elemKQuad
!Computes the local source vector for a force f
PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: detJ, f
INTEGER:: l, m
localF = 0.D0
Xi = 0.D0
DO l=1, 3
Xi(1) = corQuad(l)
DO m = 1, 3
Xi(2) = corQuad(m)
detJ = self%detJac(Xi, 4)
fPsi = self%fPsi(Xi, 4)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ
END DO
END DO
END FUNCTION elemFQuad
END FUNCTION partialDerQuad
PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array)
IMPLICIT NONE
@ -494,6 +417,75 @@ MODULE moduleMesh2DCart
END FUNCTION gatherMFQuad
!Computes element local stiffness matrix
PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:4)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: r
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l, m
localK = 0.D0
Xi = 0.D0
!Start 2D Gauss Quad Integral
DO l = 1, 3
Xi(2) = corQuad(l)
DO m = 1, 3
Xi(1) = corQuad(m)
dPsi = self%dPsi(Xi, 4)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
wQuad(l)*wQuad(m)/detJ
END DO
END DO
END FUNCTION elemKQuad
!Computes the local source vector for a force f
PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: dPsi(1:3, 1:4)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: detJ, f
INTEGER:: l, m
localF = 0.D0
Xi = 0.D0
DO l = 1, 3
Xi(1) = corQuad(l)
DO m = 1, 3
Xi(2) = corQuad(m)
dPsi = self%dPsi(Xi, 4)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ
END DO
END DO
END FUNCTION elemFQuad
!Checks if a particle is inside a quad element
PURE FUNCTION insideQuad(Xi) RESULT(ins)
IMPLICIT NONE
@ -506,18 +498,6 @@ MODULE moduleMesh2DCart
END FUNCTION insideQuad
!Gets nodes from quadrilateral element
PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
END FUNCTION getNodesQuad
!Transforms physical coordinates to element coordinates
PURE FUNCTION phy2logQuad(self,r) RESULT(Xi)
IMPLICIT NONE
@ -527,6 +507,7 @@ MODULE moduleMesh2DCart
REAL(8):: Xi(1:3)
REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: conv
!Iterative newton method to transform coordinates
@ -535,7 +516,9 @@ MODULE moduleMesh2DCart
DO WHILE(conv > 1.D-2)
dPsi = self%dPsi(XiO, 4)
invJ = self%invJac(XiO, 4, dPsi)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
fPsi = self%fPsi(XiO, 4)
f = (/ DOT_PRODUCT(fPsi,self%x), &
DOT_PRODUCT(fPsi,self%y), &
@ -550,31 +533,56 @@ MODULE moduleMesh2DCart
END FUNCTION phy2logQuad
!Gets the next element for a logical position Xi
SUBROUTINE nextElementQuad(self, Xi, nextElement)
SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
REAL(8):: XiArray(1:4)
INTEGER:: nextInt
XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
nextInt = MAXLOC(XiArray,1)
!Selects the higher value of directions and searches in that direction
NULLIFY(nextElement)
NULLIFY(neighbourElement)
SELECT CASE (nextInt)
CASE (1)
nextElement => self%e1
neighbourElement => self%e1
CASE (2)
nextElement => self%e2
neighbourElement => self%e2
CASE (3)
nextElement => self%e3
neighbourElement => self%e3
CASE (4)
nextElement => self%e4
neighbourElement => self%e4
END SELECT
END SUBROUTINE nextElementQuad
END SUBROUTINE neighbourElementQuad
!Computes element area
PURE SUBROUTINE areaQuad(self)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout):: self
REAL(8):: Xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:4)
REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
Xi = 0.D0
dPsi = self%dPsi(Xi, 4)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4)
!Computes total volume of the cell
self%volume = detJ
!Computes volume per node
self%arNodes = fPsi*detJ
END SUBROUTINE areaQuad
!TRIA ELEMENT
!Init tria element
@ -617,6 +625,18 @@ MODULE moduleMesh2DCart
END SUBROUTINE initCellTria2DCart
!Gets node indexes from triangular element
PURE FUNCTION getNodesTria(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/self%n1%n, self%n2%n, self%n3%n /)
END FUNCTION getNodesTria
!Random position in quadrilateral volume
FUNCTION randPosCellTria(self) RESULT(r)
USE moduleRandom
@ -639,31 +659,10 @@ MODULE moduleMesh2DCart
END FUNCTION randPosCellTria
!Calculates area for triangular element
PURE SUBROUTINE areaTria(self)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(inout):: self
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, 4)/2.D0
fPsi = self%fPsi(Xi, 4)
self%volume = detJ
self%arNodes = fPsi*detJ
END SUBROUTINE areaTria
!Shape functions for triangular element
PURE FUNCTION fPsiTria(self, Xi, nNodes) RESULT(fPsi)
PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
@ -675,10 +674,9 @@ MODULE moduleMesh2DCart
END FUNCTION fPsiTria
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiTria(self, Xi, nNodes) RESULT(dPsi)
PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
@ -690,76 +688,22 @@ MODULE moduleMesh2DCart
END FUNCTION dPsiTria
PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dx, dy)
PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy
REAL(8):: pDer(1:3, 1:3)
dx = (/ DOT_PRODUCT(dPsi(1,:),self%x), &
DOT_PRODUCT(dPsi(2,:),self%x) /)
dy = (/ DOT_PRODUCT(dPsi(1,:),self%y), &
DOT_PRODUCT(dPsi(2,:),self%y) /)
pDer = 0.D0
END SUBROUTINE partialDerTria
pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%x(1:3)), &
DOT_PRODUCT(dPsi(2,1:3),self%x(1:3)) /)
pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), &
DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /)
!Computes element local stiffness matrix
PURE FUNCTION elemKTria(self, nNodes) RESULT(localK)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3), dPsi(1:3,1:3)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
localK=0.D0
Xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 4
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
dPsi = self%dPsi(Xi, 3)
detJ = self%detJac(Xi, 3, dPsi)
invJ = self%invJac(Xi, 3, dPsi)
fPsi = self%fPsi(Xi, 3)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
END DO
END FUNCTION elemKTria
!Computes element local source vector
PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
INTEGER:: l
localF = 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, 3)
fPsi = self%fPsi(Xi, 3)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wTria(l)*detJ
END DO
END FUNCTION elemFTria
END FUNCTION partialDerTria
PURE FUNCTION gatherEFTria(self, Xi) RESULT(array)
IMPLICIT NONE
@ -799,6 +743,66 @@ MODULE moduleMesh2DCart
END FUNCTION gatherMFTria
!Computes element local stiffness matrix
PURE FUNCTION elemKTria(self, nNodes) RESULT(localK)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3,1:3)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
localK=0.D0
Xi=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 4
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
dPsi = self%dPsi(Xi, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
END DO
END FUNCTION elemKTria
!Computes element local source vector
PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: fPsi(1:3)
REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, f
INTEGER:: l
localF = 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, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 3)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wTria(l)*detJ
END DO
END FUNCTION elemFTria
PURE FUNCTION insideTria(Xi) RESULT(ins)
IMPLICIT NONE
@ -811,18 +815,6 @@ MODULE moduleMesh2DCart
END FUNCTION insideTria
!Gets node indexes from triangular element
PURE FUNCTION getNodesTria(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/self%n1%n, self%n2%n, self%n3%n /)
END FUNCTION getNodesTria
!Transforms physical coordinates to element coordinates
PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
IMPLICIT NONE
@ -830,96 +822,94 @@ MODULE moduleMesh2DCart
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: r(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:3)
REAL(8):: dPsi(1:3, 1:3)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3, 1:3), detJ
!Direct method to convert coordinates
Xi = 0.D0
deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /)
dPsi = self%dPsi(Xi, 3)
invJ = self%invJac(Xi, 3, dPsi)
detJ = self%detJac(Xi, 3, dPsi)
pDer = self%partialDer(3, dPsi)
invJ = self%invJac(pDer)
detJ = self%detJac(pDer)
Xi = MATMUL(invJ,deltaR)/detJ
END FUNCTION phy2logTria
SUBROUTINE nextElementTria(self, Xi, nextElement)
SUBROUTINE neighbourElementTria(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
REAL(8):: XiArray(1:3)
INTEGER:: nextInt
XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /)
nextInt = MINLOC(XiArray,1)
NULLIFY(nextElement)
NULLIFY(neighbourElement)
SELECT CASE (nextInt)
CASE (1)
nextElement => self%e1
neighbourElement => self%e1
CASE (2)
nextElement => self%e2
neighbourElement => self%e2
CASE (3)
nextElement => self%e3
neighbourElement => self%e3
END SELECT
END SUBROUTINE nextElementTria
END SUBROUTINE neighbourElementTria
!Calculates area for triangular element
PURE SUBROUTINE areaTria(self)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(inout):: self
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
dPsi = self%dPsi(Xi, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4)
!Computes total volume of the cell
self%volume = detJ
!Computes volume per node
self%arNodes = fPsi*detJ
END SUBROUTINE areaTria
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D
!Computes element Jacobian determinant
PURE FUNCTION detJ2DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ)
PURE FUNCTION detJ2DCart(pDer) RESULT(dJ)
IMPLICIT NONE
CLASS(meshCell2DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: dJ
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dx(1:2), dy(1:2)
IF(PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(Xi, 4)
END IF
CALL self%partialDer(nNodes, dPsi, dx, dy)
dJ = dx(1)*dy(2)-dx(2)*dy(1)
dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1)
END FUNCTION detJ2DCart
!Computes element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ2DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ)
PURE FUNCTION invJ2DCart(pDer) RESULT(invJ)
IMPLICIT NONE
CLASS(meshCell2DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3)
REAL(8):: dPsi(1:3,1:nNodes)
REAL(8):: dx(1:2), dy(1:2)
IF(PRESENT(dPsi_in)) THEN
dPsi=dPsi_in
ELSE
dPsi = self%dPsi(Xi, 4)
END IF
invJ = 0.D0
CALL self%partialDer(nNodes, dPsi, dx, dy)
invJ(1,1:2) = (/ dy(2), -dx(2) /)
invJ(2,1:2) = (/ -dy(1), dx(1) /)
invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /)
invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /)
invJ(3, 3) = 1.D0
END FUNCTION invJ2DCart

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

@ -67,7 +67,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: neighbourElement => neighbourElementQuad
!PARTICLUAR PROCEDURES
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, PASS, PRIVATE:: area => areaQuad
END TYPE meshCell2DCylQuad
@ -99,7 +99,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: neighbourElement => neighbourElementTria
!PARTICULAR PROCEDURES
PROCEDURE, PASS:: area => areaTria
PROCEDURE, PASS, PRIVATE:: area => areaTria
END TYPE meshCell2DCylTria
@ -256,8 +256,13 @@ MODULE moduleMesh2DCyl
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3, r4
!Assign node index
self%n = n
!Assign number of nodes of cell
self%nNodes = SIZE(p)
!Assign nodes to element
self%n1 => nodes(p(1))%obj
self%n2 => nodes(p(2))%obj
self%n3 => nodes(p(3))%obj
@ -428,7 +433,7 @@ MODULE moduleMesh2DCyl
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: dPsi(1:3, 1:4)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: r
REAL(8):: invJ(1:3,1:3), detJ
@ -445,8 +450,7 @@ MODULE moduleMesh2DCyl
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
fPsi = self%fPsi(Xi, 4)
r = DOT_PRODUCT(fPsi,self%r)
r = self%gatherF(Xi, 4, self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
r*wQuad(l)*wQuad(m)/detJ
@ -467,7 +471,8 @@ MODULE moduleMesh2DCyl
REAL(8), INTENT(in):: source(1:nNodes)
REAL(8):: localF(1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: fPsi(1:4)
REAL(8):: dPsi(1:3, 1:4)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: r
REAL(8):: detJ, f
@ -483,7 +488,7 @@ MODULE moduleMesh2DCyl
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4)
r = DOT_PRODUCT(fPsi,self%r)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi,source)
localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ
@ -539,7 +544,7 @@ MODULE moduleMesh2DCyl
END FUNCTION phy2logQuad
!Gets the next element for a logical position Xi
!Get the next element for a logical position Xi
SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement)
IMPLICIT NONE
@ -572,7 +577,8 @@ MODULE moduleMesh2DCyl
IMPLICIT NONE
CLASS(meshCell2DCylQuad), INTENT(inout):: self
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: detJ
REAL(8):: fPsi(1:4)
REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
@ -583,24 +589,24 @@ MODULE moduleMesh2DCyl
Xi = 0.D0
dPsi = self%dPsi(Xi, 4)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)*PI8 !4*2*pi
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
self%volume = r*detJ*PI8 !4*2*pi
!Computes volume per node
Xi = (/-5.D-1, -5.D-1, 0.D0/)
r = self%gatherF(Xi, 4, self%r)
self%arNodes(1) = fPsi(1)*r*detJ
self%arNodes(1) = fPsi(1)*self%volume
Xi = (/ 5.D-1, -5.D-1, 0.D0/)
r = self%gatherF(Xi, 4, self%r)
self%arNodes(2) = fPsi(2)*r*detJ
self%arNodes(2) = fPsi(2)*self%volume
Xi = (/ 5.D-1, 5.D-1, 0.D0/)
r = self%gatherF(Xi, 4, self%r)
self%arNodes(3) = fPsi(3)*r*detJ
self%arNodes(3) = fPsi(3)*self%volume
Xi = (/-5.D-1, 5.D-1, 0.D0/)
r = self%gatherF(Xi, 4, self%r)
self%arNodes(4) = fPsi(4)*r*detJ
self%arNodes(4) = fPsi(4)*self%volume
END SUBROUTINE areaQuad
@ -619,7 +625,7 @@ MODULE moduleMesh2DCyl
!Assign node index
self%n = n
!Assign nomber of nodes to cell
!Assign number of nodes of cell
self%nNodes = SIZE(p)
!Assign nodes to element
@ -736,7 +742,7 @@ MODULE moduleMesh2DCyl
self%n2%emData%phi, &
self%n3%emData%phi /)
array = -self%gatherDF(Xi, 4, phi)
array = -self%gatherDF(Xi, 3, phi)
END FUNCTION gatherEFTria
@ -772,8 +778,8 @@ MODULE moduleMesh2DCyl
INTEGER, INTENT(in):: nNodes
REAL(8):: localK(1:nNodes,1:nNodes)
REAL(8):: Xi(1:3)
REAL(8):: dPsi(1:3,1:3)
REAL(8):: r
REAL(8):: fPsi(1:3), dPsi(1:3,1:3)
REAL(8):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
@ -788,8 +794,7 @@ MODULE moduleMesh2DCyl
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
invJ = self%invJac(pDer)
fPsi = self%fPsi(Xi, 3)
r = DOT_PRODUCT(fPsi,self%r)
r = self%gatherF(Xi, 3, self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
END DO
@ -809,22 +814,23 @@ MODULE moduleMesh2DCyl
REAL(8):: fPsi(1:3)
REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: detJ, f
REAL(8):: r
INTEGER:: l
localF = 0.D0
Xi = 0.D0
!Start 2D Gauss Quad Integral
DO l=1, 4
DO l = 1, 4
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
dPsi = self%dPsi(Xi, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 3)
r = DOT_PRODUCT(fPsi,self%r)
f = DOT_PRODUCT(fPsi,source)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi, source)
localF = localF + r*f*fPsi*wTria(l)*detJ
END DO
@ -897,10 +903,10 @@ MODULE moduleMesh2DCyl
CLASS(meshCell2DCylTria), INTENT(inout):: self
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:3)
REAL(8):: r
self%volume = 0.D0
self%arNodes = 0.D0
@ -908,13 +914,13 @@ MODULE moduleMesh2DCyl
Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
dPsi = self%dPsi(Xi, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)*PI !2PI*1/2
fPsi = self%fPsi(Xi, 4)
detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 3)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
r = DOT_PRODUCT(fPsi, self%r)
self%volume = r*detJ*PI !2PI*1/2
!Computes volume per node
self%arNodes = fPsi*r*detJ
self%arNodes = fPsi*self%volume
END SUBROUTINE areaTria
@ -939,9 +945,9 @@ MODULE moduleMesh2DCyl
invJ = 0.D0
invJ(1,1:2) = (/ pDer(2,2), -pDer(1,2) /)
invJ(2,1:2) = (/ -pDer(2,1), pDer(1,1) /)
invJ(3,3) = 1.D0
invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /)
invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /)
invJ(3, 3) = 1.D0
END FUNCTION invJ2DCyl

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

@ -30,7 +30,7 @@ MODULE moduleMesh3DCart
PROCEDURE, PASS:: intersection => intersection3DCartTria
PROCEDURE, PASS:: randPos => randPosEdgeTria
!PARTICULAR PROCEDURES
PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria
PROCEDURE, NOPASS, PRIVATE:: fPsi => fPsiEdgeTria
END TYPE meshEdge3DCartTria
@ -60,7 +60,7 @@ MODULE moduleMesh3DCart
PROCEDURE, PASS:: phy2log => phy2logTetra
PROCEDURE, PASS:: neighbourElement => neighbourElementTetra
!PARTICULAR PROCEDURES
PROCEDURE, PASS:: calcVol => volumeTetra
PROCEDURE, PASS, PRIVATE:: calcVol => volumeTetra
END TYPE meshCell3DCartTetra