wip

fluid/fluid_problem.c

@@ -121,12 +121,10 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
     }
     double f[D],dfdu,dfddp;
     particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
-
     if (!local_vector)
       continue;
     double rho,nu;
     fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
-    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
     double f0[n_fields], f1[n_fields][D], f00[n_fields*n_fields],f11[n_fields*n_fields][D][D],f10[n_fields*n_fields][D],f01[n_fields*n_fields][D];
     for (int i = 0; i < n_fields; ++i) {
       f0[i] = 0;
@@ -145,12 +143,15 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
         }
       }
     }
+    double pspg = drag_in_stab?(problem->stability?problem->taup[iel]/rho:0):0;
     for (int d = 0; d < D; ++d) {
       f0[U+d] += f[d];
-      f1[P][d] += pspg*f[d];
-      for (int e = 0; e < D; ++e) {
-        double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
-        f1[U+d][e] += supg*f[d];
+      if (problem->stability){
+        f1[P][d] += pspg*f[d];
+        for (int e = 0; e < D; ++e) {
+          double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
+          f1[U+d][e] += supg*f[d];
+        }
       }
     }
     fe_fields_add_to_local_vector(problem->fields, f0, f1, sf, dsf, 1, local_vector);
@@ -159,12 +160,14 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
     for (int d = 0; d < D; ++d) {
       f00[(U+d)*n_fields+U+d] = dfdu;
       f01[(U+d)*n_fields+P][d] = dfddp;
-      f10[P*n_fields+U+d][d] = pspg*dfdu;
-      f11[P*n_fields+P][d][d] = pspg*dfddp;
-      for (int e = 0; e < D; ++e) {
-        double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
-        f10[(U+d)*n_fields+U+d][e] += supg*dfdu;
-        f11[(U+d)*n_fields+P][e][d] += supg*dfddp;
+      if(problem->stability){
+        f10[P*n_fields+U+d][d] = pspg*dfdu;
+        f11[P*n_fields+P][d][d] = pspg*dfddp;
+        for (int e = 0; e < D; ++e) {
+          double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
+          f10[(U+d)*n_fields+U+d][e] += supg*dfdu;
+          f11[(U+d)*n_fields+P][e][d] += supg*dfddp;
+        }
       }
     }
     fe_fields_add_to_local_matrix(problem->fields, f00, f01, f10, f11, sf, dsf, 1, local_matrix);
@@ -243,56 +246,135 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   }
 }
 
+
+#if 0
+static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double ds[][D], const double sold[D], const double dsold[][D], const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double f01[][D])
+{
+  const int vid = wbnd->vid;
+  const int pid = wbnd->pid;
+  const int cid = wbnd->cid;
+  const int aid = wbnd->aid;
+  double p = s[P];
+  const int n_fields = fluid_problem_n_fields(problem);
+  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;
+  for (int iD = 0; iD < D; ++iD) {
+    u[iD] = s[U+iD];
+    dp[iD] = ds[P][iD];
+    unold += sold[U+iD]*n[iD];
+    unmesh += mesh_velocity[iD]*n[iD];
+    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]*c);
+    unext += uext[iD]*n[iD];
+    dcn += dc[iD]*n[iD];
+    s_c += (u[iD]-uext[iD])*n[iD];
+  }
+  double h = problem->element_size[eid];
+  if (problem->stability){
+    if (wbnd->type == BND_WALL && pid >= 0) {
+        f0[P] = -(p-data[pid])/h*problem->taup[eid];
+        f00[P*n_fields+P] += -1/h*problem->taup[eid];
+    }
+  }
+  if(wbnd->type == BND_OPEN) {
+    f0[P] = unext;
+    for (int i = 0; i < D; ++i) {
+      f00[P*n_fields+U+i] += vid<0?n[i]:0;
+    }
+  }
+  for (int id = 0; id < D; ++id) {
+    f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
+    f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
+  }
+  double c_du_o_c[D][D];
+  double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
+  for (int iD = 0; iD < D; ++iD) {
+    for (int jD = 0; jD < D; ++jD)
+      c_du_o_c[iD][jD] = ds[U+iD][jD]-u[iD]/c*dc[jD];
+  }
+  if(wbnd->compute_viscous_term == 1){
+    for (int id = 0; id < D; ++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];
+        f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
+        f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
+        f01[(U+id)*n_fields+U+id][jd] -= mu*n[jd];
+        f01[(U+id)*n_fields+U+jd][id] -= mu*n[jd];
+      }
+    }
+  }
+  if (problem->advection) {
+    double unbnd = unold;
+    if (wbnd->type == BND_WALL) unbnd = unold/2;
+    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:uext[id])-(unold - unmesh)*u[id])*rho/c;
+        f00[(U+id)*n_fields+U+id] -= (unold-unmesh)*rho/c;
+      }
+    }
+  }
+}
+#endif
+
 static double pow2(double a) {return a*a;}
 void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
   const Mesh *mesh = problem->mesh;
