!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 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):: 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 conv = 1.D0 XiO = 0.D0 f(3) = 0.D0 DO WHILE(conv > 1.D-4) dPsi = self%dPsi(XiO, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) fPsi = self%fPsi(XiO, 4) f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & DOT_PRODUCT(fPsi,self%y) /) - r(1:2) Xi = XiO - MATMUL(invJ, f)/detJ conv = MAXVAL(DABS(Xi-XiO),1) XiO = Xi 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.D0, 1.D0, 0.D0 /) dPsi(2,:) = (/ -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):: deltaR(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) pDer = self%partialDer(3, dPsi) invJ = self%invJac(pDer) detJ = self%detJac(pDer) Xi = MATMUL(invJ,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(1,2) /) invJ(2, 1:2) = (/ -pDer(2,1), 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