Lecture Mesh via gmsh_mesh pour scontact

fluid/gmsh_mesh.c

@@ -750,3 +750,174 @@ void gmsh_mesh_write_partition(const GmshMesh *mesh, const char *filename, int p
 void gmsh_mesh_write(const GmshMesh *mesh, const char *filename) {
   gmsh_mesh_write_partition(mesh, filename, -1);
 }
+
+
+
+
+
+static int intcmp(const void *p0, const void *p1) {
+  return *(int*)p0 - *(int*)p1;
+}
+static int get_node_id(int id, int *nodes_map, int n_nodes) {
+  int *found = bsearch(&id, nodes_map, n_nodes, sizeof(int), intcmp);
+  return found ? found-nodes_map : -1;
+}
+
+GmshMesh *gmsh_mesh_load_msh(const char *filename){
+  GmshMesh *mesh = gmsh_mesh_new();
+  gmsh_mesh_read(mesh,filename);
+  return mesh;
+}
+
+int gmsh_mesh_n_physical_names(GmshMesh *mesh) {
+  return mesh->n_physical_names;
+}
+int* gmsh_mesh_physical_name_tags(GmshMesh *mesh) {
+  return mesh->physical_name_tag;
+}
+int* gmsh_mesh_physical_name_dims(GmshMesh *mesh) {
+  return mesh->physical_name_dim;
+}
+char* gmsh_mesh_physical_name(GmshMesh *mesh, int i) {
+  return mesh->physical_names[i];
+}
+size_t gmsh_mesh_n_elements(GmshMesh* mesh, const int dim) {
+  int n_elements = 0;
+  for(int ie = 0; ie < mesh->n_entities[dim]; ++ie){
+    const Entity *e = mesh->entities[dim][ie];
+    if (e->n_nodes_by_element != dim +1) {
+      exit(1);
+    }
+    n_elements += e->n_elements;
+  }
+  return n_elements;
+}
+int* gmsh_mesh_elements(GmshMesh* mesh, const int dim) {
+  int n_elements = gmsh_mesh_n_elements(mesh, dim);
+  int* elements = malloc(sizeof(int)*n_elements*(dim+1));
+  int i_elt = 0;
+  for(int ie = 0; ie < mesh->n_entities[dim]; ++ie){
+    const Entity *e = mesh->entities[dim][ie];
+    for(size_t i = 0; i < e->n_elements; ++i){
+      for(size_t j = 0; j < e->n_nodes_by_element; ++j){
+        elements[i_elt++] = e->elements[i*e->n_nodes_by_element+j];
+      }
+    }
+  }
+  return elements;
+}
+int* gmsh_mesh_nodes_map(GmshMesh *mesh, const int dim) {
+  int n_elements = gmsh_mesh_n_elements(mesh, dim);
+  int* elements = gmsh_mesh_elements(mesh, dim);
+  int n_enodes = n_elements*(dim+1);
+  int *nodes_map = malloc(sizeof(int)*n_enodes);
+  memcpy(nodes_map, elements, sizeof(int)*n_enodes);
+  qsort(nodes_map, n_enodes, sizeof(int), intcmp);
+  int n_nodes = 0;
+  for(int i = 0; i < n_enodes; ++i){
+    if(i == 0 || nodes_map[i] != nodes_map[n_nodes-1]){
+      nodes_map[n_nodes++] = nodes_map[i];
+    }
+  }
+  free(elements);
+  return nodes_map;
+}
+int gmsh_mesh_n_nodes(GmshMesh *mesh, const int dim) {
+  int n_elements = gmsh_mesh_n_elements(mesh, dim);
+  int* elements = gmsh_mesh_elements(mesh, dim);
+  int n_enodes = n_elements*(dim+1);
+  int *nodes_map = malloc(sizeof(int)*n_enodes);
+  memcpy(nodes_map, elements, sizeof(int)*n_enodes);
+  qsort(nodes_map, n_enodes, sizeof(int), intcmp);
+  int n_nodes = 0;
+  for(int i = 0; i < n_enodes; ++i){
+    if(i == 0 || nodes_map[i] != nodes_map[n_nodes-1]){
+      nodes_map[n_nodes++] = nodes_map[i];
+    }
+  }
+  free(elements);
+  free(nodes_map);
+  return n_nodes;
+}
+double* gmsh_mesh_nodes(GmshMesh *mesh, const int dim) {
+  int n_nodes = gmsh_mesh_n_nodes(mesh, dim);
+  int* nodes_map  = gmsh_mesh_nodes_map(mesh, dim);
+
+  double* x = (double*) malloc(sizeof(double)*3*n_nodes);
+  for (int d = 0; d < 4; ++d){
+    for(int ie = 0; ie < mesh->n_entities[d]; ++ie){
+      const Entity *e = mesh->entities[d][ie];
+      for(size_t i = 0; i < e->n_nodes; ++i){
+        int nid = get_node_id(e->nodes_tag[i], nodes_map, n_nodes);
+        if (nid >= 0) {
+          for (int id = 0; id < 3; ++id) {
+            x[nid*3+id] = e->x[i*3+id];
+          }
+        }
+      }
+    }
+  }
+
+  free(nodes_map);
+  return x;
+}
+size_t* gmsh_mesh_nodes_tag(GmshMesh *mesh, const int dim) {
+  int n_nodes = gmsh_mesh_n_nodes(mesh, dim);
+  int* nodes_map  = gmsh_mesh_nodes_map(mesh, dim);
+
+  size_t* tag = (size_t*) malloc(sizeof(size_t)*n_nodes);
+  for (int d = 0; d < 4; ++d){
+    for(int ie = 0; ie < mesh->n_entities[d]; ++ie){
+      const Entity *e = mesh->entities[d][ie];
+      for(size_t i = 0; i < e->n_nodes; ++i){
+        int nid = get_node_id(e->nodes_tag[i], nodes_map, n_nodes);
+        if (nid >= 0) {
+          tag[nid] = e->nodes_tag[i];
+        }
+      }
+    }
+  }
+  free(nodes_map);
+  return tag;
+}
+size_t gmsh_mesh_n_entities(GmshMesh *mesh, int dim){
+  return mesh->n_entities[dim];
+}
+int gmsh_mesh_entity_tag(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->tag;
+}
+int gmsh_mesh_entity_dim(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->dim;
+}
+size_t gmsh_mesh_entity_n_physicals(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->n_physicals;
+}
+int* gmsh_mesh_entity_physicals(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->physicals;
+}
+size_t gmsh_mesh_entity_n_elements(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->n_elements;
+}
+int gmsh_mesh_element_type(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->element_type;
+}
+size_t gmsh_mesh_element_tag(GmshMesh *mesh, int dim, int ie, int id) {
+  return mesh->entities[dim][ie]->elements_tag[id];
+}
+size_t gmsh_mesh_element_n_partitions(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->n_partitions;
+}
+int* gmsh_mesh_element_partitions(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->partitions;
+}
+size_t gmsh_mesh_element_n_nodes(GmshMesh *mesh, int dim, int ie) {
+  return mesh->entities[dim][ie]->n_nodes_by_element;
+}
+size_t* gmsh_mesh_element_nodes_tag(GmshMesh *mesh, int dim, int ie, int ielt) {
+  const Entity* e = mesh->entities[dim][ie];
+  size_t* element = (size_t*) malloc(sizeof(size_t)*e->n_nodes_by_element);
+  for (int j = 0; j < e->n_nodes_by_element; ++j) {
+    element[j] = e->elements[ielt*e->n_nodes_by_element+j];
+  }
+  return element;
+}

fluid/gmsh_mesh.h

@@ -54,4 +54,27 @@ 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);
+
+GmshMesh *gmsh_mesh_load_msh(const char *filename);
+int gmsh_mesh_n_physical_names(GmshMesh *mesh);
+char* gmsh_mesh_physical_name(GmshMesh *mesh, int i);
+int* gmsh_mesh_physical_name_tags(GmshMesh *mesh);
+int* gmsh_mesh_physical_name_dims(GmshMesh *mesh);
+int gmsh_mesh_n_nodes(GmshMesh *mesh, const int dimension);
+int* gmsh_mesh_nodes_map(GmshMesh *mesh, const int dimension);
+double* gmsh_mesh_nodes(GmshMesh *mesh, const int dimension);
+size_t* gmsh_mesh_nodes_tag(GmshMesh *mesh, const int dimension);
+
+size_t gmsh_mesh_n_entities(GmshMesh *mesh, int dimension);
+int gmsh_mesh_entity_tag(GmshMesh *mesh, int dimension, int ie);
+int gmsh_mesh_entity_dim(GmshMesh *mesh, int dimension, int ie);
+size_t gmsh_mesh_entity_n_physicals(GmshMesh *mesh, int dim, int ie);
+int* gmsh_mesh_entity_physicals(GmshMesh *mesh, int dimension, int ie);
+size_t gmsh_mesh_entity_n_elements(GmshMesh *mesh, int dimension, int ie);
+int gmsh_mesh_element_type(GmshMesh *mesh, int dimension, int ie);
+size_t gmsh_mesh_element_tag(GmshMesh *mesh, int dim, int ie, int id);
+size_t gmsh_mesh_element_n_partitions(GmshMesh *mesh, int dim, int ie);
+int* gmsh_mesh_element_partitions(GmshMesh *mesh, int dim, int ie);
+size_t gmsh_mesh_element_n_nodes(GmshMesh *mesh, int dim, int ie);
+size_t* gmsh_mesh_element_nodes_tag(GmshMesh *mesh, int dim, int ie, int ielt);
 #endif