-
   problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
-  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
-  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
-
-  double maxtaup = 0;
-  double mintaup = DBL_MAX;
-  const int n_fields = fluid_problem_n_fields(problem);
-  double maxtaa = 0;
-  double mintaa = DBL_MAX;
   for (int iel = 0; iel < mesh->n_elements; ++iel) {
-    double normu = {0.};
-    const int *el = &mesh->elements[iel*N_N];
-    double a = 0;
-    double amax = 0.5;
-    double amin = 0.5;
-    double c = 0;
-    for (int i=0; i< N_N; ++i) {
-      double normup = 0;
-      c += problem->porosity[el[i]]/N_N;
-      if(problem->n_fluids == 2){
-        a += problem->concentration[iel*N_N+i];
-        amax = fmax(amax,problem->concentration[iel*N_N+i]);
-        amin = fmin(amin,problem->concentration[iel*N_N+i]);
-      }
-      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);//to change
-      normu = fmax(normu,sqrt(normup));
-    }
     double dxidx[DIMENSION][DIMENSION];
     double detj = mesh_dxidx(mesh,iel,dxidx);
-#if DIMENSION == 2
-    const double h = 2*sqrt(detj/(2*M_PI));
-#else
-    const double h = 2*cbrt(detj/(8*M_PI));
-#endif
+  #if DIMENSION == 2
+      const double h = 2*sqrt(detj/(2*M_PI));
+  #else
+      const double h = 2*cbrt(detj/(8*M_PI));
+  #endif
     problem->element_size[iel] = h;
-    double rho, mu;
-    a /= N_N;
-    a = fmin(1.,fmax(0.,a));
-    fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
-    double nu = mu/rho;
-    problem->taup[iel] = 1/sqrt(
-        (problem->temporal?pow2(2/dt):0.)
-        +(problem->advection ?pow2(2*normu/h):0.)
-        +pow2(4*nu/pow2(h)));
-    problem->tauc[iel] = problem->advection?h*normu*fmin(h*normu/(6*nu),0.5):0.;
-    problem->taup[iel] = problem->tauc[iel] = 0;
+  }
+  if (problem->stability){
+    problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
+    problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
+    double maxtaup = 0;
+    double mintaup = DBL_MAX;
+    const int n_fields = fluid_problem_n_fields(problem);
+    double maxtaa = 0;
+    double mintaa = DBL_MAX;
+    for (int iel = 0; iel < mesh->n_elements; ++iel) {
+      double normu = {0.};
+      double h = problem->element_size[iel];
+      const int *el = &mesh->elements[iel*N_N];
+      double a = 0;
+      double amax = 0.5;
+      double amin = 0.5;
+      double c = 0;
+      for (int i=0; i< N_N; ++i) {
+        double normup = 0;
+        c += problem->porosity[el[i]]/N_N;
+        if(problem->n_fluids == 2){
+          a += problem->concentration[iel*N_N+i];
+          amax = fmax(amax,problem->concentration[iel*N_N+i]);
+          amin = fmin(amin,problem->concentration[iel*N_N+i]);
+        }
+        for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);//to change
+        normu = fmax(normu,sqrt(normup));
+      }
+      double rho, mu;
+      a /= N_N;
+      a = fmin(1.,fmax(0.,a));
+      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
+      double nu = mu/rho;
+      problem->taup[iel] = 1/sqrt(
+          (problem->temporal?pow2(2/dt):0.)
+          +(problem->advection?pow2(2*normu/h):0.)
+          +pow2(4*nu/pow2(h)));
+      problem->tauc[iel] = problem->advection?h*normu*fmin(h*normu/(6*nu),0.5):0.;
+    }
   }
 }
 
@@ -414,6 +496,8 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
   double *g = problem->g;
   double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
   double drag = problem->volume_drag + problem->quadratic_drag*nold;
