Periodic case for grains in the drag experiment

python/marblesbag/scontact2.py

@@ -97,6 +97,9 @@ class ParticleProblem :
 
     def write_vtk(self, odir, i, t) :
         lib.particleProblemWriteVtk(self._p, odir.encode(), i)
+        
+    def periodic(self, idnum, zonexmin,zoneymin,zonexmax,zoneymax,newX,newY) :
+        lib.checkFree(self._p, c_int(idnum), c_double(zonexmin), c_double(zoneymin), c_double(zonexmax), c_double(zoneymax), c_double(newX), c_double(newY), c_void_p(), c_int(0))
 
     def read_vtk(self,dirname,i) :
         lib.particleProblemReadVtk(self._p, dirname.encode(),i)

scontact/Untitled

+      double *omega = &p->omega;
+      //printf("%f\n",omega[1]);
+      //break;
+      double vt = (v0[0] - v1[0]) * c->t[0] + (v0[1] - v1[1]) * c->t[1] +c->r0*omega[c->o0] +c->r1*omega[c->o1]; 
+      double ct;
+      //double vt = 1.0;
+      if(vt>0.0){
+        if(vt>=0.2666*dp){ct = 0.1333*dp;}
+        else{ct = vt;}
+        //coordAdd2(&p->omega[c->o0], -(2.5*c->r0)*ct*c->a0); //*c->a0
+        //coordAdd2(&p->omega[c->o1], -(2.5*c->r1)*ct*c->a1); //*c->a0
+        coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, c->t); //*c->a0
+        coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, c->t); //*c->a1
+        
+      }
+      else if(vt<0.0){
+        if(vt<=-0.2666*dp){ct = -0.1333*dp;}
+        else{ct = vt;}
+        //coordAdd2(&p->omega[c->o0], -(2.5*c->r0)*ct*c->a0); //*c->a0
+        //coordAdd2(&p->omega[c->o1], -(2.5*c->r1)*ct*c->a1); //*c->a0
+        coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*c->a0, c->t); //*c->a0
+        coordAdd(&p->velocity[c->o1 * DIMENSION], ct*c->a1, c->t); //*c->a1
+        
+      }
\ No newline at end of file

