porosity_periodic

fluid/fluid_problem.c

@@ -1212,32 +1212,49 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
   free(problem->node_volume);
   Mesh *mesh = problem->mesh;
   problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
-  for (int i = 0; i < problem->mesh->n_nodes; ++i){
-    problem->node_volume[i] = 0;
+  int independent_nodes = mesh->n_nodes - mesh->n_periodic;
+  double* periodic_node_volume = (double*)malloc(sizeof(double)*independent_nodes);
+  for (int i = 0; i < independent_nodes; ++i){
+    periodic_node_volume[i] = 0;
   }
+  int* periodic_mapping = mesh->periodic_mapping;
   for (int iel = 0; iel < mesh->n_elements; ++iel) {
     double dxidx[D][D];
     const int *el = &mesh->elements[iel*N_N];
     const double detj = mesh_dxidx(mesh, iel, dxidx);
     const double vol = node_volume_from_detj(detj);
     for (int i = 0; i < N_N; ++i)
-      problem->node_volume[el[i]] += vol;
+      // problem->node_volume[el[i]] += vol;
+      periodic_node_volume[periodic_mapping[el[i]]] += vol;
   }
+  for(int i = 0; i < problem->mesh->n_nodes; ++i){
+    problem->node_volume[i] = periodic_node_volume[periodic_mapping[i]];
+  }
+  free(periodic_node_volume);
 }
 
 static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw)
 {
-  for (int i = 0; i < mesh->n_nodes; ++i){
-    porosity[i] = 0;
+  int independent_nodes = mesh->n_nodes - mesh->n_periodic;
+  double* periodic_porosity = (double*)malloc(sizeof(double)*independent_nodes);
+  for (int i = 0; i < independent_nodes; ++i){
+    periodic_porosity[i] = 0;
   }
+
+  int* periodic_mapping = mesh->periodic_mapping;
   for (int i = 0; i < n_particles; ++i) {
     if(particle_element_id[i] == -1) continue;
     double sf[N_SF];
     shape_functions(&particle_uvw[i*D],sf);
     const int *el = &mesh->elements[particle_element_id[i]*N_N];
     for (int j=0; j<N_N;++j)
-      porosity[el[j]] += sf[j]*particle_volume[i];
+      periodic_porosity[periodic_mapping[el[j]]] += sf[j]*particle_volume[i];
   }
+  for(int i = 0; i < mesh->n_nodes; ++i){
+    porosity[i] = periodic_porosity[periodic_mapping[i]];
+  }
+  free(periodic_porosity);
+
   for (int i = 0; i < mesh->n_nodes; ++i){
     if(node_volume[i]==0.){
       porosity[i] = 0.;
@@ -1914,9 +1931,9 @@ double *fluid_problem_element_size(const FluidProblem *p)
 {
   return p->element_size;
 }
-int *fluid_problem_parent_tag(const FluidProblem *p)
+int *fluid_problem_periodic_mapping(const FluidProblem *p)
 {
-  return p->mesh->parent_tag;
+  return p->mesh->periodic_mapping;
 }
 
 void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals)

fluid/mesh.c

@@ -267,26 +267,26 @@ Mesh *mesh_load_msh(const char *filename)
   }
   // Periodic
   m->n_periodic = gmsh_m->n_periodic;
-  m->parent_tag = malloc(sizeof(int)*m->n_nodes);
+  m->periodic_mapping = malloc(sizeof(int)*m->n_nodes);
   for(int i = 0; i < m->n_nodes; ++i){
-    m->parent_tag[i] = i;
+    m->periodic_mapping[i] = i;
   }
   for(int i = 0; i < m->n_periodic; ++i){
     Periodic *p = gmsh_m->periodic[i];
     for(int j = 0; j < p->n_nodes; ++j){
       int nid = get_node_id(p->nodes[2*j+1], nodes_map, m->n_nodes);
       int nidSlave = get_node_id(p->nodes[2*j+0], nodes_map, m->n_nodes);
-      m->parent_tag[nid] = m->parent_tag[nidSlave];
+      m->periodic_mapping[nid] = m->periodic_mapping[nidSlave];
     }
   }
   int *renumber = malloc(sizeof(int)*m->n_nodes);
   int c = 0;
   for (int i = 0; i < m->n_nodes; ++i) {
-    if (m->parent_tag[i] == i) 
+    if (m->periodic_mapping[i] == i)
       renumber[i] = c++;
   }
   for (int i = 0; i < m->n_nodes; ++i) {
-    m->parent_tag[i] = renumber[m->parent_tag[i]];
+    m->periodic_mapping[i] = renumber[m->periodic_mapping[i]];
   }
   free(renumber);
 
@@ -308,7 +308,7 @@ void mesh_free(Mesh *m)
   }
   free(m->boundary_names);
   free(m->interfaces);
-  free(m->parent_tag);
+  free(m->periodic_mapping);
   free(m);
 }
 

fluid/mesh.h

