accès aux iterations du solveur fluide

marblesbag/fluid.py

@@ -54,7 +54,7 @@ class fluid_problem :
         return forces
 
     def implicit_euler(self, dt) :
-        lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
+        return lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
 
     def set_particles(self, mass, volume, position, velocity):
         def np2c(a) :

marblesbag/fluid3.py

@@ -57,7 +57,7 @@ class fluid_problem :
         return forces
 
     def implicit_euler(self, dt) :
-        lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
+        return lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
 
     def set_particles(self, mass, volume, position, velocity, elid=None):
         def np2c(a) :

marblesbag/scontact/scontact.c

@@ -14,6 +14,19 @@ typedef struct _Disk Disk;
 typedef struct _Segment Segment;
 typedef struct _Triangle Triangle;
 
+static int check_word_in_file(FILE *f, const char *w) {
+  char word[256];
+  if(fscanf(f,"%255s", word) != 1){
+    printf("Invalid mesh file (\"%s\" expected).\n", w);
+    return -1;
+  }
+  if(strcmp(word, w) != 0){
+    printf("Invalid mesh file (\"%s\" expected).\n", w);
+    return -1;
+  }
+  return 0;
+}
+
 struct _ParticleProblem{
   int jacobi;
   double relax;
@@ -990,10 +1003,10 @@ static int particleProblemWriteContactVtp(ParticleProblem *p, const char *dirnam
     printf("cannot open file %s\n", filename);
     return -1;
   }
-  size_t n_contacts[DIMENSION+2] = {0};
-  for (size_t i = 0; i < vectorSize(p->contacts); ++i) {
-    n_contacts[p->contacts[i].type]++;
-  }
+/*  size_t n_contacts[DIMENSION+2] = {0};*/
+/*  for (size_t i = 0; i < vectorSize(p->contacts); ++i) {*/
+/*    n_contacts[p->contacts[i].type]++;*/
+/*  }*/
   fprintf(f, "<VTKFile type=\"PolyData\">\n");
   fprintf(f, "<PolyData>\n");
   fprintf(f, "<FieldData>\n");
@@ -1039,6 +1052,11 @@ static int particleProblemWriteParticleVtp(ParticleProblem *p, const char *dirna
     fprintf(f, "%.16g ", p->particles[i].r);
   }
   fprintf(f, "\n</DataArray>\n");
+  fprintf(f, "<DataArray Name=\"Mass\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
+  for (size_t i = 0; i < vectorSize(p->particles); ++i) {
+    fprintf(f, "%.16g ", p->particles[i].m);
+  }
+  fprintf(f, "\n</DataArray>\n");
   fprintf(f, "<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float32\" format=\"ascii\">\n");
   for (size_t i = 0; i < vectorSize(p->particles); ++i) {
 #if DIMENSION == 2
@@ -1063,6 +1081,7 @@ static int particleProblemWriteParticleVtp(ParticleProblem *p, const char *dirna
   fprintf(f, "</PolyData>\n");
   fprintf(f, "</VTKFile>\n");
   fclose(f);
+  free(filename);
   return 0;
 }
 
@@ -1077,7 +1096,6 @@ static int particleProblemWriteBndVtu(ParticleProblem *p, const char *dirname, i
   size_t nt = 0;
 #endif
   FILE *f = fopen(filename, "w");
-  f = fopen(filename, "w");
   fprintf(f, "<VTKFile type=\"UnstructuredGrid\">\n");
   fprintf(f, "<UnstructuredGrid>\n");
   fprintf(f, "<FieldData>\n");
@@ -1179,6 +1197,226 @@ int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id) {
   return 0;
 }
 
+
+
+
+
+static int particleProblemReadBndVtu(ParticleProblem *p, const char *dirname, int id) {
+  char *filename = malloc((strlen(dirname)+50)*sizeof(char));
+  sprintf(filename,"%s/boundaries_%05d.vtu",dirname, id);
+  size_t nd;
+  size_t ns;
+  size_t nt = 0;
+  size_t scratch;
+  FILE *f = fopen(filename, "r");
+  check_word_in_file(f, "<VTKFile type=\"UnstructuredGrid\">\n");
+  check_word_in_file(f, "<UnstructuredGrid>\n");
+  check_word_in_file(f, "<FieldData>\n");
+  check_word_in_file(f, "<DataArray Name=\"sizes\" NumberOfTuples=\"3\" NumberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n");
+  if (fscanf(f, "%lu %lu %lu\n", &nd, &ns, &scratch)!=3){
+    printf("Cannot read number of disks and segments\n");
+    return 1;
+  }
+  check_word_in_file(f, "</DataArray>\n");
+  check_word_in_file(f, "</FieldData>\n");
+  if (fscanf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"%lu\">\n", &scratch,&scratch)!=2){
+    printf("Cannot read number of points/cells\n");
+    return 1;
+  }
+  nt = scratch-nd-ns;
+  check_word_in_file(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
+  double scratchd;
+  for (size_t i = 0; i < nd; ++i) {
+    Disk *d = vectorPush(&p->disks);
+#if DIMENSION == 2
+    if (fscanf(f, "%lf %lf %lf ", &d->x[0], &d->x[1], &scratchd)!=3){
+      printf("Cannot read disks\n");
+      return 1;
+    }
+#else
+    if (fscanf(f, "%lf %lf %lf ", &d->x[0], &d->x[1], &d->x[2])!=3){
+      printf("Cannot read disks\n");
+      return 1;
+    }
+#endif
+  }
+  for (size_t i = 0; i < ns; ++i) {
+    Segment *s = vectorPush(&p->segments);
+#if DIMENSION == 2
+    if (fscanf(f, "%lf %lf %lf ", &s->p[0][0], &s->p[0][1],  &scratchd)!=3){
+       printf("Cannot read segments\n");
+       return 1;
+    }
+#else
+    if (fscanf(f, "%lf %lf %lf ", &s->p[0][0], &s->p[0][1], &s->p[0][2])!=3){
+       printf("Cannot read segments\n");
+       return 1;
+    }
+#endif
+    if (fscanf(f, "%lf %lf %lf ", &s->p[1][0], &s->p[1][1], DIMENSION == 2 ? &scratchd : &s->p[1][2])!=3){
+       printf("Cannot read segments\n");
+       return 1;
+    }
+  }
+#if DIMENSION == 3
+  for (int i = 0; i < nt; ++i) {
+    Triangle *t = vectorPush(&p->triangles);
+    if (fscanf(f, "%.16g %.16g %.16g ", &(t->p[0][0]), &(t->p[0][1]), &(t->p[0][2]))!=3||
+    fscanf(f, "%.16g %.16g %.16g ", &(t->p[1][0]), &(t->p[1][1]), &(t->p[1][2]))!=3||
+    fscanf(f, "%.16g %.16g %.16g ", &(t->p[2][0]), &(t->p[2][1]), &(t->p[2][2]))!=3){
+      printf("Cannot read triangles\n");
+      return 1;
+    }
+  }
+#endif
+  fclose(f);
+  free(filename);
+  return 0;
+}
+
+
+
+
+
+static int particleProblemReadParticleVtp(ParticleProblem *p, const char *dirname, int id) {
+  char *filename = malloc((strlen(dirname)+50)*sizeof(char));
+  sprintf(filename,"%s/particles_%05d.vtp",dirname, id);
+  FILE *f = fopen(filename, "r");
+  if (!f){
+    printf("cannot open file %s\n", filename);
+    return -1;
+  }
+  size_t n_pp_contacts = 0;
+  size_t n_p;
+  check_word_in_file(f, "<VTKFile type=\"PolyData\">\n");
+  check_word_in_file(f, "<PolyData>\n");
+  if (fscanf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"0\">\n", &n_p)!=1){
+    printf("Cannot read number of particles\n");
+    return 1;
+  }
+  check_word_in_file(f, "<PointData Scalars=\"Radius\" Vectors=\"Velocity\">\n");
+  check_word_in_file(f, "<DataArray Name=\"Radius\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
+  for (size_t i = 0; i < n_p; ++i) {
+    Particle *particle = vectorPush(&p->particles);
+    if (fscanf(f, "%lf ", &particle->r)!=1){
+      printf("Cannot read radii\n");
+      return 1;
+    }
+  }
+  check_word_in_file(f, "\n</DataArray>\n");
+  check_word_in_file(f, "<DataArray Name=\"Mass\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
+  for (size_t i = 0; i < n_p; ++i) {
+    if (fscanf(f, "%lf ", &p->particles[i].m)!=1){
+      printf("Cannot read radii\n");
+      return 1;
+    }
+  }
+  check_word_in_file(f, "\n</DataArray>\n");
+  check_word_in_file(f, "<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float32\" format=\"ascii\">\n");
+  for (size_t i = 0; i < n_p; ++i) {
+#if DIMENSION == 2
+    vectorPushN(&p->velocity, DIMENSION);
+    if (fscanf(f, "%lf %lf 0 ", &p->velocity[i * 2], &p->velocity[i * 2 + 1])!=2){
+      printf("Cannot read particles velocities\n");
+      return 1;
+    }
+#else
+    vectorPushN(&p->velocity, DIMENSION);
+    if (fscanf(f, "%lf %lf %lf ", &(p->velocity[i * 3]), &(p->velocity[i * 3 + 1]), &(p->velocity[i * 3 + 2]))!=3){
+      printf("Cannot read particle velocities\n");
+      return 1;
+    }
+#endif
+  }
+  check_word_in_file(f, "\n</DataArray>\n");
+  check_word_in_file(f, "</PointData>\n");
+  check_word_in_file(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
+  for (size_t i = 0; i < n_p; ++i) {
+#if DIMENSION == 2
+    vectorPushN(&p->position, DIMENSION);
+    if (fscanf(f, "%lf %lf 0 ", &(p->position[i * 2]), &(p->position[i * 2 + 1]))!=2){
+      printf("Cannot read particles positions\n");
+      return 1;
+    }
+#else
+    vectorPushN(&p->position, DIMENSION);
+    if (fscanf(f, "%lf %lf %lf ", &(p->position[i * 3]), &(p->position[i * 3 + 1]), &(p->position[i * 3 + 2]))!=3){
+      printf("Cannot read particles positions\n");
+      return 1;
+    }
+#endif
+  }
+  fclose(f);
+  return 0;
+}
+
+
+
+static int _readContactVtkSingle(ParticleProblem *p, FILE *f, int ctype)
+{
+  size_t n = 0;
+  char scratch;
+  if (fscanf(f, "<DataArray Name=\"%s\" NumberOfTuples=\"%zu\" NumberOfComponents=\"1\" type=\"Float32\" format=\"ascii\">\n",&scratch, &n)!=2){
+     printf("Cannot read number of contacts for impulsion\n");
+     return 1;
+  }
+  for (size_t i = 0; i < n; ++i) {
+    Contact *c =  vectorPush(&p->contacts);
+    c[i].type = ctype;
+    if (fscanf(f, "%lf ", &(c[i].dv))!=1){
+      printf("Cannot read impulsion\n");
+      return 1;
+    }
+  }
+  check_word_in_file(f, "\n</DataArray>\n");
+  if (fscanf(f, "<DataArray Name=\"%s_idx\" NumberOfTuples=\"%zu\" NumberOfComponents=\"2\" type=\"UInt64\" format=\"ascii\">\n",&scratch,&n)!=2){
+    printf("Cannot read number of contacts for objects\n");
+    return 1;
+  }
+  for (size_t i = 0; i < n; ++i) {
+    if (fscanf(f, "%zu %zu ", &(p->contacts[i].o0), &(p->contacts[i].o1))!=2){
+      printf("Cannot read objects in contacts\n");
+      return 1;
+    }
+  }
+  check_word_in_file(f, "\n</DataArray>\n");
+  return 0;
+}
+
+static int particleProblemReadContactVtp(ParticleProblem *p, const char *dirname, int id) {
+  char *filename = malloc((strlen(dirname)+50)*sizeof(char));
+  sprintf(filename,"%s/contacts_%05d.vtp",dirname, id);
+  FILE *f = fopen(filename, "r");
+  if (!f){
+    printf("cannot open file %s\n", filename);
+    return -1;
+  }
+  check_word_in_file(f, "<VTKFile type=\"PolyData\">\n");
+  check_word_in_file(f, "<PolyData>\n");
+  check_word_in_file(f, "<FieldData>\n");
+  _readContactVtkSingle(p,f,PARTICLE_PARTICLE);
+  _readContactVtkSingle(p,f,PARTICLE_DISK);
+  _readContactVtkSingle(p,f,PARTICLE_SEGMENT);
+#if DIMENSION == 3
+  _readContactVtkSingle(p,f,PARTICLE_TRIANGLE);
+#endif
+  fclose(f);
+  free(filename);
+  return 0;
+}
+
+
+int particleProblemReadVtk(ParticleProblem *p, const char *dirname, int id) {
+  if (particleProblemReadParticleVtp(p,dirname,id) != 0)
+    return -1;
+  if (particleProblemReadBndVtu(p,dirname,id) != 0)
+    return -1;
+  if (particleProblemReadContactVtp(p,dirname,id) != 0)
+    return -1;
+
+  return 0;
+}
+
 void particleProblemLoad(ParticleProblem *p, const char *filename)
 {
   FILE *input = fopen(filename, "r");
@@ -1249,6 +1487,7 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
       exit(0);
     }
   }
+  fclose(input);
 }
 
 

marblesbag/src/fluid_problem.c

@@ -98,6 +98,18 @@ static double invDD(const double m[3][3], double inv[3][3]) {
 
 #endif
 
+static int check_word_in_file(FILE *f, const char *w) {
+  char word[256];
+  if(fscanf(f,"%255s", word) != 1){
+    printf("Invalid mesh file (\"%s\" expected).\n", w);
+    return -1;
+  }
+  if(strcmp(word, w) != 0){
+    printf("Invalid mesh file (\"%s\" expected).\n", w);
+    return -1;
+  }
+  return 0;
+}
 
 int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
 {
@@ -606,7 +618,7 @@ static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nB
   }
 }
 
-void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
+int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
 {
   static double totaltime = 0;
   struct timespec tic, toc;
@@ -618,7 +630,8 @@ void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
     solution_old[i] = problem->solution[i];
   double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
   double norm0 = 0;
-  for (int i=0; i<newton_max_it; ++i) {
+  int i;
+  for (i=0; i<newton_max_it; ++i) {
     double norm = 0;
     hxtLinearSystemZeroMatrix(problem->linear_system);
     for (int i=0; i<mesh->n_nodes*n_fields; ++i)
@@ -644,6 +657,7 @@ void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
   free(dx);
   free(rhs);
   free(solution_old);
+  return i;
 }
 
 static void fluid_problem_compute_node_volume(FluidProblem *problem) {

marblesbag/src/fluid_problem.h

@@ -39,7 +39,7 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
 int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir, double t, int iter);
 // complete force on the particle (including gravity)
 void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
-void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
+int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
 void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid);
 void fluid_problem_free(FluidProblem *problem);
 FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);