Changement format

fluid/CMakeLists.txt

@@ -21,6 +21,8 @@
 
 
 set(SRC
+  file_reader.c
+  gmsh_mesh.c
   fluid_problem.c
   fluid_problem_concentration.c
   mesh.c

fluid/file_reader.c

+#include "file_reader.h"
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+#include <errno.h>
+
+struct FileReaderStruct {
+  FILE *f;
+  char *file_name;
+  char line[4096];
+  char *pos;
+  char *prev;
+  long int iline;
+};
+
+void file_reader_delete(FileReader *r) {
+  fclose(r->f);
+  free(r->file_name);
+  free(r);
+}
+
+void file_reader_error(FileReader *r, const char *msg, ...) {
+  printf("Error in file  '%s' : ", r->file_name);
+  va_list args;
+  va_start(args, msg);
+  vprintf(msg,args);
+  printf("\n");
+  file_reader_delete(r);
+  exit(1);
+}
+
+void file_reader_error_at(FileReader *r, const char *msg,...) {
+  printf("Error in file '%s' line %lu column %lu : ", r->file_name, r->iline+1, r->prev-r->line+1);
+  va_list args;
+  va_start(args, msg);
+  vprintf(msg,args);
+  printf("\n");
+  printf("%s", r->line);
+  while(isblank(*r->prev))r->prev++;
+  for (int i = 0; i < (r->prev-r->line);++i) printf(" ");
+  printf("^\n");
+  file_reader_delete(r);
+  exit(1);
+}
+
+int file_reader_next_line_or_eof(FileReader *r) {
+  r->prev = r->pos;
+  if (r->pos){
+    while(isblank(*r->pos))r->pos++;
+    if (*r->pos != 0 && *r->pos != '\n')
+      file_reader_error_at(r,"end of line expected");
+  }
+
+  r->iline += 1;
+  char *l = fgets(r->line, 4096, r->f);
+  if (l == NULL) {
+    return 0;
+  }
+  r->pos = r->line;
+  return 1;
+}
+
+void file_reader_next_line_skip(FileReader *r) {
+  r->prev = r->pos;
+  r->iline += 1;
+  char *l = fgets(r->line, 4096, r->f);
+  if (l == NULL) {
+    file_reader_error(r, "Unexpected end of file.");
+  }
+  r->pos = r->line;
+}
+
+void file_reader_next_line(FileReader *r) {
+  r->prev = r->pos;
+  if (r->pos){
+    while(isblank(*r->pos))r->pos++;
+    if (*r->pos != 0 && *r->pos != '\n')
+      file_reader_error_at(r,"end of line expected");
+  }
+
+  r->iline += 1;
+  char *l = fgets(r->line, 4096, r->f);
+  if (l == NULL) {
+    file_reader_error(r, "Unexpected end of file.");
+  }
+  r->pos = r->line;
+}
+
+double file_reader_get_double(FileReader *r) {
+  r->prev = r->pos;
+  errno = 0;
+  double d = strtod(r->pos,&r->pos);
+  if (errno != 0)
+    file_reader_error_at(r,"expecting a double precision number");
+  return d;
+}
+
+long file_reader_get_int(FileReader *r) {
+  r->prev = r->pos;
+  errno = 0;
+  long d = strtol(r->pos,&r->pos,0);
+  if (errno != 0)
+    file_reader_error_at(r,"expecting an integer");
+  return d;
+}
+
+unsigned long file_reader_get_unsigned_int(FileReader *r) {
+  r->prev = r->pos;
+  errno = 0;
+  unsigned long d = strtoul(r->pos,&r->pos,0);
+  if (errno != 0)
+    file_reader_error_at(r,"expecting an unsigned integer");
+  return d;
+}
+
+void file_reader_get_word(FileReader *r, char *txt, int maxl) {
+  r->prev = r->pos;
+  char *next = r->pos;
+  while (isblank(*next)) next += 1;
+  int p = 0;
+  if (isspace(*next))
+    file_reader_error_at(r,"expecting a string");
+  while (!isspace(*next)) {
+    if (*next == 0)break;
+    txt[p++] = *next;
+    if(p+1 >= maxl) file_reader_error_at(r,"string is longer than expected");
+    next++;
+  }
+  txt[p] = '\0';
+  r->pos = next;
+}
+
+
+void file_reader_get_quoted_string(FileReader *r, char *txt, int maxl) {
+  r->prev = r->pos;
+  char *next = r->pos;
+  while (isblank(*next)) next += 1;
+  if (*next != '"') file_reader_error_at(r,"expecting a quoted string");
+  int p = 0;
+  next++;
+  while (*next != '"') {
+    if (*next == '\n' || *next == 0) file_reader_error_at(r,"expecting a quoted string");
+    txt[p++] = *next;
+    if(p+1 >= maxl) file_reader_error_at(r,"quoted string is longer than expected");
+    next++;
+  }
+  next++;
+  txt[p] = '\0';
+  r->pos = next;
+}
+
+void file_reader_assert(FileReader *r, const char *txt) {
+  r->prev = r->pos;
+  size_t l = strlen(txt);
+  if (strncmp(txt, r->pos, l) == 0) {
+    r->pos += l;
+    return;
+  }
+  file_reader_error_at(r,"expecting '%s'",txt);
+}
+
+FileReader *file_reader_new(const char *file_name) {
+  FileReader *r = malloc(sizeof(FileReader));
+  r->file_name = strdup(file_name);
+  r->f = fopen(file_name,"r");
+  if (!r->f) file_reader_error(r,"Cannot open file.");
+  r->iline = -1;
+  r->pos = NULL;
+  file_reader_next_line(r);
+  return r;
+}

fluid/file_reader.h

+#ifndef FILE_READER_H
+#define FILE_READER_H
+typedef struct FileReaderStruct FileReader;
+void file_reader_delete(FileReader *r);
+void file_reader_error(FileReader *r, const char *msg, ...);
+void file_reader_error_at(FileReader *r, const char *msg,...);
+int file_reader_next_line_or_eof(FileReader *r);
+void file_reader_next_line(FileReader *r);
+void file_reader_next_line_skip(FileReader *r);
+double file_reader_get_double(FileReader *r);
+long file_reader_get_int(FileReader *r);
+unsigned long file_reader_get_unsigned_int(FileReader *r);
+void file_reader_get_word(FileReader *r, char *txt, int maxl);
+void file_reader_get_quoted_string(FileReader *r, char *txt, int maxl);
+void file_reader_assert(FileReader *r, const char *txt);
+FileReader *file_reader_new(const char *file_name);
+#endif

fluid/gmsh_mesh.h

+#ifndef SLIM_GMSH_MESH_H
+#define SLIM_GMSH_MESH_H
+#include <stdlib.h>
+
+typedef struct {
+  int tag;
+  int dim;
+  int parent_tag;
+  int parent_dim;
+  size_t n_physicals;
+  int *physicals;
+  size_t n_partitions;
+  int *partitions;
+  size_t n_nodes;
+  int element_type;
+  double *x;
+  size_t *nodes_tag;
+  size_t n_elements;
+  size_t *elements; // lines in 1d, triangles in 2d
+  size_t *elements_tag;
+  size_t n_nodes_by_element;
+  size_t n_bounds;
+  double bbmin[3],bbmax[3];
+  int *bounds;
+} Entity;
+
+typedef struct {
+  int n_affine;
+  double *affine;
+  int dim;
+  int tag;
+  int parent_tag;
+  size_t n_nodes;
+  size_t *nodes; // 2*n_nodes
+}Periodic;
+
+typedef struct {
+  int n_partitions;
+  size_t n_entities[4];
+  Entity **entities[4];
+  int n_physical_names;
+  char **physical_names;
+  int *physical_name_tag;
+  int *physical_name_dim;
+  int n_periodic;
+  Periodic **periodic;
+} GmshMesh;
+
+GmshMesh *gmsh_mesh_new();
+void gmsh_mesh_read(GmshMesh *mesh, const char *file_name);
+void gmsh_mesh_write(const GmshMesh *mesh, const char *filename);
+void gmsh_mesh_free(GmshMesh *mesh);
+void gmsh_mesh_partition_2d(GmshMesh *mesh, int npart);
+Entity *gmsh_mesh_new_entity(GmshMesh *mesh, int dim, int tag,int parent_dim,int parent_tag);
+Entity *gmsh_mesh_get_entity(const GmshMesh *mesh, int dim, int tag);
+void gmsh_periodic_free(Periodic *p);
+#endif

fluid/mesh.c

@@ -2,26 +2,27 @@
  * MigFlow - Copyright (C) <2010-2018>
  * <Universite catholique de Louvain (UCL), Belgium
  *  Universite de Montpellier, France>
- * 	
+ *
  * List of the contributors to the development of MigFlow: see AUTHORS file.
  * Description and complete License: see LICENSE file.
- * 	
- * This program (MigFlow) is free software: 
- * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+ *
+ * This program (MigFlow) is free software:
+ * you can redistribute it and/or modify it under the terms of the GNU Lesser General
  * Public License as published by the Free Software Foundation, either version
  * 3 of the License, or (at your option) any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU Lesser General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU Lesser General Public License
- * along with this program (see COPYING and COPYING.LESSER files).  If not, 
+ * along with this program (see COPYING and COPYING.LESSER files).  If not,
  * see <http://www.gnu.org/licenses/>.
  */
 
 #include "mesh.h"
+#include "gmsh_mesh.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -184,22 +185,98 @@ static void sortDi(int i[DIMENSION]) {
 #else
    int n[3] = {i[0],i[1],i[2]};
    if (n[0]>n[1]) {
-      i[1]=n[0];    
-      i[0]=n[1];   
+      i[1]=n[0];
+      i[0]=n[1];
    } else {
-      i[1]=n[1];  
-      i[0]=n[0];  
-   } 
-   if (i[1]>n[2]) { 
-      i[2]=i[1];    
-      if(i[0]>n[2]){         
-         i[1]=i[0];                
+      i[1]=n[1];
+      i[0]=n[0];
+   }
+   if (i[1]>n[2]) {
+      i[2]=i[1];
+      if(i[0]>n[2]){
+         i[1]=i[0];
          i[0]=n[2];
-      } else i[1]=n[2];      
-   } else i[2]=n[2]; 
+      } else i[1]=n[2];
+   } else i[2]=n[2];
 #endif
 }
