JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
fpakc/src/modules/mesh/2DCart/moduleMesh2DCart.f90
JGonzalez 3d7b1ce476 Fixed! 
So it seems that rectangles and triangles are now working properly.

I have also checked the phy2log routine, now it seems a bit more complicated, but it is much clearer.
Maybe in the future is worth rethinking to improve speed (specially for quad elements)
2025-08-03 22:14:19 +02:00

1441 lines
40 KiB
Fortran
Raw Blame History

!moduleMesh2DCart: 2D Cartesian coordinate system
! x == x
! y == y
! z == unused
MODULE moduleMesh2DCart
USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE
!Values for Gauss integral
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:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCart
!Element coordinates
REAL(8):: x = 0.D0, y = 0.D0
CONTAINS
!meshNode DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initNode2DCart
PROCEDURE, PASS:: getCoordinates => getCoord2DCart
END TYPE meshNode2DCart
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCart
!Element coordinates
REAL(8):: x(1:2) = 0.D0, y(1:2) = 0.D0
!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
PROCEDURE, PASS:: randPos => randPosEdge
END TYPE meshEdge2DCart
!Quadrilateral volume element
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartQuad
!Element coordinates
REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
!Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCellQuad2DCart
PROCEDURE, PASS:: getNodes => getNodesQuad
PROCEDURE, PASS:: randPos => randPosCellQuad
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:: phy2log => phy2logQuad
PROCEDURE, PASS:: neighbourElement => neighbourElementQuad
!PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad
END TYPE meshCell2DCartQuad
!Triangular volume element
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartTria
!Element coordinates
REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
!Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
CONTAINS
!meshCell DEFERRED PROCEDURES
PROCEDURE, PASS:: init => initCellTria2DCart
PROCEDURE, PASS:: getNodes => getNodesTria
PROCEDURE, PASS:: randPos => randPosCellTria
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:: phy2log => phy2logTria
PROCEDURE, PASS:: neighbourElement => neighbourElementTria
!PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria
END TYPE meshCell2DCartTria
CONTAINS
!NODE FUNCTIONS
!Inits node element
SUBROUTINE initNode2DCart(self, n, r)
USE moduleSpecies
USE moduleRefParam
USE OMP_LIB
IMPLICIT NONE
CLASS(meshNode2DCart), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%x = r(1)/L_ref
self%y = r(2)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output:
ALLOCATE(self%output(1:nSpecies))
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initNode2DCart
!Get coordinates from node
PURE FUNCTION getCoord2DCart(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNode2DCart), INTENT(in):: self
REAL(8):: r(1:3)
r = (/self%x, self%y, 0.D0/)
END FUNCTION getCoord2DCart
!EDGE FUNCTIONS
!Init edge element
SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface)
USE moduleSpecies
USE moduleBoundary
USE moduleErrors
USE moduleRefParam, ONLY: L_ref
IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
REAL(8), DIMENSION(1:3):: r1, r2
INTEGER:: s
self%n = n
self%nNodes = SIZE(p)
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
self%x = (/r1(1), r2(1)/)
self%y = (/r1(2), r2(2)/)
self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref
!Normal vector
self%normal = (/ -(self%y(2)-self%y(1)), &
self%x(2)-self%x(1) , &
0.D0 /)
self%normal = self%normal/NORM2(self%normal)
!Boundary index
self%boundary => boundary(bt)
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
CALL pointBoundaryFunction(self, s)
END DO
!Physical surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdge2DCart
!Get nodes from edge
PURE FUNCTION getNodes2DCart(self, nNodes) RESULT(n)
IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(in):: self
INTEGER, INTENT(in):: nNodes
INTEGER:: n(1:nNodes)
n = (/self%n1%n, self%n2%n /)
END FUNCTION getNodes2DCart
!Calculate intersection between position and edge
PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0
REAL(8), DIMENSION(1:3):: r
REAL(8), DIMENSION(1:3):: edge0, edgeV
REAL(8):: tI
edge0 = (/self%x(1), self%y(1), 0.D0 /)
edgeV = (/self%x(2), self%y(2), 0.D0 /) - edge0
tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)
r = edge0 + tI*edgeV
END FUNCTION intersection2DCartEdge
!Calculate a random position in edge
FUNCTION randPosEdge(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(in):: self
REAL(8):: rnd
REAL(8):: r(1:3)
REAL(8):: p1(1:2), p2(1:2)
rnd = random()
p1 = (/self%x(1), self%y(1) /)
p2 = (/self%x(2), self%y(2) /)
r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
r(3) = 0.D0
END FUNCTION randPosEdge
!VOLUME FUNCTIONS
!QUAD FUNCTIONS
!Init element
SUBROUTINE initCellQuad2DCart(self, n, p, nodes)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
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
self%n4 => nodes(p(4))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
r3 = self%n3%getCoordinates()
r4 = self%n4%getCoordinates()
self%x = (/r1(1), r2(1), r3(1), r4(1)/)
self%y = (/r1(2), r2(2), r3(2), r4(2)/)
!Assign node volume
CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initCellQuad2DCart
!Get 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
!Random position in quadrilateral volume
FUNCTION randPosCellQuad(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
Xi(1) = random(-1.D0, 1.D0)
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(Xi, 4)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randPosCellQuad
!Compute 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 = 0.D0
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), INTENT(in):: dPsi(1:3,1:nNodes)
REAL(8):: pDer(1:3, 1:3)
pDer = 0.D0
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 partialDerQuad
PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: array(1:3)
REAL(8):: phi(1:4)
phi = (/ self%n1%emData%phi, &
self%n2%emData%phi, &
self%n3%emData%phi, &
self%n4%emData%phi /)
array = -self%gatherDF(Xi, 4, phi)
END FUNCTION gatherEFQuad
PURE FUNCTION gatherMFQuad(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: array(1:3)
REAL(8):: B(1:4,1:3)
B(:,1) = (/ self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1), &
self%n4%emData%B(1) /)
B(:,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2), &
self%n4%emData%B(2) /)
B(:,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3), &
self%n4%emData%B(3) /)
array = self%gatherF(Xi, 4, B)
END FUNCTION gatherMFQuad
!Compute 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):: 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
!Check if Xi is inside the element
PURE FUNCTION insideQuad(Xi) RESULT(ins)
IMPLICIT NONE
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)
END FUNCTION insideQuad
!Transform physical coordinates to element coordinates with a Taylor series
PURE FUNCTION phy2logQuad(self,r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
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
conv = 1.D0
Xi0 = 0.D0
Xi(3) = 0.D0
DO WHILE(conv > 1.D-4)
fPsi = self%fPsi(Xi0, 4)
x0(1) = dot_product(fPsi, self%x)
x0(2) = dot_product(fPsi, self%y)
deltaR = r(1:2) - x0
dPsi = self%dPsi(Xi0, 4)
pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer)
pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /)
pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /)
Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
Xi0(1:2) = Xi(1:2)
END DO
END FUNCTION phy2logQuad
!Get the neighbour element for a logical position Xi
SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
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)
!Select the higher value of directions and searches in that direction
NULLIFY(neighbourElement)
SELECT CASE (nextInt)
CASE (1)
neighbourElement => self%e1
CASE (2)
neighbourElement => self%e2
CASE (3)
neighbourElement => self%e3
CASE (4)
neighbourElement => self%e4
END SELECT
END SUBROUTINE neighbourElementQuad
!Compute element volume
PURE SUBROUTINE volumeQuad(self)
USE moduleRefParam, ONLY: L_ref
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
!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)
!Compute total volume of the cell
self%volume = detJ*4.D0/L_ref
!Compute volume per node
self%n1%v = self%n1%v + fPsi(1)*self%volume
self%n2%v = self%n2%v + fPsi(2)*self%volume
self%n3%v = self%n3%v + fPsi(3)*self%volume
self%n4%v = self%n4%v + fPsi(4)*self%volume
END SUBROUTINE volumeQuad
!TRIA FUNCTIONS
!Init element
SUBROUTINE initCellTria2DCart(self, n, p, nodes)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3
!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
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
r3 = self%n3%getCoordinates()
self%x = (/r1(1), r2(1), r3(1)/)
self%y = (/r1(2), r2(2), r3(2)/)
!Assign node volume
CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initCellTria2DCart
!Random position in cell
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 cell
FUNCTION randPosCellTria(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:3)
Xi(1) = random( 0.D0, 1.D0)
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
Xi(3) = 0.D0
fPsi = self%fPsi(Xi, 4)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randPosCellTria
!Compute element functions in point Xi
PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi(1:3)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = 1.D0 - Xi(1) - Xi(2)
fPsi(2) = Xi(1)
fPsi(3) = Xi(2)
END FUNCTION fPsiTria
!Compute element derivative functions in point Xi
PURE FUNCTION dPsiTria(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:3) = (/ -1.D0, 1.D0, 0.D0 /)
dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
END FUNCTION dPsiTria
!Compute the derivatives in global coordinates
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):: pDer(1:3, 1:3)
pDer = 0.D0
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)) /)
END FUNCTION partialDerTria
!Gather electric field at position Xi
PURE FUNCTION gatherEFTria(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: array(1:3)
REAL(8):: phi(1:3)
phi = (/ self%n1%emData%phi, &
self%n2%emData%phi, &
self%n3%emData%phi /)
array = -self%gatherDF(Xi, 3, phi)
END FUNCTION gatherEFTria
!Gather magnetic field at position Xi
PURE FUNCTION gatherMFTria(self, Xi) RESULT(array)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: array(1:3)
REAL(8):: B(1:3,1:3)
B(:,1) = (/ self%n1%emData%B(1), &
self%n2%emData%B(1), &
self%n3%emData%B(1) /)
B(:,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2), &
self%n3%emData%B(2) /)
B(:,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3), &
self%n3%emData%B(3) /)
array = self%gatherF(Xi, 3, B)
END FUNCTION gatherMFTria
!Compute cell 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
!Compute 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
!Check if Xi is inside the element
PURE FUNCTION insideTria(Xi) RESULT(ins)
IMPLICIT NONE
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
END FUNCTION insideTria
!Transform physical coordinates to element coordinates
PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
REAL(8):: dPsi(1:3,1:4)
REAL(8):: pDer(1:3, 1:3)
!Direct method to convert coordinates
Xi(3) = 0.D0
deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /)
dPsi = self%dPsi(Xi, 3)
pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer)
pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /)
pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /)
Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ
END FUNCTION phy2logTria
!Get the neighbour cell for a logical position Xi
SUBROUTINE neighbourElementTria(self, Xi, neighbourElement)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
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(neighbourElement)
SELECT CASE (nextInt)
CASE (1)
neighbourElement => self%e1
CASE (2)
neighbourElement => self%e2
CASE (3)
neighbourElement => self%e3
END SELECT
END SUBROUTINE neighbourElementTria
!Calculate volume for triangular element
PURE SUBROUTINE volumeTria(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
!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, 3)
!Computes total volume of the cell
self%volume = detJ
!Computes volume per node
self%n1%v = self%n1%v + fPsi(1)*self%volume
self%n2%v = self%n2%v + fPsi(2)*self%volume
self%n3%v = self%n3%v + fPsi(3)*self%volume
END SUBROUTINE volumeTria
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D
!Compute element Jacobian determinant
PURE FUNCTION detJ2DCart(pDer) RESULT(dJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: dJ
dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1)
END FUNCTION detJ2DCart
!Compute element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJ2DCart(pDer) RESULT(invJ)
IMPLICIT NONE
REAL(8), INTENT(in):: pDer(1:3, 1:3)
REAL(8):: invJ(1:3,1:3)
invJ = 0.D0
invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /)
invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /)
invJ(3, 3) = 1.D0
END FUNCTION invJ2DCart
SUBROUTINE connectMesh2DCart(self)
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self
INTEGER:: e, et
DO e = 1, self%numCells
!Connect Cell-Cell
DO et = 1, self%numCells
IF (e /= et) THEN
CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)
END IF
END DO
SELECT TYPE(self)
TYPE IS(meshParticles)
!Connect Cell-Edge
DO et = 1, self%numEdges
CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)
END DO
END SELECT
END DO
END SUBROUTINE connectMesh2DCart
!Selects type of elements to build connection
SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshCell), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshCell2DCartQuad)
!Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshCell2DCartQuad)
!Element B is a quadrilateral
CALL connectQuadQuad(elemA, elemB)
TYPE IS(meshCell2DCartTria)
!Element B is a triangle
CALL connectQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshCell2DCartTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshCell2DCartQuad)
!Element B is a quadrilateral
CALL connectQuadTria(elemB, elemA)
TYPE IS(meshCell2DCartTria)
!Element B is a triangle
CALL connectTriaTria(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectCellCell
SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemB)
CLASS IS(meshEdge2DCart)
SELECT TYPE(elemA)
TYPE IS(meshCell2DCartQuad)
!Element A is a quadrilateral
CALL connectQuadEdge(elemA, elemB)
TYPE IS(meshCell2DCartTria)
!Element A is a triangle
CALL connectTriaEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectCellEdge
SUBROUTINE connectQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n4%n) THEN
elemA%e1 => elemB
elemB%e4 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n4%n) THEN
elemA%e2 => elemB
elemB%e4 => elemA
ELSEIF (elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n4%n) THEN
elemA%e3 => elemB
elemB%e4 => elemA
ELSEIF (elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n4%n) THEN
elemA%e4 => elemB
elemB%e4 => elemA
ELSEIF (elemA%n4%n == elemB%n4%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e4 => elemB
elemB%e3 => elemA
END IF
END IF
END SUBROUTINE connectQuadQuad
SUBROUTINE connectQuadTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e4 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectQuadTria
SUBROUTINE connectTriaTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemA
CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
END SUBROUTINE connectTriaTria
SUBROUTINE connectQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
END SUBROUTINE connectQuadEdge
SUBROUTINE connectTriaEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell2DCartTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF
END SUBROUTINE connectTriaEdge
END MODULE moduleMesh2DCart
Reference in a new issue View git blame Copy permalink