scontact/Untitled4

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) 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 marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+#grains creation and placement + physical boundaries definition
+'''r is the radius of the grains
+rout is the outer radius of the geometry
+rin is the inner radius of the geometry
+rhop is the particles density''' 
+def genInitialPosition(filename, r, rin, rhop) :
+    p = scontact2.ParticleProblem()
+    #Loading of the mesh.msh file specifying physical boundaries name
+    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
+    #Definition of the points where the grains are located
+    x = np.arange(-0.025+r, 0.025-r, 2*r)
+    y = np.arange(-0.03+r, 0.035-r, 2*r)
+    x, y = np.meshgrid(x, y)
+    R2 = x**2 + y**2
+    #condition to be outside the inner boundary
+    keep = R2 > (rin + r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+      rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
+      #Addition of an particle object at each point
+      p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
+    p.write_vtk(filename,0,0)
+
+
+t = 0
+ii = 0
+#physical parameters
+g = 9.81                                # gravity
+rho = 1.253e3                           # fluid density
+rhop = 1000                            # grains density
+nu = 1e-1                               # kinematic viscosity
+tEnd = 50                               # final time
+
+#numerical parameters
+h = 0.002                               # mesh size
+dt = 1e-2                               # time step
+epsilon = 1e-5                          # stabilisation parametre
+
+#geometry parameters
+rin = 0.0064                            # inner radius
+r = 397e-6                           # grains radius
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#Object particles creation
+genInitialPosition(outputdir, r, rin, rhop)
+p = scontact2.ParticleProblem()
+p.read_vtk(outputdir,0)
+
+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)
+#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
+strong_boundaries = [("Outer",2,0.),("Outer",1,0),("Outer",0,0),("Inner",1,0),("Inner",2,0),("Inner",1,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+p.velocity()[:,1] += -0.05
+tic = time.clock()
+wut = p.position()[:,1]<-0.02
+Displace = np.nonzero(p.position()[:,1]<-0.02)
+print(lel)
+#fluid.export_vtk(outputdir, 0, 0)
+#Computation loop
+while t < tEnd :
+    #Fluid solver
+    #fluid.implicit_euler(dt)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    pmax = np.max(p.position()[:,1])
+    lowY = pmax+r1
+    highY = lowY + 5*r1
+    #number of sub time step
+    nsub = max(10, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        tol = 1e-6
+        p.iterate(dt/nsub, forces, tol)
+    #Localisation of the grains in the fluid
+    #fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

scontact/Untitled5

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) 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 marblesbag import fluid as fluid
+from marblesbag import scontact2
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+#grains creation and placement + physical boundaries definition
+'''r is the radius of the grains
+rout is the outer radius of the geometry
+rin is the inner radius of the geometry
+rhop is the particles density''' 
+def genInitialPosition(filename, r, rin, rhop) :
+    p = scontact2.ParticleProblem()
+    #Loading of the mesh.msh file specifying physical boundaries name
+    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
+    #Definition of the points where the grains are located
+    x = np.arange(-0.025+r, 0.025-r, 2*r)
+    y = np.arange(-0.03+r, 0.035-r, 2*r)
+    x, y = np.meshgrid(x, y)
+    R2 = x**2 + y**2
+    #condition to be outside the inner boundary
+    keep = R2 > (rin + r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+      rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
+      #Addition of an particle object at each point
+      p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
+    p.write_vtk(filename,0,0)
+
+
+t = 0
+ii = 0
+#physical parameters
+g = 9.81                                # gravity
+rho = 1.253e3                           # fluid density
+rhop = 1000                            # grains density
+nu = 1e-1                               # kinematic viscosity
+tEnd = 50                               # final time
+
+#numerical parameters
+h = 0.002                               # mesh size
+dt = 1e-2                               # time step
+epsilon = 1e-5                          # stabilisation parametre
+
+#geometry parameters
+rin = 0.0064                            # inner radius
+r = 397e-6                           # grains radius
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#Object particles creation
+genInitialPosition(outputdir, r, rin, rhop)
+p = scontact2.ParticleProblem()
+p.read_vtk(outputdir,0)
+
+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)
+#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
+strong_boundaries = [("Outer",2,0.),("Outer",1,0),("Outer",0,0),("Inner",1,0),("Inner",2,0),("Inner",1,0.)]
+fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+fluid.export(outputdir,0,0)
+p.velocity()[:,1] += -0.05
+tic = time.clock()
+r1 = np.min(p.r())
+#fluid.export_vtk(outputdir, 0, 0)
+#Computation loop
+while t < tEnd :
+    #Fluid solver
+    #fluid.implicit_euler(dt)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    pmax = np.max(p.position()[:,1])
+    lowY = pmax+r1
+    lowX = -0.025 + r1
+    highY = lowY + 5*r1
+    highX = 0.025 - r1
+    #number of sub time step
+    nsub = max(10, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        tol = 1e-6
+        p.iterate(dt/nsub, forces, tol)
+        if t>10 :
+          Displace = np.nonzero(p.position()[:,1]<-0.02)
+          for j in range(len(Displace))
+            Xnew = 
+            p.periodic(Displace[i],lowX,lowY,highX,highY,Xnew,Ynew)
+    #Localisation of the grains in the fluid
+    #fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))

scontact/scontact.h

@@ -24,6 +24,7 @@
 #ifndef _SCONTACT_H_
 #define _SCONTACT_H_
 #include <stdlib.h>
+#include "quadtree.h"
 
 typedef struct _ParticleProblem ParticleProblem;
 
@@ -33,6 +34,7 @@ void particleProblemLoad(ParticleProblem *p, const char *filename);
 void particleProblemWrite(const ParticleProblem *p, const char *filename);
 void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit);
 void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
+void checkFree(ParticleProblem *p, int id,double zonexmin, double zoneymin, double zonexmax, double zoneymax, double newX, double newY, Cell *tree, int iter);
 double particleProblemMaxDt(const ParticleProblem *p, double alert);
 void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m);
 size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag);
@@ -45,6 +47,7 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIME
 double *particleProblemDisk(ParticleProblem *p);
 double *particleProblemSegment(ParticleProblem *p);
 double *particleProblemParticle(ParticleProblem *p);
+double *particleProblemOmega(ParticleProblem *p);
 double *particleProblemVelocity(ParticleProblem *p);
 double *particleProblemPosition(ParticleProblem *p);
 double *particleProblemExternalForces(ParticleProblem *p);

testcases/couette-2d/melangeur.py

@@ -43,7 +43,7 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
     #Loading of the mesh.msh file specifying physical boundaries name
     p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
     #Definition of the points where the grains are located
-    x = np.arange(rout, -rout, 2 * -r)
+    x = np.arange(rout, -rout, -2*r)
     x, y = np.meshgrid(x, x)
     R2 = x**2 + y**2
     #condition to be inside the outer boundary
@@ -65,7 +65,6 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
 
 t = 0
 ii = 0
-
 #physical parameters
 g =  0                                  # gravity
 rho = 1.253e3                           # fluid density

testcases/couette-2d/mesh.geo

@@ -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;
+

testcases/depot-2d/Classical/mesh.geo

