Lecture du maillage dans scontact + update CMakeLists

fluid/gmsh_mesh.h

-#ifndef SLIM_GMSH_MESH_H
-#define SLIM_GMSH_MESH_H
+#ifndef GMSH_MESH_H
+#define GMSH_MESH_H
 #include <stdlib.h>
 
 typedef struct {

python/mesh.py

@@ -28,6 +28,7 @@ import sys
 import numpy as np
 
 dir_path = os.path.dirname(os.path.realpath(__file__))
+print(dir_path)
 lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
 lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
 

python/scontact.py

@@ -80,6 +80,12 @@ class ParticleState(np.ndarray) :
         state = cls((nparticles, 5 if dim == 2 else 9))
         return state
 
+class entity():
+    def __init__(self, dim, physicals, elements) :
+        self.dimension = dim
+        self.physicals = physicals
+        self.elements = elements
+
 
 class ParticleProblem :
     """Create the numerical structure containing all the physical particles that appear in the problem"""
@@ -503,54 +509,144 @@ class ParticleProblem :
         """
         if shift == None :
             shift = [0]*self._dim
-        m = mesh.Mesh(filename)
+        _mesh = self._lib.gmsh_mesh_load_msh(filename.encode())
+        n_phys = self.n_physical_names(_mesh)
+        phys_names = self.physical_names(_mesh, n_phys)
+        phys_tag = self.physical_tags(_mesh, n_phys)
+        phys_dims = self.physical_dims(_mesh,n_phys)
+        physicals = [{},{},{},{}]
+        for i in range(n_phys):
+            physicals[int(phys_dims[i])][phys_names[i].decode()] = int(phys_tag[i])
+        entities = self.get_entities(_mesh)
         self._addv = set()
         if self._dim == 2 :
-            pnum = [m.getPhysicalNumber(1, i) for i in tags]
-            for ent in m.entities :
+            pnum = [self.getPhysicalNumber(physicals,1,i) for i in tags]
+            for ent in entities :
                 if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
                     continue
                 tag = 0 if not ent.physicals else ent.physicals[0]
-                stag = m.getPhysicalName(ent.dimension, tag) or str(tag)
+                stag = self.getPhysicalName(physicals,ent.dimension, tag) or str(tag)
                 for i, el in enumerate(ent.elements) :
-                    for v in el.vertices :
+                    for v in el:
                         if v[3] in self._addv :
                             continue
                         self._addv.add(v[3])
                         self.add_boundary_disk((v[0] + shift[0], v[1] + shift[1]), 0.000, stag,material)
-                    v0 = el.vertices[1]
-                    v1 = el.vertices[0]
+                    v0 = el[1]
+                    v1 = el[0]
                     self.add_boundary_segment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), stag,material)
         else :
             self._adde = set()
-            pnum = [m.getPhysicalNumber(2, i) for i in tags]
-            for ent in m.entities :
+            pnum = [self.getPhysicalNumber(physicals,2, i) for i in tags]
+            for ent in entities :
                 if ent.dimension !=2 or not [i for i in pnum if i in ent.physicals]:
                     continue
                 tag = 0 if not ent.physicals else ent.physicals[0]
-                stag = m.getPhysicalName(ent.dimension, tag) or str(tag)
+                stag = self.getPhysicalName(physicals,ent.dimension, tag) or str(tag)
                 for i, el in enumerate(ent.elements) :
                     for j in range(3):
-                        v0 = el.vertices[j]
+                        v0 = el[j]
                         if not v0[3] in self._addv :
                             self._addv.add(v0[3])
                             self.add_boundary_disk((v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),0.,stag,material)
-                        v1 = el.vertices[(j+1)%3]
+                        v1 = el[(j+1)%3]
                         if (not (v0[3],v1[3]) in self._adde) and (not (v1[3],v0[3]) in self._adde) :
                             self._adde.add((v0[3],v1[3]))
                             self.add_boundary_segment(
                                 (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)
-
-                    v0 = el.vertices[1]
-                    v1 = el.vertices[0]
-                    v2 = el.vertices[2]
+                    v0 = el[1]
+                    v1 = el[0]
+                    v2 = el[2]
                     self.add_boundary_triangle(
                         (v0[0]+shift[0],v0[1]+shift[1],v0[2]+shift[2]),
                         (v1[0]+shift[0],v1[1]+shift[1],v1[2]+shift[2]),
                         (v2[0]+shift[0],v2[1]+shift[1],v2[2]+shift[2]),
                           stag,material)
-
     def dim(self) :
         return self._dim
+
+    def getPhysicalNumber(self,physicals, dim, name):
+        t = physicals[dim].get(name, None)
+        if not t :
+            physicals[dim][name] = max(physicals[dim].values()) + 1 if physicals[dim] else 1
+            t = physicals[dim][name]
+        return t
+    def getPhysicalName(self, physicals, dim, number):
+      for name, num in physicals[dim].items() :
+          if num == number :
+              return name
+      return ""
+    def get_entities(self,_mesh, dim=2):
+        entities = []
+        _nnodes = self.n_nodes(_mesh, dim)
+        _nodes = self.nodes(_mesh, dim, _nnodes)
+        _nodes_tag = self.nodes_tag(_mesh, dim, _nnodes)
+        vmap = {}
+        for i, node in enumerate(_nodes):
+            vmap[int(_nodes_tag[i])] = [*_nodes[i], _nodes_tag[i]]
+        for d in range(4):
+            for ie in range(self.n_entities(_mesh,d)):
+                edim = self.entity_dim(_mesh,d,ie)
+                ephys = self.entity_physicals(_mesh,d,ie)
+                nelt = self.n_elements(_mesh,d,ie)
+                elts = [None]*nelt
+                elts = [[vmap[int(i)] for i in self.element_nodes_tag(_mesh,d,ie,id)] for id in range(nelt)]
+                entities.append(entity(edim, ephys, elts))
+        return entities
+    def n_physical_names(self,_mesh):
+        self._lib.gmsh_mesh_n_physical_names.restype = c_int
+        return self._lib.gmsh_mesh_n_physical_names(_mesh)
+    def physical_names(self,_mesh,n_physical_names):
+        self._lib.gmsh_mesh_physical_name.restype = c_char_p
+        return [self._lib.gmsh_mesh_physical_name(_mesh, c_int(i)) for i in range(n_physical_names)]
+    def physical_tags(self,_mesh,n_physical_names):
+        f = getattr(self._lib,"gmsh_mesh_physical_name_tags")
+        f.restype = POINTER(c_int)
+        return np.ctypeslib.as_array(f(_mesh), shape=(n_physical_names,)).flatten()
+    def physical_dims(self, _mesh, n_physical_names):
+        f = getattr(self._lib,"gmsh_mesh_physical_name_dims")
+        f.restype = POINTER(c_int)
+        return np.ctypeslib.as_array(f(_mesh), shape=(n_physical_names,)).flatten()
+    def n_nodes(self,_mesh, dim):
+        self._lib.gmsh_mesh_n_nodes.restype = c_int
+        return self._lib.gmsh_mesh_n_nodes(_mesh, c_int(dim))
+    def nodes(self,_mesh, dim, nnodes):
+        f = getattr(self._lib, "gmsh_mesh_nodes")
+        f.restype = POINTER(c_double)
+        return np.ctypeslib.as_array(f(_mesh, c_int(dim)),shape=(nnodes,3))
+    def nodes_tag(self,_mesh, dim, nnodes):
+        f = getattr(self._lib, "gmsh_mesh_nodes_tag")
+        f.restype = POINTER(c_size_t)
+        return np.ctypeslib.as_array(f(_mesh, c_int(dim)),shape=(nnodes,))
+    def n_entities(self,_mesh, dim):
+        self._lib.gmsh_mesh_n_entities.restype = c_size_t
+        return self._lib.gmsh_mesh_n_entities(_mesh, c_int(dim))
+    def entity_dim(self,_mesh,dim,ie):
+        self._lib.gmsh_mesh_entity_dim.restype = c_int
+        return self._lib.gmsh_mesh_entity_dim(_mesh, c_int(dim), c_int(ie))
+    def n_physicals(self,_mesh,dim,ie):
+        self._lib.gmsh_mesh_entity_n_physicals.restype = c_size_t
+        return self._lib.gmsh_mesh_entity_n_physicals(_mesh,c_int(dim),c_int(ie))
+    def entity_n_physicals(self,_mesh,dim,ie):
+        self._lib.gmsh_mesh_entity_n_physicals.restype = c_size_t
+        return self._lib.gmsh_mesh_entity_n_physicals(_mesh, c_int(dim), c_int(ie))
+    def entity_physicals(self,_mesh,dim,ie):
+        n_physicals = self.entity_n_physicals(_mesh,dim,ie)
+        if n_physicals == 0 :
+            return np.array([])
+        f = getattr(self._lib, "gmsh_mesh_entity_physicals")
+        f.restype = POINTER(c_int)
+        return np.ctypeslib.as_array(f(_mesh, c_int(dim), c_int(ie)),shape=(n_physicals,))
+    def n_elements(self,_mesh,dim,ie):
+        self._lib.gmsh_mesh_entity_n_elements.restype = c_size_t
+        return self._lib.gmsh_mesh_entity_n_elements(_mesh, c_int(dim), c_int(ie))
+    def element_n_nodes(self,_mesh,dim,ie):
+        self._lib.gmsh_mesh_element_n_nodes.restype = c_size_t
+        return self._lib.gmsh_mesh_element_n_nodes(_mesh, c_int(dim), c_int(ie))
+    def element_nodes_tag(self,_mesh,dim,ie, ielt):
+        nodes_by_elt = self.element_n_nodes(_mesh,dim,ie)
+        f = getattr(self._lib, "gmsh_mesh_element_nodes_tag")
+        f.restype = POINTER(c_size_t)
+        return None if nodes_by_elt < 1 else np.ctypeslib.as_array(f(_mesh, c_int(dim), c_int(ie), c_int(ielt)),shape=(nodes_by_elt,))

scontact/CMakeLists.txt

@@ -22,6 +22,8 @@
 set(SRC
   scontact.c
   ../tools/quadtree.c
+  ../fluid/file_reader.c
+  ../fluid/gmsh_mesh.c
 )
 
 include_directories(. ../tools)

testcases/periodicTest/poiseuille.py

@@ -84,7 +84,8 @@ fluid.export_vtk(outputdir,0,0)
 
 ii = 0
 tic = time.time()
-while t < tEnd :
+# while t < tEnd :
+while ii < 10:
     #Fluid solver
     time_integration.iterate(fluid,None,dt)
     t += dt

testcases/periodicTest/poiseuilleGrain.py

@@ -104,6 +104,7 @@ G = p.mass()*g
 ii = 0
 tic = time.time()
 while t < tEnd :
+# while ii < 10:
     #Fluid solver
     time_integration.iterate(fluid,p,dt,external_particles_forces=G)
     t += dt