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

Fixing some issues in 3D Cartesian coordinates. Included 3D pusher.

Browse source 
Still there are issues linking a volume to an edge.
This commit is contained in:
Jorge Gonzalez 2021-03-15 10:00:34 +01:00
commit c236c5e0e2
8 changed files with 123 additions and 148 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/modules/mesh/1DCart/moduleMesh1DCartRead.f90
View file

@ -218,6 +218,9 @@ MODULE moduleMesh1DCartRead
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. &

3
src/modules/mesh/1DRad/moduleMesh1DRadRead.f90
View file

@ -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. &

21
src/modules/mesh/2DCart/moduleMesh2DCartRead.f90
View file

@ -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

18
src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
View file

@ -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

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

@ -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()

208
src/modules/mesh/3DCart/moduleMesh3DCartRead.f90
View file

@ -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

7
src/modules/moduleInject.f90
View file

@ -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(), &

7
src/modules/moduleSolver.f90
View file

@ -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