+L = 0.025;
+H = 0.06;
+y = 0;
+lc = 0.0015;
+
+Point(1) = {-L, H, 0, lc};
+Point(2) = {-L, -H, 0, lc};
+Point(3) = {L, -H, 0, lc};
+Point(4) = {L, H, 0, lc};
+
+Point(5) = {-8*L/9.,-H,0,lc};
+Point(6) = {-7*L/9.,-H,0,lc};
+Point(7) = {-6*L/9.,-H,0,lc};
+Point(8) = {-5*L/9.,-H,0,lc};
+Point(9) = {-4*L/9.,-H,0,lc};
+Point(10) = {-3*L/9.,-H,0,lc};
+Point(11) = {-2*L/9.,-H,0,lc};
+Point(12) = {-L/9.,-H,0,lc};
+
+Point(21) = {0.,-H,0,lc};
+Point(22) = {L/9.,-H,0,lc};
+Point(23) = {2*L/9.,-H,0,lc};
+Point(24) = {3*L/9.,-H,0,lc};
+Point(25) = {4*L/9.,-H,0,lc};
+Point(26) = {5*L/9.,-H,0,lc};
+Point(27) = {6*L/9.,-H,0,lc};
+Point(28) = {7*L/9.,-H,0,lc};
+Point(29) = {8*L/9.,-H,0,lc};
+
+Point(13) = {-8*L/9.,H,0,lc};
+Point(14) = {-7*L/9.,H,0,lc};
+Point(15) = {-6*L/9.,H,0,lc};
+Point(16) = {-5*L/9.,H,0,lc};
+Point(17) = {-4*L/9.,H,0,lc};
+Point(18) = {-3*L/9.,H,0,lc};
+Point(19) = {-2*L/9.,H,0,lc};
+Point(20) = {-L/9.,H,0,lc};
+
+Point(30) = {0.,H,0,lc};
+Point(31) = {L/9.,H,0,lc};
+Point(32) = {2*L/9.,H,0,lc};
+Point(33) = {3*L/9.,H,0,lc};
+Point(34) = {4*L/9.,H,0,lc};
+Point(35) = {5*L/9.,H,0,lc};
+Point(36) = {6*L/9.,H,0,lc};
+Point(37) = {7*L/9.,H,0,lc};
+Point(38) = {8*L/9.,H,0,lc};
+
+Line(1) = {1, 2};
+
+Line(2) = {2, 5};
+Line(3) = {5, 6};
+Line(4) = {6, 7};
+Line(5) = {7, 8};
+Line(6) = {8, 9};
+Line(7) = {9, 10};
+Line(8) = {10, 11};
+Line(9) = {11, 12};
+Line(10) = {12, 21};
+Line(11) = {21, 22};
+Line(12) = {22, 23};
+Line(13) = {23, 24};
+Line(14) = {24, 25};
+Line(15) = {25, 26};
+Line(16) = {26, 27};
+Line(17) = {27, 28};
+Line(18) = {28, 29};
+Line(19) = {29, 3};
+
+Line(20) = {3, 4};
+
+Line(21) = {4, 38};
+Line(22) = {38, 37};
+Line(23) = {37, 36};
+Line(24) = {36, 35};
+Line(25) = {35, 34};
+Line(26) = {34, 33};
+Line(27) = {33, 32};
+Line(28) = {32, 31};
+Line(29) = {31, 30};
+Line(30) = {30, 20};
+Line(31) = {20, 19};
+Line(32) = {19, 18};
+Line(33) = {18, 17};
+Line(34) = {17, 16};
+Line(35) = {16, 15};
+Line(36) = {15, 14};
+Line(37) = {14, 13};
+Line(38) = {13, 1};
+
+
+
+Line Loop(1) = {1:38};
+
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2,4,6,8,10,12,14,16,18};
+Physical Line("BottomOut") = {3,5,7,9,11,13,15,17,19};
+Physical Line("Lateral") = {1,20};
+Physical Line("Top") = {38,36,34,32,30,28,26,24,22};
+Physical Line("TopOut") = {21,23,25,27,29,31,33,35,37};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+Mesh.MshFileVersion = 2;
+Mesh.MshFileVersion = 2;

testcases/depot-2d/dep.py

-# Marblesbag - Copyright (C) <2010-2018>
+# MigFlow - Copyright (C) <2010-2018>
 # <Universite catholique de Louvain (UCL), Belgium
 #  Universite de Montpellier, France>
 # 	
-# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# List of the contributors to the development of MigFlow: see AUTHORS file.
 # Description and complete License: see LICENSE file.
 # 	
-# This program (Marblesbag) is free software: 
+# 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.
@@ -28,97 +28,84 @@ import os
 import time
 import shutil
 import random
