better arrays

fluid/fluid_problem.c

@@ -133,8 +133,8 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
         dc[k] += problem->porosity[el[i]]*dphi[i][k];
       }
     }
-    double f[D],dfdu[D][D],dfddp[D][D];
-    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
+    double f[D],dfdu[D*D],dfddp[D*D];
+    particle_force_f(problem,f,dfdu,dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
     if (!local_vector)
       continue;
     double rho,nu;
@@ -156,14 +156,14 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
       for (int iphi = 0; iphi < N_SF; ++iphi){
         for (int jphi = 0; jphi < N_SF; ++jphi){
             for(int jD = 0; jD < D; jD++){
-              LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += phi[iphi]*phi[jphi]*dfdu[iD][jD];
-              LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][jD]*dfddp[iD][jD];
-              LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu[iD][jD];
-              LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][jD]*dfddp[iD][jD];    
+              LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += phi[iphi]*phi[jphi]*dfdu[iD*D+jD];
+              LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][jD]*dfddp[iD*D+jD];
+              LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu[iD*D+jD];
+              LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][jD]*dfddp[iD*D+jD];    
           for (int kD = 0; kD < D; ++kD) {
             double supg = drag_in_stab?sold[U+kD]*problem->taup[iel]:0;
-            LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += supg*dphi[iphi][kD]*phi[jphi]*dfdu[iD][jD];
-            LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][kD]*dphi[jphi][jD]*dfddp[iD][jD];
+            LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += supg*dphi[iphi][kD]*phi[jphi]*dfdu[iD*D+jD]; 
+            LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][kD]*dphi[jphi][jD]*dfddp[iD*D+jD];
           }
         }
         }
@@ -363,7 +363,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
   }
   else{
     gamma = particle_drag_coeff(Due,mu,rho,vol,c);
-  }
+  } 
   gamma = gamma*problem->drag_coeff_factor;
   double mat[D][D] = {[0][0] = 1, [1][1] = 1};
 #if DIMENSION==3