+  // f0[P] = p;
+  // f00[P*n_fields+P] = 1;
   for (int i = 0; i < D; ++i) {
     f0[U+i] =
       + c*dp[i]
@@ -424,7 +508,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
     f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
     f01[(U+i)*n_fields+P][i] += c;
     if(problem->temporal) {
-      f0[U+i] += rho*(u[i]-uold[i])/dt;
+      f0[U+i] += rho*(u[i]-uold[i])/dt;// - rho;
       f00[(U+i)*n_fields+U+i] += rho/dt;
     }
     if (problem->advection) {
@@ -472,6 +556,115 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
   }
 }
 
+#if 0
+static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double f1[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {
+  double p = sol[P];
+  double taup = problem->stability?problem->taup[iel]:0;
+  double tauc = problem->stability?problem->tauc[iel]:0;
+  if(problem->stab_param != 0)
+    taup = tauc = 0;
+  const int n_fields = fluid_problem_n_fields(problem);
+  double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
+  for (int i = 0; i < n_fields; ++i) {
+    for (int j = 0; j < n_fields; ++j) {
+      int r = i*n_fields+j;
+      f00[r] = 0;
+      for (int id = 0; id < D; ++id) {
+        f01[r][id] = 0;
+        f10[r][id] = 0;
+        for (int jd = 0; jd < D; ++jd) {
+          f11[r][id][jd] = 0;
+        }
+      }
+    }
+  }
+  double divu = 0;
+  double nold = 0;
+  for (int iD = 0; iD < D; ++iD) {
+    u[iD] = sol[U+iD];
+    uold[iD] = solold[U+iD];
+    nold += uold[iD]*uold[iD];
+    umesh[iD] = mesh_velocity[iD];
+    dp[iD] = dsol[P][iD];
+    for (int jD = 0; jD < D; ++jD) {
+      du[iD][jD] = dsol[U+iD][jD];
+      duold[iD][jD] = dsolold[U+iD][jD];
+    }
+    divu += du[iD][iD];
+  }
+  nold = sqrt(nold);
+  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
+  //f00[P*n_fields+P] = 1/dt*0.1;
+  double *g = problem->g;
+  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
+  double drag = problem->volume_drag + problem->quadratic_drag*nold;
+  for (int i = 0; i < D; ++i) {
+    f1[P][i] = -u[i];
+    f10[P*n_fields+U+i][i] += -1;
+
+    f0[U+i] =
+      + c*dp[i]
+      - g[i]*rhoreduced*c
+      - bf[i]*c
+      + drag*u[i];
+      // + 5.3e5*u[i];
+    f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
+    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;
+    }
+    if (problem->advection) {
+      // advection :
+      // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
+      //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
+      for (int j = 0; j < D; ++j) {
+        f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
+        f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
+        f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
+        f01[(U+i)*n_fields+(U+i)][j] += rho/c*uold[j];
+        f01[(U+i)*n_fields+(U+i)][j] -= rho/c*umesh[j];
+      }
+    }
+    for (int j = 0; j < D; ++j) {
+      f1[U+i][j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
+      f10[(U+i)*n_fields+U+j][j] += -mu*dc[i]/c;
+      f10[(U+i)*n_fields+U+i][i] += -mu*dc[j]/c;
+      f11[(U+i)*n_fields+U+j][j][i] += mu;
+      f11[(U+i)*n_fields+U+i][j][j] += mu;
+      }
+    }
+    if (problem->stability){
+    for (int i = 0; i < D; ++i) {
+      // SUPG
+      for (int j = 0; j < D; ++j) {
+        double supg = uold[j]*taup;
+        f1[U+i][j] += supg*f0[U+i];
+        f10[(U+i)*n_fields+U+i][j] += supg*f00[(U+i)*n_fields+U+i];
+        f11[(U+i)*n_fields+P][j][i] += supg*f01[(U+i)*n_fields+P][i];
+        for (int k = 0; k < D; ++k)
+          f11[(U+i)*n_fields+(U+i)][j][k] += supg*f01[(U+i)*n_fields+(U+i)][k];
+      }
+      // least square incompressible constraint
+      double lsic = tauc*rho;
+      f1[U+i][i] += ((c-cold)/dt+divu)*lsic;
+      for (int j = 0; j < D; ++j) {
+        f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
+      }
+      // PSPG
+      double pspg = taup/rho;
+      f1[P][i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
+      f11[P*n_fields+P][i][i] += pspg*(f01[(U+i)*n_fields+P][i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
+      f10[P*n_fields+U+i][i] += pspg*f00[(U+i)*n_fields+(U+i)];
+      for (int j = 0; j < D; ++j) {
+        f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
+      }
+    }
+  }
+}
+#endif
+
 void weak_boundary_values(const FEFields *fields, const Mesh *mesh, const MeshBoundary *bnd, const WeakBoundary *wbnd, double *data);
 static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
 {
@@ -860,7 +1053,7 @@ static void fluid_problem_surface_tension(FluidProblem *problem, const double *s
     double c = 0;
     double kappa = 0;
     double uold[D] = {0};
-    double taup = problem->taup[iel];
+    double taup = problem->stability?problem->taup[iel]:0;
     for (int i = 0; i < N_N; ++i){
       c += problem->porosity[el[i]]/N_N;
       nda = 0;
@@ -903,12 +1096,14 @@ static void fluid_problem_surface_tension(FluidProblem *problem, const double *s
           //local_vector[(U+id)*N_N+iphi] += c*sigma*jw*(divnda[id]-divext_da[id])*phi[iphi];
           local_vector[(U+id)*element->nlocal+iphi] += -c*sigma*jw*kappa*da[id]*phi[iphi];
          // if (nda>1e-8) printf("tension surface = %.8e\n",c*10000000000000000*sigma*jw*(kappa[iphi])*grad_a_cg[el[iphi]*D+id]*phi[iphi]);
-          for (int jd = 0; jd < D; jd++){
-            double supg = uold[jd]*taup;
-            local_vector[(U+id)*element->nlocal+iphi] += dphi[iphi][jd]*c*sigma*jw*kappa*da[id]*supg;
+          if(problem->stability){
+            for (int jd = 0; jd < D; jd++){
+              double supg = uold[jd]*taup;
+              local_vector[(U+id)*element->nlocal+iphi] += dphi[iphi][jd]*c*sigma*jw*kappa*da[id]*supg;
+            }
+            double pspg = taup/rho;
+            local_vector[P*D+iphi] += jw*dphi[iphi][id]*c*sigma*kappa*da[id]*pspg;
           }
-          double pspg = taup/rho;
-          local_vector[P*D+iphi] += jw*dphi[iphi][id]*c*sigma*kappa*da[id]*pspg;
         }
       }
     }
@@ -1234,7 +1429,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
   }
 }
 
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f) {
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f, int stability) {
   if (n_fluids != 1 && n_fluids != 2) {
     printf("error : only 1 or 2 fluids are supported\n");
   }
@@ -1265,6 +1460,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
   problem->reduced_gravity = 0;
   problem->strong_boundaries = NULL;
   problem->n_strong_boundaries = 0;
+  problem->stability = stability;
   problem->temporal = temporal;
   problem->bulk_force = NULL;
   problem->advection = advection;
@@ -1278,7 +1474,6 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
     problem->rho[i] = rho[i];
     problem->mu[i] = mu[i];
   }
-  problem->coeffStab = coeffStab;
   for (int i = 0; i < D; ++i){
     problem->g[i] = g[i];
   }
@@ -1338,8 +1533,10 @@ void fluid_problem_free(FluidProblem *problem)
   free(problem->particle_volume);
   free(problem->particle_velocity);
   free(problem->particle_contact);
-  free(problem->taup);
-  free(problem->tauc);
+  if(problem->taup) 
+    free(problem->taup);
+  if(problem->tauc)
+     free(problem->tauc);
   free(problem->element_size);
   free(problem->grad_a_cg);
   for (int i = 0; i < problem->n_strong_boundaries; ++i)
@@ -1634,41 +1831,22 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
   problem->bulk_force = malloc(sizeof(double)*mesh->n_nodes*DIMENSION);
   for (int i=0; i<mesh->n_nodes*DIMENSION; ++i)
     problem->bulk_force[i] = 0;
-  // free(problem->porosity);
-  // problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
   free(problem->grad_a_cg);
   problem->grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
   problem->all_kappa = malloc(sizeof(double)*mesh->n_elements);
-
-  // free(problem->oldporosity);
-  // problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
-  // for(int i = 0; i < mesh->n_nodes; ++i){
-    // problem->porosity[i] = 1.;
-    // problem->oldporosity[i] = 1.;
-  // }
-
   free(problem->boundary_force);
-  problem->boundary_force = (double*)malloc(DIMENSION*mesh->n_boundaries*sizeof(double));
+  problem->boundary_force = malloc(DIMENSION*mesh->n_boundaries*sizeof(double));
   for(int i = 0; i < mesh->n_boundaries*DIMENSION; i++) problem->boundary_force[i] = 0;
-
-
-  // free(problem->solution);
-  // int n_fields = fluid_problem_n_fields(problem);
-  // int ndof = fluid_problem_get_ndof(problem->fields, problem->mesh);
-  // problem->solution = calloc(ndof, sizeof(double));
-  // problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));// to change
   if (problem->n_fluids == 2) {
     free(problem->concentration);
-    problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double));
+    problem->concentration = malloc(mesh->n_elements*N_N*sizeof(double));
   }
   free(problem->mesh_velocity);
-  problem->mesh_velocity = (double*)calloc(mesh->n_nodes*D,sizeof(double));
-  // for(int i = 0; i < mesh->n_nodes*D; ++i) problem->mesh_velocity[i] = 0;
-  problem->taup = malloc(sizeof(double)*mesh->n_elements);
-  for (int i = 0; i < problem->mesh->n_elements; ++i) problem->taup[i] = 0;
-  // for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){ // to change
-    // problem->solution[i] = 0.;
-// }
+  problem->mesh_velocity = calloc(mesh->n_nodes*D,sizeof(double));
+  if(problem->stability){
+    problem->taup = calloc(mesh->n_elements, sizeof(double));
+    problem->tauc = calloc(mesh->n_elements, sizeof(double));
+  }
   fluid_problem_compute_node_volume(problem);
   if (problem->boundaries) mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
   problem->boundaries = mesh_boundaries_create(problem->mesh);

fluid/fluid_problem.h

@@ -92,6 +92,7 @@ struct FluidProblem {
   double *all_kappa;
   int n_fluids;
   //stabilisation coefficients
+  int stability;
   double *taup;
   double *tauc;
   int *particle_element_id;
@@ -111,7 +112,7 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
 void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
 void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
 void fluid_problem_free(FluidProblem *problem);
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f);
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f, int stability);
 void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume);
 void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh);
 int *fluid_problem_particle_element_id(FluidProblem *problem);

