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

Slight improvement in using triangles in electrostaic pusher, but still

Browse source 
no final solution.
This commit is contained in:
Jorge Gonzalez 2020-11-29 12:09:54 +01:00
commit bf5310c2c3
1 changed files with 23 additions and 37 deletions
Download patch file Download diff file Expand all files Collapse all files

60
src/modules/moduleMeshCyl.f95
View file

@ -108,6 +108,8 @@ MODULE moduleMeshCyl
!Connectivity to adjacent elements
CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
REAL(8):: arNodes(1:3) = 0.D0
!Derivatives in z,r real space
REAL(8), DIMENSION(1:3):: dPsiZ, dPsiR
CONTAINS
PROCEDURE, PASS:: init => initVolTriaCyl
@ -618,39 +620,16 @@ MODULE moduleMeshCyl
CLASS(meshVolCylTria), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8):: ldCorner(1:3) !left down (ld) corner of square enclousing the triangle
INTEGER:: minP
INTEGER:: pOrder(1:3)
REAL(8), DIMENSION(1:3):: r1, r2, r3
REAL(8):: A
!Assign node index
self%n = n
!Find the first node
r1 = mesh%nodes(p(1))%obj%getCoordinates()
r2 = mesh%nodes(p(2))%obj%getCoordinates()
r3 = mesh%nodes(p(3))%obj%getCoordinates()
ldCorner = (/ MINVAL((/r1(1), r2(1), r3(1)/), 1), &
MINVAL((/r1(2), r2(2), r3(2)/), 1), &
0.D0 /)
minP = MINLOC((/ NORM2(r1-ldCorner), NORM2(r2-ldCorner), NORM2(r3-ldCorner) /), 1)
pOrder(1) = p(minP)
SELECT CASE(minP)
CASE(1)
pOrder(2) = p(2)
pOrder(3) = p(3)
CASE(2)
pOrder(2) = p(3)
pOrder(3) = p(1)
CASE(3)
pOrder(2) = p(1)
pOrder(3) = p(2)
END SELECT
!Assign nodes to element
self%n1 => mesh%nodes(pOrder(1))%obj
self%n2 => mesh%nodes(pOrder(2))%obj
self%n3 => mesh%nodes(pOrder(3))%obj
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
self%n3 => mesh%nodes(p(3))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
@ -663,6 +642,20 @@ MODULE moduleMeshCyl
self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3)
!Derivatives in z/r for shape functions (node independent)
!TODO: This is used because invJ.dPsi does not produce the right Electric field
A = self%z(2)*self%r(3) - self%z(3)*self%r(2) + &
self%z(3)*self%r(1) - self%z(1)*self%r(3) + &
self%z(1)*self%r(2) - self%z(2)*self%r(1)
self%dPsiZ = (/ self%r(2)-self%r(3), &
self%r(3)-self%r(1), &
self%r(1)-self%r(2) /)/A
self%dPsiR = (/ self%z(3)-self%z(2), &
self%z(1)-self%z(3), &
self%z(2)-self%z(1) /)/A
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
@ -882,9 +875,6 @@ MODULE moduleMeshCyl
CLASS(meshVolCylTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:2,1:3)
REAL(8):: dPsiR(1:2,1:3)
REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: phi(1:3)
REAL(8):: EF(1:3)
@ -892,12 +882,8 @@ MODULE moduleMeshCyl
self%n2%emData%phi, &
self%n3%emData%phi/)
dPsi = self%dPsi(xi)
detJ = self%detJac(xi,dPsi)
invJ = self%invJac(xi,dPsi)
dPsiR = MATMUL(invJ, dPsi)/detJ
EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
EF(1) = -DOT_PRODUCT(self%dPsiZ(:), phi)
EF(2) = -DOT_PRODUCT(self%dPsiR(:), phi)
EF(3) = 0.D0
END FUNCTION gatherEFTria
@ -926,7 +912,7 @@ MODULE moduleMeshCyl
REAL(8):: dPsi(1:2,1:3)
!Direct method to convert coordinates
xi = (/ 0.D0, 0.D0, 0.D0 /) !Irrelevant, required for input
xi = 0.D0 !Irrelevant, required for input
deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)