python/mesh.py

@@ -20,7 +20,16 @@
 # see <http://www.gnu.org/licenses/>.
 
 from . import gmshType
+from ctypes import *
 import time
+import signal
+import os
+import sys
+import numpy as np
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
+lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
 
 class Mesh() :
   class element():
@@ -47,17 +56,24 @@ class Mesh() :
       self.physicals = physicals
       self.elements = []
 
-  def __init__(self, filename = None) :
+  def __init__(self, filename = None, dim = 2) :
     self.vertices = []
     self.physicals = [{}, {}, {}, {}]
     self.entities = []
     self.maxeid = 0
     self.useFormat3 = False
-    # self.useFormat4 = False
+    # A generaliser : Comment connaitre la dimension du probleme
+    self.dim = dim
+    self._mesh = None
+    if self.dim == 2 :
+        self._lib = lib2
+    elif self.dim == 3 :
+        self._lib = lib3
     if filename :
       self.read(filename)
 
   def write(self, filename, format3=None) :
+
     with open(filename, "w") as output :
         if format3 == None :
           format3 = self.useFormat3
@@ -104,67 +120,170 @@ class Mesh() :
           output.write("$EndElements\n")
 
   def read(self, filename):
-    with open(filename, "r") as fin :
-        l = fin.readline()
-        vmap = {}
-        entitymap = {}
-        while l != "" :
-          w = l.split()
-          if w[0] == "$MeshFormat":
-            l = fin.readline().split()
-            if float(l[0]) == 3.:
-              self.useFormat3 = True
-            elif int(float(l[0])) == 2 :
-              self.useFormat3 = False
-            else :
-              print("error : cannot read mesh format " + l[0])
-            l = fin.readline()
-          elif w[0] == "$PhysicalNames" :
-            n = int(fin.readline())
-            for i in range(n) :
-              dim, tag, name = fin.readline().split()
-              self.physicals[int(dim)][name[1:-1]] = int(tag)
-            fin.readline()
-          elif w[0] == "$Entities" and self.useFormat3:
-            n = int(fin.readline())
-            print(n)
-            for i in range(n) :
-              l = fin.readline().split()
-              j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
-              self.entities.append(Mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
-              entitymap[(dim, j)] = self.entities[-1]
-            fin.readline()
-          elif w[0] == "$Nodes" :
-            n = int(fin.readline())
-            for i in range(n) :
-              if self.useFormat3 :
-                (j, x, y, z, t) = fin.readline().split()
-              else :
-                (j, x, y, z) = fin.readline().split()
-              self.vertices.append([float(x), float(y), float(z), int(j)])
-              vmap[int(j)] = self.vertices[-1]
-          elif w[0] == "$Elements" :
-            n = int(fin.readline())
-            for i in range(n) :
-              l = fin.readline().split()
-              if self.useFormat3 :
-                j, t, e, nf = int(l[0]), int(l[1]), int(l[2]), int(l[3])
-                nv = gmshType.Type[t].numVertices
-                vertices = [vmap[int(i)] for i in l[4:4+nv]]
-                partition = [int(i) for i in l[5 + nv : 5 + nv + int(l[4 + nv])]] if nf > nv else []
-              else :
-                j, t, nf, p, e = int(l[0]), int(l[1]), int(l[2]), int(l[3]), int(l[4])
-                vertices = [vmap[int(i)] for i in l[3 + nf:]]
-                partition = [int(i) for i in l[6 : 6 + int(l[5])]] if nf > 2 else []
-              edim = gmshType.Type[t].baseType.dimension
-              entity = entitymap.get((edim, e), None)
-              if not self.useFormat3 and not entity:
-                entity = Mesh.entity(e, edim, [p])
-                self.entities.append(entity)
-                entitymap[(edim, e)] = entity
-              entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
-              self.maxeid = max(self.maxeid, j)
-          l = fin.readline()
+    self._mesh = self._lib.gmsh_mesh_load_msh(filename.encode())
+    phys_dims = self.physical_name_dims()
+    phys_names = self.physical_names()
+    phys_tag = self.physical_name_tags()
+    for i in range(self.n_physical_names()):
+        self.physicals[int(phys_dims[i])][phys_names[i].decode()] = int(phys_tag[i])
+    nodes = self.nodes()
+    nodes_tag = self.nodes_tag()
+    vmap = {}
+    self.vertices = [[0]*4]*self.n_nodes()
+    for i, node in enumerate(nodes):
+        self.vertices[i] = [*node, nodes_tag[i]]
+        vmap[int(nodes_tag[i])] = self.vertices[i]
+    self.entities = []
+    for d in range(4):
+        for ie in range(self.n_entities(d)):
+            tag = self.entity_tag(d,ie)
+            dim = self.entity_dim(d,ie)
+            etype = gmshType.Type[self.element_type(d,ie)]
+            phys = self.entity_physicals(d,ie) # return a RunTime Warning
+            e = Mesh.entity(tag,dim,phys)
+            n_elt = self.n_elements(d,ie)
+            e.elements = [None]*n_elt
+            for id in range(n_elt):
+                elt_tag = self.element_tag(d,ie,id)
+                partition = self.element_partitions(d,ie)
+                elt_nodes_tag = self.element_nodes_tag(d,ie,id)
+                vertices = [vmap[int(i)] for i in elt_nodes_tag]
+                e.elements[id] = Mesh.element(etype, vertices, partition, elt_tag)
+                self.maxeid = max(self.maxeid, elt_tag)
+            self.entities.append(e)
+
+    # with open(filename, "r") as fin :
+    #     l = fin.readline()
+    #     vmap = {}
+    #     entitymap = {}
+    #     while l != "" :
+    #       w = l.split()
+    #       if w[0] == "$MeshFormat":
+    #         l = fin.readline().split()
+    #         if float(l[0]) == 3.:
+    #           self.useFormat3 = True
+    #         elif int(float(l[0])) == 2 :
+    #           self.useFormat3 = False
+    #         else :
+    #           print("error : cannot read mesh format " + l[0])
+    #         l = fin.readline()
+    #       elif w[0] == "$PhysicalNames" :
+    #         n = int(fin.readline())
+    #         for i in range(n) :
+    #           dim, tag, name = fin.readline().split()
+    #           self.physicals[int(dim)][name[1:-1]] = int(tag)
+    #         fin.readline()
+    #       elif w[0] == "$Entities" and self.useFormat3:
+    #         n = int(fin.readline())
+    #         for i in range(n) :
+    #           l = fin.readline().split()
+    #           j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
+    #           self.entities.append(Mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
+    #           entitymap[(dim, j)] = self.entities[-1]
+    #         fin.readline()
+    #       elif w[0] == "$Nodes" :
+    #         n = int(fin.readline())
+    #         for i in range(n) :
+    #           if self.useFormat3 :
+    #             (j, x, y, z, t) = fin.readline().split()
+    #           else :
+    #             (j, x, y, z) = fin.readline().split()
+    #           self.vertices.append([float(x), float(y), float(z), int(j)])
+    #           vmap[int(j)] = self.vertices[-1]
+    #       elif w[0] == "$Elements" :
+    #         n = int(fin.readline())
+    #         for i in range(n) :
+    #           l = fin.readline().split()
+    #           if self.useFormat3 :
+    #             j, t, e, nf = int(l[0]), int(l[1]), int(l[2]), int(l[3])
+    #             nv = gmshType.Type[t].numVertices
+    #             vertices = [vmap[int(i)] for i in l[4:4+nv]]
+    #             partition = [int(i) for i in l[5 + nv : 5 + nv + int(l[4 + nv])]] if nf > nv else []
+    #           else :
+    #             j, t, nf, p, e = int(l[0]), int(l[1]), int(l[2]), int(l[3]), int(l[4]) #nf : number of tags
+    #             vertices = [vmap[int(i)] for i in l[3 + nf:]]
+    #             partition = [int(i) for i in l[6 : 6 + int(l[5])]] if nf > 2 else []
+    #           edim = gmshType.Type[t].baseType.dimension
+    #           entity = entitymap.get((edim, e), None)
+    #           if not self.useFormat3 and not entity:
+    #             entity = Mesh.entity(e, edim, [p])
+    #             self.entities.append(entity)
+    #             entitymap[(edim, e)] = entity
+    #           # print(vertices)
+    #           entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
+    #           self.maxeid = max(self.maxeid, j)
+    #       l = fin.readline()
+
+  def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
+    f = getattr(self._lib,"gmsh_mesh_"+f_name)
+    f.restype = POINTER(typec)
+    return np.ctypeslib.as_array(f(self._mesh),shape=(nrow,ncol))
+
+
+  def n_physical_names(self):
+      self._lib.gmsh_mesh_n_physical_names.restype = c_int
+      return self._lib.gmsh_mesh_n_physical_names(self._mesh)
+  def physical_names(self):
+      self._lib.gmsh_mesh_physical_name.restype = c_char_p
+      return [self._lib.gmsh_mesh_physical_name(self._mesh, c_int(i)) for i in range(self.n_physical_names())]
+  def physical_name_tags(self):
+      return self._get_matrix("physical_name_tags", self.n_physical_names(), 1, typec=c_int).flatten()
+  def physical_name_dims(self):
+      return self._get_matrix("physical_name_dims", self.n_physical_names(), 1, typec=c_int).flatten()
+
+  def n_nodes(self):
+      self._lib.gmsh_mesh_n_nodes.restype = c_int
+      return self._lib.gmsh_mesh_n_nodes(self._mesh, c_int(self.dim))
+  def nodes(self):
+      f = getattr(self._lib, "gmsh_mesh_nodes")
+      f.restype = POINTER(c_double)
+      return np.ctypeslib.as_array(f(self._mesh, c_int(self.dim)),shape=(self.n_nodes(),3))
+  def nodes_tag(self):
+      f = getattr(self._lib,"gmsh_mesh_nodes_tag")
+      f.restype = POINTER(c_size_t)
+      return np.ctypeslib.as_array(f(self._mesh, c_int(self.dim)), shape=(self.n_nodes(),))
+  def n_entities(self,dim):
+      self._lib.gmsh_mesh_n_entities.restype = c_size_t
+      return self._lib.gmsh_mesh_n_entities(self._mesh, c_int(dim))
+  def entity_tag(self,dim,ie):
+      self._lib.gmsh_mesh_entity_tag.restype = c_int
+      return self._lib.gmsh_mesh_entity_tag(self._mesh, c_int(dim), c_int(ie))
+  def entity_dim(self,dim,ie):
+      self._lib.gmsh_mesh_entity_dim.restype = c_int
+      return self._lib.gmsh_mesh_entity_dim(self._mesh, c_int(dim), c_int(ie))
+  def entity_n_physicals(self,dim,ie):
+      self._lib.gmsh_mesh_entity_n_physicals.restype = c_size_t
+      return self._lib.gmsh_mesh_entity_n_physicals(self._mesh, c_int(dim), c_int(ie))
+  def entity_physicals(self,dim,ie):
+      n_physicals = self.entity_n_physicals(dim,ie)
+      f = getattr(self._lib, "gmsh_mesh_entity_physicals")
+      f.restype = POINTER(c_int)
+      return np.ctypeslib.as_array(f(self._mesh, c_int(dim), c_int(ie)),shape=(n_physicals,))
+  def n_elements(self,dim,ie):
+      self._lib.gmsh_mesh_entity_n_elements.restype = c_size_t
+      return self._lib.gmsh_mesh_entity_n_elements(self._mesh, c_int(dim), c_int(ie))
+  def element_type(self,dim,ie):
+      self._lib.gmsh_mesh_element_type.restype = c_int
+      return self._lib.gmsh_mesh_element_type(self._mesh, c_int(dim), c_int(ie))
+  def element_tag(self,dim,ie,id):
+      self._lib.gmsh_mesh_element_tag.restype = c_size_t
+      return self._lib.gmsh_mesh_element_tag(self._mesh, c_int(dim), c_int(ie),c_int(id))
+  def element_n_partitions(self,dim,ie):
+      self._lib.gmsh_mesh_element_n_partitions.restype = c_size_t
+      return self._lib.gmsh_mesh_element_n_partitions(self._mesh, c_int(dim), c_int(ie))
+  def element_partitions(self,dim,ie):
+      n_partitions = self.element_n_partitions(dim,ie)
+      f = getattr(self._lib, "gmsh_mesh_element_partitions")
+      f.restype = POINTER(c_int)
+      return np.ctypeslib.as_array(f(self._mesh, c_int(dim), c_int(ie)),shape=(n_partitions,))
+  def element_n_nodes(self,dim,ie):
+      self._lib.gmsh_mesh_element_n_nodes.restype = c_size_t
+      return self._lib.gmsh_mesh_element_n_nodes(self._mesh, c_int(dim), c_int(ie))
+  def element_nodes_tag(self,dim,ie, ielt):
+      nodes_by_elt = self.element_n_nodes(dim,ie)
+      f = getattr(self._lib, "gmsh_mesh_element_nodes_tag")
+      f.restype = POINTER(c_size_t)
+      return np.ctypeslib.as_array(f(self._mesh, c_int(dim), c_int(ie), c_int(ielt)),shape=(nodes_by_elt,))
 
   def getPhysicalNumber(self, dim, name) :
     t = self.physicals[dim].get(name, None)

python/scontact.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/>.
 
 """Model for Immersed Granular Flow -- Particles user interface
@@ -24,7 +24,7 @@
     Contact: jonathan.lambrechts@uclouvain.be
     Webpage: www.migflow.be
 
-    MigFlow computes immersed granular flows using an unresolved FEM-DEM model. 
+    MigFlow computes immersed granular flows using an unresolved FEM-DEM model.
     The discrete phase is computed by updating iteratively the particle velocities until a set of velocities respecting the non-interpenetration contact law for the next time step is found
 """
 
@@ -103,7 +103,7 @@ class ParticleProblem :
         size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//8
         buf = (size*c_double).from_address(ptr)
         return np.ctypeslib.as_array(buf).reshape([-1,ncol])
-    
+
     def _get_idx(self, fName) :
         f = getattr(self._lib,"particleProblem"+fName)
         f.restype = c_void_p
@@ -132,7 +132,7 @@ class ParticleProblem :
               self._lib = lib2f if rotation_enabled else lib2fwr
             else:
               self._lib = lib2
-            
+
             self._coord_type = c_double*2
         elif dim == 3 :
             self._lib = (lib3f if rotation_enabled else lib3fwr) if (friction_enabled) else lib3
@@ -219,7 +219,7 @@ class ParticleProblem :
         if self._dim == 2:
             return np.ndarray((0),self._triangletype)
         return self._get_array("Triangle",self._triangletype)
-    
+
     def contact_forces(self):
         return (self._get_matrix("ContactForces",self._dim))
 
@@ -292,7 +292,7 @@ class ParticleProblem :
         dt -- numerical time step
         forces -- list of vectors containing the forces applied on the particles
         tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
-        """ 
+        """
         self._lib.particleProblemIterate.restype = c_int
         return self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1),c_int(force_motion),_np2c(forces))
 
