integration by part on pressure term

fluid/fluid_problem.c

@@ -233,6 +233,8 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
         f01[(U+id)*n_fields+U+id][jd] -= mu*n[jd];
         f01[(U+id)*n_fields+U+jd][id] -= mu*n[jd];
       }
+      f0[U+id] += c*p*n[id]; 
+      f00[(U+id)*n_fields+P] += c*n[id];
     }
   }
   if (problem->advection) {
@@ -418,13 +420,13 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
     f11[P*n_fields+P][i][i] += stab_param;
 
     f0[U+i] =
-      + c*dp[i]
       - g[i]*rhoreduced*c
       - bf[i]*c
       + drag*u[i];
+      // + c*dp[i]
       // + 5.3e5*u[i];
     f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
-    f01[(U+i)*n_fields+P][i] += c;
+    // f01[(U+i)*n_fields+P][i] += c;
     if(problem->temporal) {
       f0[U+i] += rho*(u[i]-uold[i])/dt;
       f00[(U+i)*n_fields+U+i] += rho/dt;
@@ -448,6 +450,8 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
       f11[(U+i)*n_fields+U+j][j][i] += mu;
       f11[(U+i)*n_fields+U+i][j][j] += mu;
       }
+    f1[U+i][i] -= c*p;
+    f10[(U+i)*n_fields+P][i] = -c;
     }
     if (problem->stability){
     for (int i = 0; i < D; ++i) {