chngement inertie

marblesbag/src/fluid_problem.c

@@ -629,12 +629,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
           double utau = 0;
           double dphisigma=0;
           for (int k = 0; k < DIMENSION; ++k) {
-            utau += u[k]*tau[j][k];
+            utau += u[k]*dphii[k]*u[j]/c;
             dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
           }
           local_vector[iphi+N_SF*j] += jw*(
               dphisigma
-              +rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
+              +rho*phii*dudt[j]-rho*utau
               -p*dphii[j]-fu[j]*phii
               );
         }
@@ -660,7 +660,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
           for (int j = 0; j < DIMENSION; ++j) {
             dtau[j] = dphij[j]-phij*dc[j]/c;
             dphiidphij += dphii[j]*dphij[j];
-            udtau += u[j]*dtau[j];
+            udtau += u[j]*phij*dphii[j]/c;
             dphiidtau += dphii[j]*dtau[j];
           }
           for (int j = 0; j < DIMENSION; ++j) {
@@ -669,9 +669,9 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
             LOCAL_MATRIX(P,U) +=  jw*phii*dphij[j];
             for (int k = 0; k < DIMENSION; ++k) {
               int V = k;
-              LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
+              LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
             }
-            LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
+            LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
           }
           LOCAL_MATRIX(P,P) +=  jw*epsilon[el[iphi]]*dphiidphij;
         }