-
-#grains creation and placement + physical boundaries definition
-'''r is the radius of the grains
-ly is the height of the box
-rhop is the particles density'''
-def genInitialPosition(filename, r, ly, rhop) :
+'''
+Creation of the deposit used for the avalanch test cases
+'''
+
+'''
+genInitialPosition is a function that sets all the grains centre positions and create the grains objects to add in the computing structure
+Args:    -filename    : name of the output file
+         -r           : mean radius of the grains
+         -ly          : domain height
+         -lx          : domain width 
+         -rhop        : grains density        
+'''
+def genInitialPosition(filename, r, rin, rhop) :
     p = scontact2.ParticleProblem()
     #Loading of the mesh.msh file specifying physical boundaries name
-    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","BottomOut","TopOut"])
+    p.load_msh_boundaries("mesh.msh", ["Outer", "Inner"])
     #Definition of the points where the grains are located
     x = np.arange(-0.025+r, 0.025-r, 2*r)
-    y = np.arange(0.06-r, .04, -2*r)
+    y = np.arange(-0.03+r, 0.12-r, 2*r)
     x, y = np.meshgrid(x, y)
-    x = x.flat
-    y = y.flat
-    for i in range(len(x)) :
-        rhop = random.choice([1100,1200,1300])
-        #Addition of an particle object at each point
-        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
+    R2 = x**2 + y**2
+    #condition to be outside the inner boundary
+    keep = R2 > (rin + r)**2
+    x = x[keep]
+    y = y[keep]
+    for i in range(x.shape[0]) :
+      rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
+      #Addition of an particle object at each point
+      p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1);
     p.write_vtk(filename,0,0)
 
-
+#Define output directory
 outputdir = "output"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
-
-t = 0
-ii = 0
+#Initial time and iteration
+t         =  0
+ii        =  0
 
 
 #physical parameters
-r = 1e-4                            # radius
-ly = 5e-2                           # box height
-g =  -9.81                          # gravity
-rho = 1000                          # fluid density
-rhop = 1500                         # grains density
-nu = 1e-6                           # kinematic viscosity
-tEnd = 100                          # final time
+g         = -9.81                                     #gravity
+rhop      =  1500                                     #grains density
+r         =  5e-4                                     #radius
+rin       =  0.0064                                   #grains area width
+tEnd      =  50                                        #final time
 
 #numerical parameters
-lcmin = 0.001                       # mesh size
-dt = 1e-3                           # time step
-alpha = 2.5e-5                      # stabilization coefficient
-epsilon = alpha*lcmin**2 /nu        # stabilization parameter
-print('epsilon',epsilon)
-
-shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+dt        =  1e-3                                     #time step
 
-#Object particles creation
-genInitialPosition(outputdir, r, ly, rhop)
-p = scontact2.ParticleProblem()
-p.read_vtk(outputdir,0)
+#Initialise grains
+genInitialPosition(outputdir, r, rin, rhop)
+p         =  scontact2.ParticleProblem()
+p.read_vtk(outputdir,20)
 
 print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
 print("RHOP = %g" % rhop)
-outf = 200                          # number of iterations between output files
-
-ii = 0
-
-#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
-#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
-strong_boundaries = [("Top",2,0.),("TopOut",1,0),("Top",1,0),("BottomOut",1,0),("Bottom",1,0),("Lateral",0,0.),("Lateral",1,0.)]
-fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
-fluid.export(outputdir,0,0)
-
+outf      =  200                                      #number of iterations between output files
 
 tic = time.clock()
+forces = np.zeros_like(p.velocity())
 #Computation loop
-fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
-while t < tEnd :
-    #Fluid solver
-    fluid.implicit_euler(dt)
-    forces = fluid.compute_node_force(dt)
+while t < tEnd : 
     #Computation of the new velocities
-    vn = p.velocity() + forces * dt / p.mass()
+    forces[:,1] = g*p.mass()[0]-p.velocity()[:,1]
+    vn = p.velocity() + forces * dt / p.mass() 
     vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
-    #number of sub time step
-    nsub = max(10, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    #number of contact sub time step
+    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
     print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
     #Contact solver
     for i in range(nsub) :
-        tol = 1e-8
-        p.iterate(dt/nsub, forces, tol)
-    #Localisation of the grains in the fluid
-    fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+        p.iterate(dt/nsub, forces) 
     t += dt
-    #Output files writting
+
+    #Output files writing
     if ii %outf == 0 :
-        ioutput = int(ii/outf) + 1
+        ioutput = int(ii/outf) + 1 + 20
         p.write_vtk(outputdir, ioutput, t)
-        fluid.export_vtk(outputdir, t, ioutput)
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))