frontiere : vitesse fluide

fluid/fluid_problem.c

@@ -181,9 +181,10 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   const int aid = wbnd->aid;
   double p = s[P];
   const int n_fields = fluid_problem_n_fields(problem);
-  double u[D], dp[D];
+  double u[D], dp[D], uext[D];
   double s_c = 0;
   double unold = 0;
+  double cext = cid<0?c:data[cid];
   double unext = 0;
   double unmesh = 0;
   double dcn = 0;
@@ -192,9 +193,10 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     dp[iD] = ds[P*D+iD];
     unold += sold[U+iD]*n[iD];
     unmesh += mesh_velocity[iD]*n[iD];
-    unext += (vid<0?s[U+iD]:data[vid+iD])*n[iD];
+    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]);
+    unext += uext[iD]*n[iD];
     dcn += dc[iD]*n[iD];
-    s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
+    s_c += (u[iD]-uext[iD])*n[iD];
   }
   double h = problem->element_size[eid];
   if (wbnd->type == BND_WALL && pid >= 0) {
@@ -219,7 +221,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   }
   if(wbnd->type != BND_SYMMETRY && wbnd->compute_viscous_term == 1){
     for (int id = 0; id < D; ++id) {
-      f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
+      f0[U+id] += sigma*(u[id]-uext[id]+s_c*n[id]);
       f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
       for (int jd = 0; jd <D; ++jd) {
         f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
@@ -236,7 +238,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     else if(vid>0) unbnd = (unold+unext)/2;
     if (unbnd - unmesh < 0) {
       for (int id = 0; id < D; ++id) {
-        f0[U+id] += ((unbnd - unmesh)*(vid<0?0:data[vid+id])-(unold - unmesh)*u[id])*rho/c;
+        f0[U+id] += ((unbnd - unmesh)*(vid<0?0:uext[id])-(unold - unmesh)*u[id])*rho/c;
         f00[(U+id)*n_fields+U+id] -= unold*rho/c;
       }
     }

scontact/scontact.c

@@ -227,6 +227,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
   c->type = PARTICLE_DISK;
   return c->D < alert;
 }
+
 struct _PeriodicSegment{
   size_t entity_id;
   double p[2][DIMENSION];