@@ -362,7 +362,7 @@ class ParticleProblem :
         x = np.array([[0.,0.,0.]])
         VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
         VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
-        
+
     def read_vtk(self,dirname,i,contact_factor=1) :
         """Read an existing output file to set the particle data."""
         self._lib.particleProblemClearBoundaries(self._p)
@@ -381,7 +381,7 @@ class ParticleProblem :
         tagmap = {}
         for j,n in enumerate(tagnames) :
             tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
-        
+
         offsets = np.hstack([[0],cells["offsets"]])
         el = cells["connectivity"]
         tags = np.vectorize(tagmap.get)(cdata["Tag"])
@@ -472,7 +472,7 @@ class ParticleProblem :
       """Return the names of the boundaries """
       self._lib.particleProblemGetTagName.restype = c_char_p
       ntagname = self._lib.particleProblemNTag(self._p)
-      tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)]) 
+      tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
       return tagnames
 
     def clear_boundaries(self):
@@ -494,7 +494,7 @@ class ParticleProblem :
 
     def load_msh_boundaries(self, filename, tags, shift=None,material="default") :
         """Load the numerical domain and set the physical boundaries the particles cannot go through.
-        
+
         Keyword arguments:
         filename -- name of the msh file
         tags -- tags of physical boundary defined in the msh file
@@ -542,7 +542,7 @@ class ParticleProblem :
                                 (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
                                 (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
                                 stag,material)
-