paraview_utils/MigFlow.xml

@@ -63,7 +63,7 @@
        name="Script"
        command="SetScript"
        number_of_elements="1"
-       default_values="import numpy as np
in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1) + np.count_nonzero(in_periodic.CellTypes == 1)
n_segments  = (np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3) 
            + np.count_nonzero(in_periodic.CellTypes == 7) + np.count_nonzero(in_periodic.CellTypes == 3))
n_triangles = (np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
            + np.count_nonzero(in_periodic.CellTypes == 13)+np.count_nonzero(in_periodic.CellTypes==5))

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
    vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    vals_t.append(in_contacts.FieldData["particle_disk_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_disk_s"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    vals_t.append(in_contacts.FieldData["particle_segment_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_segment_s"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    vals_t.append(in_contacts.FieldData["particle_triangle_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_triangle_s"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    d0 = in_bnd.Points[triangles[contacts[:,0],:]][:,1,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
    d1 = in_bnd.Points[triangles[contacts[:,0],:]][:,2,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
    N = np.cross(d0,d1)
    N /= np.linalg.norm(N,axis=1)[:,newaxis]
    dd = in_bnd.Points[triangles[contacts[:,0],:]][:,0,:] - in_particles.Points[contacts[:,1],:]
    dist = einsum('ij,ij->i', N, dd)
    points_t[:,1,:] = in_particles.Points[contacts[:,1],:] + N*dist[:,newaxis]
    points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
    vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
    +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
    output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
+       default_values="import numpy as np
in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1) + np.count_nonzero(in_periodic.CellTypes == 1)
n_segments  = (np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3) 
            + np.count_nonzero(in_periodic.CellTypes == 7) + np.count_nonzero(in_periodic.CellTypes == 3))
n_triangles = (np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
            + np.count_nonzero(in_periodic.CellTypes == 13)+np.count_nonzero(in_periodic.CellTypes==5))

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
    vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    vals_t.append(in_contacts.FieldData["particle_disk_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_disk_s"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    vals_t.append(in_contacts.FieldData["particle_segment_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_segment_s"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    vals_t.append(in_contacts.FieldData["particle_triangle_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_triangle_s"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    d0 = in_bnd.Points[triangles[contacts[:,0],:]][:,1,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
    d1 = in_bnd.Points[triangles[contacts[:,0],:]][:,2,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
    N = np.cross(d0,d1)
    N /= np.linalg.norm(N,axis=1)[:,newaxis]
    dd = in_bnd.Points[triangles[contacts[:,0],:]][:,0,:] - in_particles.Points[contacts[:,1],:]
    dist = einsum('ij,ij->i', N, dd)
    points_t[:,1,:] = in_particles.Points[contacts[:,1],:] + N*dist[:,newaxis]
    points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
    vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
    +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
    output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.uint32)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
        panel_visibility="advanced">
        <Hints>
         <Widget type="multi_line"/>

paraview_utils/contacts-extract.py

@@ -130,7 +130,7 @@ def RequestData():
     if vals_s is not None :
         output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
     n = points.shape[0]
-    pattern = np.ndarray([nf,4], np.int)
+    pattern = np.ndarray([nf,4], np.uint32)
     for i in range(nf) :
         j = (i+1)%nf
         pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)

python/fluid.py

@@ -126,8 +126,7 @@ def _load_msh(mesh_file_name, lib, dim) :
 
 class FluidProblem :
     """Creates the numerical representation of the fluid."""
-
-    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., solver=None, solver_options="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True, order=[1,1,1]):
+    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., solver=None, solver_options="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True, order=[1,1,1], stability=True):
         """Builds the fluid structure.
 
         Keyword arguments:
@@ -171,7 +170,8 @@ class FluidProblem :
         self._fp = c_void_p(self._lib.fluid_problem_new(
             _np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), 
             c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),
-            c_double(ip_factor), c_int(temporal), c_int(advection),_np2c(np.array(self.order_fields, dtype=np.int32), dtype=np.int32)
+            c_double(ip_factor), c_int(temporal), c_int(advection),
+            _np2c(np.array(self.order_fields, dtype=np.int32), dtype=np.int32), c_int(stability)
             ))
         if self._fp == None :
             raise NameError("Cannot create fluid problem.")
@@ -193,11 +193,11 @@ class FluidProblem :
         """
         mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
         self._lib.fluid_problem_set_mesh(self._fp, mesh)
-        idx = self.set_numbering([1]*self.dimension())
+        idx = self.get_numbering([1]*self.dimension())
         self._lib.fluid_problem_mesh_local_map(self._fp, _np2c(idx, dtype=np.int32))
-        idx = self.set_numbering(self.order_fields)
-        idx_concentration = self.set_numbering(self.order_concentration)
-        idx_porosity = self.set_numbering(self.order_porosity)
+        idx = self.get_numbering(self.order_fields)
+        idx_concentration = self.get_numbering(self.order_concentration)
+        idx_porosity = self.get_numbering(self.order_porosity)
         self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
         self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
         self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
@@ -378,7 +378,7 @@ class FluidProblem :
         isp2p1 = np.max(self.order_fields) == 2
         if isp2p1 :
             field_data.append((b"velocity", self.velocity()))
-            idx = self.set_numbering([2])
+            idx = self.get_numbering([2])
             cell_data.append(("velocity_map", idx.reshape(self.n_elements(),-1)))
         else:
             data.append(("velocity", self.velocity()))
@@ -407,7 +407,7 @@ class FluidProblem :
                 _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
                 _np2c(data["parent_node_id"],np.int32) if "parent_node_id"  in data else None))
         self._lib.fluid_problem_set_mesh(self._fp, _mesh)
-        ndof = np.max(self.set_numbering(self.order_fields))+1
+        ndof = np.max(self.get_numbering(self.order_fields))+1
         id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
             + self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
             np.arange(self.n_fields()*self.n_nodes(), ndof)]
@@ -458,9 +458,9 @@ class FluidProblem :
         self._lib.fluid_problem_apply_boundary_conditions(self._fp)
         self._lib.fluid_problem_reset_boundary_force(self._fp)
         self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
-        idx = self.set_numbering(self.order_fields)
-        idx_concentration = self.set_numbering(self.order_concentration)
-        idx_porosity = self.set_numbering(self.order_porosity)
+        idx = self.get_numbering(self.order_fields)
+        idx_concentration = self.get_numbering(self.order_concentration)
+        idx_porosity = self.get_numbering(self.order_porosity)
         self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
         self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
         self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
@@ -559,12 +559,12 @@ class FluidProblem :
 
     def solution(self) :
         """Gives access to the fluid field value at each degree of freedom"""