From 31483b8b5082d80ead4fc456d76a61d1afe8999f Mon Sep 17 00:00:00 2001
From: Cyrille <91067824+cmosbeux@users.noreply.github.com>
Date: Thu, 20 Jul 2023 16:18:43 +0200
Subject: [PATCH] Update ParticleUtils.F90 with correct permutation of Mesh
 velocities

I noticed a permutation error when considering a 3D stokes Flow Solution, the solution is 4D (vx,vy,vz,p).

The changes allow to account for Flow Solution that have more than the velocity direction and should not break anything.
---
 fem/src/ParticleUtils.F90 | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/fem/src/ParticleUtils.F90 b/fem/src/ParticleUtils.F90
index bdb07a7c17..df08e2b94f 100644
--- a/fem/src/ParticleUtils.F90
+++ b/fem/src/ParticleUtils.F90
@@ -4631,6 +4631,9 @@ SUBROUTINE GetVectorFieldInMesh(Var, CurrentElement, Basis, Velo, dBasisdx, Grad
         
     
     IF( npos == 0 ) RETURN
+
+    ! MAX 3 velocity direction (4th variable is pressure)
+    ! This has to be accounted for in the Var % Values permutations
     dofs = MIN(3,Var % Dofs)
     
     !-----------------------------------------------------------------
@@ -4644,7 +4647,9 @@ SUBROUTINE GetVectorFieldInMesh(Var, CurrentElement, Basis, Velo, dBasisdx, Grad
       DO i=1,n
         j = LocalPerm(i)
 	DO k=1,dofs
-          LocalVelo(i,k) = Var % Values( Dofs*(j-1)+k)
+          ! For correct permutation we have to account for the full dimension
+          ! of the Flow Solution (e.g., Stokes 3D : vx, vy,vz, p)
+          LocalVelo(i,k) = Var % Values( Var % Dofs * (j-1) + k)
         END DO
       END DO
     ELSE    
@@ -4656,7 +4661,7 @@ SUBROUTINE GetVectorFieldInMesh(Var, CurrentElement, Basis, Velo, dBasisdx, Grad
         IF( j > 0 ) THEN
           SumBasis = SumBasis + Basis(i)
           DO k=1,dofs
-            LocalVelo(i,k) = Var % Values( Dofs*(j-1)+k)
+            LocalVelo(i,k) = Var % Values( Var % Dofs * (j-1) + k)
           END DO
         ELSE
           Basis(i) = 0.0_dp