Add shear testcase

testcases/shear-2d/mesh.geo

+a = 0.1;
+L = a;
+H = a;
+y = 0;
+lc = 0.1*a;
+
+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};

testcases/shear-2d/shear.py

+# MigFlow - Copyright (C) <2010-2020>
+# <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/>.
+
+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
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+os.chdir(dir_path)
+outputdir = "outputShear2D_2"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
+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
+rhop = 2000
+r = 2e-3
+
+# Geometrical parameters
+H = 0.1
+L = 0.1
+
+# Numerical parameters
+# tEnd = 2                                         # final time
+tEnd = 20                                         # final time
+dt = 1e-3                                                # time step
+# dt = 10                                             # time step
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+outf = 5                                           # 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,mu,rho,petsc_solver_type="-pc_type lu")
+fluid.load_msh("mesh.msh")
+fluid.set_wall_boundary("Bottom",velocity=[0,0], pressure = 0)
+fluid.set_wall_boundary("Top",velocity=[1.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)
+
+
+# Particle Problem
+p = scontact.ParticleProblem(2)
+p.load_msh_boundaries("mesh.msh", ["Top", "Bottom", "Left", "Right"])
+for x in np.arange(r,L-r,2.1*r):
+    for y in np.arange(r,H-r,2.1*r):
+        R = r*(1-np.random.random()/4)
+        p.add_particle((x+R,y),R,R**2*np.pi*rhop)
+
+
+
+ii = 0
+t = 0
+
+#set initial_condition
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
+fluid.export_vtk(outputdir,0,0)
+p.write_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.time()
+# while t < tEnd :
+while t < tEnd:
+    #Fluid solver
+    time_integration.iterate(fluid,p,dt,min_nsub=5)
+
+    #Ensure particle teleportation
+    p.position()[:,0] = np.fmod(p.position()[:,0],L)
+
+    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))
\ No newline at end of file

testcases/shear-2d/shearFluid.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/>.
+
+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
+
+dir_path = os.path.dirname(os.path.realpath(__file__))
+os.chdir(dir_path)
+outputdir = "outputShearFluid"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
+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 = 2                                         # final time
+tEnd = 100                                         # final time
+
+#numerical parameters
+# dt = 1e-3                                                # time step
+dt = 10                                             # time step
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+outf = 2                                           # 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,mu,rho,petsc_solver_type="-pc_type lu")
+fluid.load_msh("mesh.msh")
+fluid.set_wall_boundary("Bottom",velocity=[0,0], pressure = 0)
+fluid.set_wall_boundary("Top",velocity=[1.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.export_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.time()
+# while t < tEnd :
+while t < tEnd:
+    #Fluid solver
+    time_integration.iterate(fluid,None,dt)
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
\ No newline at end of file