diff --git a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 b/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 index 720e61d..d69cfdc 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 @@ -217,6 +217,9 @@ MODULE moduleMesh1DCartRead elemA%e1 => elemB elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 b/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 index 1a9f25c..4f14b5f 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 @@ -218,6 +218,9 @@ MODULE moduleMesh1DRadRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + END IF IF (.NOT. ASSOCIATED(elemA%e2) .AND. & diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 index b91627a..b0f2ee3 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -439,6 +439,9 @@ MODULE moduleMesh2DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -455,6 +458,9 @@ MODULE moduleMesh2DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -471,6 +477,9 @@ MODULE moduleMesh2DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB @@ -487,6 +496,9 @@ MODULE moduleMesh2DCartRead elemA%e4 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB @@ -511,6 +523,9 @@ MODULE moduleMesh2DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -527,6 +542,9 @@ MODULE moduleMesh2DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -543,6 +561,9 @@ MODULE moduleMesh2DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 index 085c3c3..379a0fb 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 @@ -439,6 +439,9 @@ MODULE moduleMesh2DCylRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -455,6 +458,9 @@ MODULE moduleMesh2DCylRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -471,6 +477,9 @@ MODULE moduleMesh2DCylRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB @@ -511,6 +520,9 @@ MODULE moduleMesh2DCylRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -527,6 +539,9 @@ MODULE moduleMesh2DCylRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -543,6 +558,9 @@ MODULE moduleMesh2DCylRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 1c713c9..0464c7f 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -146,7 +146,7 @@ MODULE moduleMesh3DCart self%n = n self%n1 => mesh%nodes(p(1))%obj - self%n3 => mesh%nodes(p(2))%obj + self%n2 => mesh%nodes(p(2))%obj self%n3 => mesh%nodes(p(3))%obj !Get element coordinates r1 = self%n1%getCoordinates() diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 index 40f81d0..c6b021a 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 @@ -207,6 +207,28 @@ MODULE moduleMesh3DCartRead END SUBROUTINE connectedVolEdge + PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) + IMPLICIT NONE + + INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB + LOGICAL:: coincident + INTEGER:: i + + coincident = .FALSE. + DO i = 1, 3 + IF (ANY(nodesA(i) == nodesB)) THEN + coincident = .TRUE. + + ELSE + coincident = .FALSE. + EXIT + + END IF + + END DO + + END FUNCTION coincidentNodes + SUBROUTINE connectedTetraTetra(elemA, elemB) IMPLICIT NONE @@ -217,54 +239,26 @@ MODULE moduleMesh3DCartRead !Check surface 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e1 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e4 => elemA @@ -275,54 +269,26 @@ MODULE moduleMesh3DCartRead !Check surface 2 IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e2 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e4 => elemA @@ -333,54 +299,26 @@ MODULE moduleMesh3DCartRead !Check surface 3 IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e3 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e4 => elemA @@ -391,54 +329,26 @@ MODULE moduleMesh3DCartRead !Check surface 4 IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e4 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) .OR. & - (elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e4 => elemA @@ -483,8 +393,12 @@ MODULE moduleMesh3DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF + END IF !Check surface 2 @@ -515,6 +429,9 @@ MODULE moduleMesh3DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF END IF @@ -547,6 +464,9 @@ MODULE moduleMesh3DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF END IF diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9a02cb1..e614cdf 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -230,12 +230,15 @@ MODULE moduleInject !Random position in edge partInj(n)%r = randomEdge%randPos() !Volume associated to the edge: - IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN + IF (ASSOCIATED(randomEdge%e1)) THEN partInj(n)%vol = randomEdge%e1%n - - ELSE + + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN partInj(n)%vol = randomEdge%e2%n + ELSE + PRINT *, "ERROR NO VOL ASSOCIATED TO EDGE" + END IF partInj(n)%v = (/ self%v(1)%obj%randomVel(), & diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index d6b10a5..c76b6c9 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -67,6 +67,13 @@ MODULE moduleSolver REAL(8):: tau, tauSp SELECT CASE(pusherType) + !3D Cartesian + CASE('3DCartNeutral') + self%pushParticle => push3DCartNeutral + + CASE('3DCartCharged') + self%pushParticle => push3DCartCharged + !2D Cylindrical CASE('2DCylNeutral') self%pushParticle => push2DCylNeutral