-
+Mesh *mesh_load_msh(const char *filename)
+{
+  GmshMesh* gmsh_m = gmsh_mesh_new();
+  printf("Try to read !\n");
+  gmsh_mesh_read(gmsh_m, filename);
+  printf("Mesh read !\n");
+  Mesh *m = malloc(sizeof(Mesh));
+  m->interfaces = NULL;
+  // Boundary
+  printf("Add Boundary ?\t");
+  m->n_boundaries = gmsh_m->n_physical_names;
+  m->boundary_names = malloc(sizeof(char*)*m->n_boundaries);
+  for(int i = 0; i < m->n_boundaries; ++i){
+    m->boundary_names[i] = strdup(gmsh_m->physical_names[i]);
+  }
+  printf("Boundary stored\n");
+  // Nodes
+  printf("Add Nodes ?\t");
+  m->n_nodes = 0;
+  for (int d = 0; d < 4; d ++){
+    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
+      const Entity *e = gmsh_m->entities[d][ie];
+      m->n_nodes += e->n_nodes;
+    }
+  }
+  printf("n_nodes : %d\n", m->n_nodes);
+  m->x = malloc(sizeof(double)*3*m->n_nodes);
+  int i_nodes = 0;
+  for (int d = 0; d < 4; ++d){
+    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
+      const Entity *e = gmsh_m->entities[d][ie];
+      for(size_t i = 0; i < e->n_nodes; ++i){
+        // e->x[i*3+0] e->x[i*3+1] e->x[i*3+2]
+        //printf("Entity[%d][%d]:\t Nodes %d :\t (%f %f %f)\n",d,ie,i,e->x[i*3+0],e->x[i*3+1],e->x[i*3+2]);
+        m->x[i_nodes+0] = e->x[i*3+0];
+        m->x[i_nodes+1] = e->x[i*3+1];
+        m->x[i_nodes+2] = e->x[i*3+2];
+        i_nodes += 3;
+      }
+    }
+  }
+  for(int i = 0; i < m->n_nodes; ++i){
+    printf("Nodes %d :\t (%f %f %f)\n",i,m->x[i*3+0],m->x[i*3+1],m->x[i*3+2]);
+  }
+  printf("Nodes stored\n");
+  // elements
+  printf("Add Elements ?\t");
+  m->n_elements = 0;
+  int n_nodes_by_element;
+  for (int d = 0; d < 4; ++d){
+    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
+      const Entity *e = gmsh_m->entities[d][ie];
+      m->n_elements += e->n_elements;
+      n_nodes_by_element = e->n_nodes_by_element;
+    }
+  }
+  printf("n_elements : %d\n", m->n_elements);
+  m->elements = malloc(sizeof(int)*m->n_elements*n_nodes_by_element);
+  int i_element = 0;
+  for (int d = 0; d < 4; d ++){
+    for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
+      const Entity *e = gmsh_m->entities[d][ie];
+      for(size_t i = 0; i < e->n_elements; ++i){
+        for(size_t j = 0; j < e->n_nodes_by_element; ++j){
+          m->elements[i_element] = e->elements[i*e->n_nodes_by_element+j]; // FIXME: the gmsh_free does not work with this line
+          printf("i_element = %d\n", i_element);
+          i_element++;
+        }
+      }
+    }
+  }
+  printf("Elements Stored\n");
+  gmsh_mesh_free(gmsh_m);
+  printf("Free GMSH Mesh Structure\n");
+  return m;
+}
+/*
 Mesh *mesh_load_msh(const char *filename)
 {
   Mesh *m = malloc(sizeof(Mesh));
@@ -359,6 +436,7 @@ Mesh *mesh_load_msh(const char *filename)
   }
   m->elements = realloc(m->elements, (DIMENSION+1)*sizeof(int)*m->n_elements);
   check_word_in_file(f,"$EndElements");
+
   fclose(f);
   free(physid);
   mesh_gen_edges(m,n_bnd_elements,bnd_elements,bnd_tags);
@@ -366,7 +444,7 @@ Mesh *mesh_load_msh(const char *filename)
   free(bnd_tags);
   return m;
 }
-
+*/
 int mesh_write_msh(Mesh *mesh, FILE *f)
 {
   fprintf(f,"$MeshFormat\n2.2 0 0\n$EndMeshFormat\n");

fluid/mesh.h

@@ -2,22 +2,22 @@
  * MigFlow - Copyright (C) <2010-2018>
  * <Universite catholique de Louvain (UCL), Belgium
  *  Universite de Montpellier, France>
- * 	
+ *
  * List of the contributors to the development of MigFlow: see AUTHORS file.
  * Description and complete License: see LICENSE file.
- * 	
- * This program (MigFlow) is free software: 
- * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+ *
+ * This program (MigFlow) is free software:
+ * you can redistribute it and/or modify it under the terms of the GNU Lesser General
  * Public License as published by the Free Software Foundation, either version
  * 3 of the License, or (at your option) any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU Lesser General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU Lesser General Public License
- * along with this program (see COPYING and COPYING.LESSER files).  If not, 
+ * along with this program (see COPYING and COPYING.LESSER files).  If not,
  * see <http://www.gnu.org/licenses/>.
  */
 
@@ -50,7 +50,7 @@ int mesh_write_pos_vector(const Mesh *mesh, const char *dir_name, const char *fi
 void mesh_set_elements(Mesh *m, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals);
 
 // not really necessary but useful to implement bnd conditions,
-// can be rebuilt from the interface list 
+// can be rebuilt from the interface list
 typedef struct {
   int n_interfaces;
   int *interfaces;
@@ -77,4 +77,3 @@ static int elbnd[24][3] = {
 };
 #endif
 #endif
-

testcases/couette-2d/mesh.geo

 rin = 0.0064;
 rout = 0.0254;
-lcin = rout/20; 
-lcout = rout/20;
+lcin = rout/10;
+lcout = rout/10;
 
 Point(1) = {0, 0, 0};
 
@@ -37,4 +37,4 @@ Physical Line("Outer") = {20, 21, 22, 23};
 Plane Surface(1) = {20, 10};
 Physical Surface("Domain") = {1};
 Physical Point("PtFix") = {13};
-Mesh.MshFileVersion = 2;
+Mesh.MshFileVersion = 4.1;

testcases/depot-2d/depot.py

 # MigFlow - Copyright (C) <2010-2018>
 # <Universite catholique de Louvain (UCL), Belgium
 #  Universite de Montpellier, France>
-# 	
+#
 # List of the contributors to the development of MigFlow: see AUTHORS file.
 # Description and complete License: see LICENSE file.
-# 	
-# This program (MigFlow) is free software: 
-# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+#
+# This program (MigFlow) is free software:
+# you can redistribute it and/or modify it under the terms of the GNU Lesser General
 # Public License as published by the Free Software Foundation, either version
 # 3 of the License, or (at your option) any later version.
-# 
+#
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Lesser General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU Lesser General Public License
-# along with this program (see COPYING and COPYING.LESSER files).  If not, 
+# along with this program (see COPYING and COPYING.LESSER files).  If not,
 # see <http://www.gnu.org/licenses/>.
 
 #!/usr/bin/env python
@@ -37,14 +37,14 @@ import random
 
 def genInitialPosition(filename, r, H, ly, lx, rhop) :
     """Set all the particles centre positions and create the particles objects to add in the computing structure
-    
-    Keyword arguments:    
+
+    Keyword arguments:
     filename -- name of the output file
     r -- max radius of the particles
     H -- domain height
     ly - particles area height
     lx -- particles area width
-    rhop -- particles density        
+    rhop -- particles density
     """
     # Particles structure builder
     p = scontact.ParticleProblem(2)

testcases/periodicTest/mesh.geo

+L = 5;
+H = 1;
+y = 0;
+lc = 1;//.05;
+
+Point(1) = {-L, H, 0, lc};
+Point(2) = {-L, H/2., 0, lc};
+Point(3) = {-L,0,0,lc};
+Point(4) = {L, 0, 0, lc};
+Point(5) = {L,H/2.,0,lc};
+Point(6) = {L, H, 0, lc};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 5};
+Line(5) = {5, 6};
+Line(6) = {6, 1};
+Line Loop(1) = {1:6};
+
+Plane Surface(1) = {1};
+
+Physical Line("LeftUp") = {1};
+Physical Line("LeftDown") = {2};
+Physical Line("RightDown") = {4};
+Physical Line("RightUp") = {5};
+Physical Line("Bottom") = {3};
+Physical Line("Top") = {6};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 4.1;

testcases/periodicTest/poiseuille.py

+# MigFlow - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+#
+# List of the contributors to the development of MigFlow: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+#
+# This program (MigFlow) is free software:
+# you can redistribute it and/or modify it under the terms of the GNU Lesser General
+# Public License as published by the Free Software Foundation, either version
+# 3 of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program (see COPYING and COPYING.LESSER files).  If not,
+# see <http://www.gnu.org/licenses/>.
+
+#!/usr/bin/env python
+from migflow import fluid as mbfluid
+from migflow import time_integration
+import numpy as np
+import os
+import subprocess
+import time
+import shutil
+import random
+import unittest
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+os.chdir(dir_path)
+# Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","2"])
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  np.array([0,0])                                      # gravity
+rho = 1000                                      # fluid density
+nu = 1e-3                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 100000                                     # final time
+
+#numerical parameters