From f86641110ce6bf86e55ac7cd600dbcc3f2f25976 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 9 Feb 2023 12:15:18 +0100 Subject: [PATCH 1/2] Issue in connectivity between Quad elements For unstructured meshes, the quad elements where not being propertly connected as not all possibilities were tested. This should be fixed now. --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 108 +++++++++++++++---- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 108 +++++++++++++++---- 2 files changed, 176 insertions(+), 40 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index c02078a..d50263b 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -1012,38 +1012,106 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartQuad), INTENT(inout), TARGET:: elemB !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA + 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) .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n4%n) THEN - elemA%e2 => elemB - elemB%e4 => elemA + 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) .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA + 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) .AND. & - elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA + 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 diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 307f71c..f6f2b26 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -1050,38 +1050,106 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemB !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA + 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) .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n4%n) THEN - elemA%e2 => elemB - elemB%e4 => elemA + 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) .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA + 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) .AND. & - elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA + 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 From de1d4567f374e8fc94349b9d619c505adf8e480f Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 9 Feb 2023 15:32:04 +0100 Subject: [PATCH 2/2] Issue for particles in quad cell Due to a high convergence value (1.0e-2) in phy2logQuad (variable conv), particles were being stuck in some elements, reaching a segmentation fault. The new limit (1.0e-4) should avoid this. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 2 +- src/modules/mesh/moduleMesh.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index f6f2b26..f8e41b3 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -524,7 +524,7 @@ MODULE moduleMesh2DCyl conv = 1.D0 XiO = 0.D0 - DO WHILE(conv > 1.D-2) + DO WHILE(conv > 1.D-4) dPsi = self%dPsi(XiO, 4) pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index cb252e8..0804059 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -639,8 +639,8 @@ MODULE moduleMesh CLASS(meshCell), INTENT(inout):: self CLASS(particle), INTENT(inout), TARGET:: part CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell - REAL(8):: Xi(1:3) - CLASS(meshElement), POINTER:: neighbourElement + REAL(8):: Xi(1:3) = 0.D0 + CLASS(meshElement), POINTER:: neighbourElement => NULL() INTEGER:: sp Xi = self%phy2log(part%r)