@@ -37,7 +37,7 @@ typedef struct {
   int n_interfaces;
   int *interfaces; // eid0, closure0, {(eid1, closure1),(-1,phys_id)}
   int n_periodic;
-  int *parent_tag; // n_nodes
+  int *periodic_mapping; // n_nodes
 } Mesh;
 
 Mesh *mesh_load_msh(const char *filename);

python/fluid.py

@@ -344,7 +344,7 @@ class FluidProblem :
         self._lib.fluid_problem_reset_boundary_force(self._fp)
         self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
 
-        periodic_map = self._get_matrix("parent_tag", self.n_nodes(), 1, c_int).flatten()
+        periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
         if self.sys is None :
             self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options)
         sys = self.sys

python/mesh.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/>.
 
 from . import gmshType
@@ -46,13 +46,14 @@ class Mesh() :
       self.dimension = dim
       self.physicals = physicals
       self.elements = []
-      
+
   def __init__(self, filename = None) :
     self.vertices = []
     self.physicals = [{}, {}, {}, {}]
     self.entities = []
     self.maxeid = 0
     self.useFormat3 = False
+    # self.useFormat4 = False
     if filename :
       self.read(filename)
 
@@ -114,25 +115,28 @@ class Mesh() :
             if float(l[0]) == 3.:
               self.useFormat3 = True
             elif int(float(l[0])) == 2 :
-              self.useFormat3 = False 
+              self.useFormat3 = False
+            # elif float(l[0]) == 4.1:
+            #   self.useFormat4 = True
             else :
               print("error : cannot read mesh format " + l[0])
             l = fin.readline()
-          elif w[0] == "$PhysicalNames" :
+          elif w[0] == "$PhysicalNames" : #Same as MSH2
             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:
+          elif w[0] == "$Entities" and self.useFormat3: #To adapt for MSH4
             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" :
+          elif w[0] == "$Nodes" : #To adapt for MSH4
             n = int(fin.readline())
             for i in range(n) :
               if self.useFormat3 :
@@ -141,7 +145,7 @@ class Mesh() :
                 (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" :
+          elif w[0] == "$Elements" : #To adapt for MSH4
             n = int(fin.readline())
             for i in range(n) :
               l = fin.readline().split()
@@ -160,7 +164,7 @@ class Mesh() :
                 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)) 
+              entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
               self.maxeid = max(self.maxeid, j)
           l = fin.readline()
 
@@ -180,7 +184,7 @@ class Mesh() :
   def newVertex(self, x, y, z) :
     self.vertices.append((x, y, z, len(self.vertices) + 1))
     return self.vertices[-1]
- 
+
   def newElement(self, t, partition, entity, vertices) :
     self.maxeid += 1
     entity.elements.append(Mesh.element(t, vertices, partition, self.maxeid))

testcases/periodicTest/grain.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 scontact
+from migflow import time_integration
+
+import numpy as np
+import os
+import subprocess
+import time
+import shutil
+import random
+import unittest
+
+outputdir = "outputGrain"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+subprocess.call(["gmsh", "-2", "meshGrain.geo","-clscale","1"])
+subprocess.call(["gmsh", "-2", "meshGrain2.geo","-clscale","1"])
+
+t = 0
+ii = 0
+
+r = 1e-2
+V = 0
+Vref = 0.7283018867924522
+#physical parameters
+g =  np.array([0,-9.81])                                      # gravity
+rho = 1000                                      # fluid density
+rhop = 2640
+nu = 1e-6                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 10                                     # final time
+
+#numerical parameters
+lcmin = .1                                  # mesh size
+dt = .01                                       # time step
+
+shutil.copy("meshGrain.msh", outputdir +"/meshGrain.msh")
+shutil.copy("meshGrain2.msh", outputdir +"/meshGrain2.msh")
+p = scontact.ParticleProblem(2)
+p.load_msh_boundaries("meshGrain2.msh", ["Top", "Bottom","Lateral"])
+p.add_particle((0,0.8),r,r**2*np.pi*rhop)
+p.write_vtk(outputdir,0,0)
+p = scontact.ParticleProblem(2)
+p.read_vtk(outputdir,0)
+p.write_vtk(outputdir, 0, 0)
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+
+outf = 5                                       # number of iterations between output files
+
+fluid = mbfluid.FluidProblem(2,g,nu*rho,rho)
+fluid.load_msh("meshGrain.msh")
+fluid.set_wall_boundary("Bottom")
+fluid.set_wall_boundary("Lateral")
+fluid.set_wall_boundary("Top",pressure=0)
+ii = 0
+t = 0
+
+
+fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces())
+
+fluid.export_vtk(outputdir,0,0)
+
+G = p.mass()*g
+ii = 0
+tic = time.time()
+while t < tEnd :
+    time_integration.iterate(fluid, p, dt,external_particles_forces=G)
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+        p.write_vtk(outputdir, ioutput, t)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
+V = p.state().v[0,1]

testcases/periodicTest/mesh2.geo

+L = 5;
+H = 1;
+y = 0;
+lc = 0.05;
+
+Point(1) = {0,0,0,lc};
+Point(2) = {L,0,0,lc};
+Point(3) = {L,H,0,lc};
+Point(4) = {0,H,0,lc};
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {4,3};
+Line(4) = {1,4};
+Line Loop(1) = {1,2,-3,-4};
+Plane Surface(1) = {1};
+Periodic Curve {2} = {4};
+
+Physical Line("Left") = {4};
+Physical Line("Right") = {2};
+Physical Line("Bottom") = {1};
+Physical Line("Top") = {3};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 2;

testcases/periodicTest/meshGrain.geo

+L = .2;
+H = 1;
+y = 0;
+lc = 0.1;
+
+Point(1) = {-L, H, 0,lc};
+Point(2) = {-L, -H, 0,lc};
+Point(3) = {L, -H, 0,lc};
+Point(4) = {L, H, 0,lc};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line Loop(1) = {1:4};
+Periodic Curve{2} = {-4};
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2};
+Physical Line("Top") = {4};
+Physical Line("Lateral") = {1,3};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 4.1;

testcases/periodicTest/meshGrain2.geo

+L = .2;
+H = 1;
+y = 0;
+lc = 0.1;
+
+Point(1) = {-L, H, 0,lc};
+Point(2) = {-L, -H, 0,lc};
+Point(3) = {L, -H, 0,lc};
+Point(4) = {L, H, 0,lc};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line Loop(1) = {1:4};
+Periodic Curve{2} = {4};
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2};
+Physical Line("Top") = {4};
+Physical Line("Lateral") = {1,3};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 2;

testcases/periodicTest/poiseuille.py

@@ -59,7 +59,7 @@ print('epsilon',epsilon)
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
-outf = 1                                       # number of iterations between output files
+outf = 10                                       # number of iterations between output files
 
 #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
 fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
@@ -69,10 +69,10 @@ fluid.load_msh("mesh.msh")
 fluid.set_wall_boundary("Bottom",velocity=[0,0])
 fluid.set_wall_boundary("Top",velocity=[0,0])
 
-fluid.set_strong_boundary("Bottom",0,0)
+# fluid.set_strong_boundary("Bottom",0,0)
 # if strong boundary on periodic line, it should be forced on both sides
-#fluid.set_strong_boundary("Right",2,0)
-#fluid.set_strong_boundary("Left",2,0)
+fluid.set_strong_boundary("Right",2,0)
+fluid.set_strong_boundary("Left",2,0)
 
 
 ii = 0

testcases/periodicTest/poiseuilleGrain.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
+from migflow import scontact
+
+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"])
+subprocess.call(["gmsh", "-2", "mesh2.geo","-clscale","2"])
+
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  np.array([0.01,0])                                      # gravity
+rho = 1000                                      # fluid density
+rhop = 9000                                     # particle density
+nu = 1e-3                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 200                                     # final time
+r = .25e-1                                        #particle radius
+
+
+#numerical parameters
+lcmin = .1                                  # mesh size
+dt = 1                                       # time step
+alpha = 1e-4                                    # stabilization coefficient
+epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
+print('epsilon',epsilon)
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+shutil.copy("mesh2.msh", outputdir +"/mesh2.msh")
+
+p = scontact.ParticleProblem(2)
+p.load_msh_boundaries("mesh2.msh", ["Bottom", "Top", "Right", "Left"])
+p.add_particle((0.5,0.5),r,r**2*np.pi*rhop)
+p.write_vtk(outputdir,0,0)
+p = scontact.ParticleProblem(2)
+p.read_vtk(outputdir,0)
+p.write_vtk(outputdir, 0, 0)
+print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("RHOP = %g" % rhop)
+
+outf = 1                                       # number of iterations between output files
+
+#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
+fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
+fluid.load_msh("mesh.msh")
+#fluid.set_open_boundary("Right",pressure = 0)
+#fluid.set_open_boundary("Left", pressure = 0)
+fluid.set_wall_boundary("Bottom",velocity=[0,0])
+fluid.set_wall_boundary("Top",velocity=[0,0])
+
+# fluid.set_strong_boundary("Bottom",0,0)
+# if strong boundary on periodic line, it should be forced on both sides
+fluid.set_strong_boundary("Right",2,0)
+fluid.set_strong_boundary("Left",2,0)
+
+
+ii = 0
+t = 0
+
+#set initial_condition
+fluid.set_particles(p.mass(), p.volume(), p.state(),p.contact_forces())
+fluid.export_vtk(outputdir,0,0)
+
+G = p.mass()*g
+ii = 0
+tic = time.time()
+while t < tEnd :
+    #Fluid solver
+    time_integration.iterate(fluid,p,dt,external_particles_forces=G)
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+        p.write_vtk(outputdir, ioutput, t)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
+
+s = fluid.solution()
+x = fluid.coordinates()
+# vel = (s[:,0]-1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2
+# vS = np.sum((1/(20*nu*rho)*x[:,1]*(1-x[:, 1]))**2)
+# print('Error', (vel.sum())**.5)
+#
+# if(vel.sum()**.5 > (vS**0.5)/50) :
+#     print("error is too large in Poiseuille")
+#     exit(1)