adapt_mesh

d627daae · Michel Henry · 5dedbeec · d627daae · d627daae · d627daae
Commit d627daae authored Dec 30, 2021 by Michel Henry
--- a/fluid/fluid_problem.c
+++ b/fluid/fluid_problem.c
--- a/fluid/fluid_problem.h
+++ b/fluid/fluid_problem.h
@@ -113,7 +113,9 @@ void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
 void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
 void fluid_problem_free(FluidProblem *problem);
 FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f, int stability);
-void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume);
+void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume,
+                              int *idx, int *idx_c, int *idx_p, int *idx_mesh, int n_idx, int n_idx_c, int n_idx_p, int n_idx_mesh);
+
 void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh);
 int *fluid_problem_particle_element_id(FluidProblem *problem);
 int fluid_problem_n_particles(FluidProblem *problem);
@@ -129,13 +131,18 @@ double *fluid_problem_concentration_dg(const FluidProblem *p);
 double *fluid_problem_element_size(const FluidProblem *p);
 void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
 int fluid_problem_local_size(FluidProblem *p);
-void fluid_problem_fields_local_map(FluidProblem *p, int *idx);
-void fluid_problem_mesh_local_map(FluidProblem *p, int *idx);
-void fluid_problem_concentration_local_map(FluidProblem *p, int *idx);
-void fluid_problem_porosity_local_map(FluidProblem *p, int *idx);
-void fluid_problem_local_map(FluidProblem *p, int *idx);
+void fluid_problem_fields_local_map(FluidProblem *p, int *idx, int n);
+void fluid_problem_mesh_local_map(FluidProblem *p, int *idx, int n);
+void fluid_problem_concentration_local_map(FluidProblem *p, int *idx, int n);
+void fluid_problem_porosity_local_map(FluidProblem *p, int *idx, int n);
+int *fluid_problem_idx(const FluidProblem *p);
 int fluid_problem_get_ndof(const FEFields *fields, const Mesh *mesh);
 void fluid_problem_set_problem(FluidProblem *problem);
+void fluid_problem_get_size_field(const FluidProblem *problem, double *mesh_size);
+
+double *fluid_problem_particle_position(FluidProblem *problem);
+double *fluid_problem_particle_velocity(FluidProblem *problem);
+double *fluid_problem_particle_contacts(FluidProblem *problem);

 // mesh data (for i/o)
 int fluid_problem_n_nodes(const FluidProblem *p);
@@ -146,7 +153,8 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
 int fluid_problem_n_mesh_boundaries(const FluidProblem *p);
 void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize);
 void fluid_problem_mesh_boundary_interfaces(const FluidProblem *p, int bid, int *binterfaces);
-
+int *mesh_elements(const Mesh *m);
+int mesh_n_elements(const Mesh *m);

 void fluid_problem_set_stab_param(FluidProblem *p, double  stab_param);
 void fluid_problem_set_reduced_gravity(FluidProblem *p, int reduced_gravity);

--- a/fluid/fluid_problem_fem.c
+++ b/fluid/fluid_problem_fem.c
@@ -201,6 +201,7 @@ void fe_fields_free(FEFields *fields) {
    fe_element_free(fields->element[i]);
  }
  free(fields->element);
+  free(fields->map);
  free(fields);
 }

@@ -220,52 +221,6 @@ void fe_fields_df(const FEFields *fields, const double *xi, const double dxidx[D
  }
 }

-void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]) { // unused
-  int c[fields->n];
-  c[0] = 0;
-  for (int i = 1; i < fields->n; ++i) {
-    c[i] = c[i-1]+fields->element[i-1]->nlocal;
-  }
-  int l = 0;
-  int n_n = DIMENSION+1;
-
-  for (int i = 0; i < fields->n; ++i) {
-    for (int j = 0; j < fields->element[i]->n[0]; ++j){
-      for (int k = 0; k < n_n; ++k) {
-        map[c[i]++] = mesh->elements[iel*n_n+k]*fields->nd[0]+l;
-      }
-      l++;
-    }
-  }
-// #if DIMENSION == 2
-  // int n_e = 3;
-// #else
-  // int n_e = 6;
-// #endif
-  // int shift = l*mesh->n_nodes;
-  // int edges[n_e  *mesh->n_elements];
-  // mesh_gen_edges_by_element(mesh, edges);
-  // todo add edges (and faces in 3D)
-  // l = 0;
-  // for (int i = 0; i < fields->n; ++i){
-    // for(int j = 0; j < fields->element[i]->n[1]; ++j){
-      // for(int k = 0; k < n_e; ++k){
-        // map[c[i]++] = shift+edges[iel*n_e+k]*fields->nd[1]+l;
-      // }
-      // l++;
-    // }
-  // }
-  // unused
-  if (fields->n == 3) {
-    int *el = mesh->elements + 3*iel;
-    int map2[9] = {
-      el[0]*3+0, el[1]*3+0, el[2]*3+0,
-      el[0]*3+1, el[1]*3+1, el[2]*3+1,
-      el[0]*3+2, el[1]*3+2, el[2]*3+2,
-    };
-  }
-}
-
 void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]) {
  int *map = &(fields->map[iel*fields->local_size]);
  int c = 0;

--- a/fluid/fluid_problem_fem.h
+++ b/fluid/fluid_problem_fem.h
@@ -33,7 +33,6 @@ void fe_fields_eval_grad_sf(const FEFields *fields, const Mesh *mesh, int iel, c
 void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]);
 void fe_fields_eval_grad(const FEFields *fields, const double *xi, const double dxidx[D][D], const Mesh *mesh, int iel, const double *v, double s[], double ds[][DIMENSION]);
 void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]);
-void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]);
 void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, const double f1[][D], const double *sf, const double dsf[][D], double jw, double *local_vector);
 static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *f00, const double f01[][D], const double f10[][D], const double f11[][D][D], const double *sf, const double dsf[][D], double jw, double *local_matrix) {
  int c = 0;

--- a/fluid/mesh_find.c
+++ b/fluid/mesh_find.c
@@ -130,6 +130,7 @@ void mesh_tree_particle_to_mesh(MeshTree *mt, size_t nParticles, double *particl
    vectorClear(found);
  }
  vectorFree(found);
+
 }



--- a/paraview_utils/Makefile
+++ b/paraview_utils/Makefile
@@ -19,6 +19,6 @@
 # along with this program (see COPYING and COPYING.LESSER files).  If not, 
 # see <http://www.gnu.org/licenses/>.

-MigFlow.xml : python_filter_generator.py contacts-extract.py ids-extract.py concentration-dg.py names.py stress-tensor.py p2.py
-	python python_filter_generator.py $@ contacts-extract.py ids-extract.py concentration-dg.py names.py stress-tensor.py p2.py
+MigFlow.xml : python_filter_generator.py contacts-extract.py ids-extract.py concentration-dg.py names.py stress-tensor.py p2.py dynamic_pressure.py
+	python3 python_filter_generator.py $@ contacts-extract.py ids-extract.py concentration-dg.py names.py stress-tensor.py p2.py dynamic_pressure.py

--- a/paraview_utils/MigFlow.xml
+++ b/paraview_utils/MigFlow.xml
@@ -407,6 +407,64 @@
     </StringVectorProperty>


+   <Hints>
+     <ShowInMenu category="MigFlow"/>
+   </Hints>
+
+   </SourceProxy>
+</ProxyGroup>
+ <ProxyGroup name="filters">
+   <SourceProxy name="MFPDYN" class="vtkPythonProgrammableFilter" label="MigFlow dynamic pressure">
+
+     <Documentation
+       long_help="Remove the hydrostatic pressure from the pressure field"
+       short_help="Remove the hydrostatic pressure from the pressure field">
+     </Documentation>
+
+
+     <InputProperty
+       name="Input"
+       command="SetInputConnection">
+         <ProxyGroupDomain name="groups">
+           <Group name="sources"/>
+           <Group name="filters"/>
+         </ProxyGroupDomain>
+
+         <DataTypeDomain name="input_type">
+           <DataType value="vtkUnstructuredGrid"/>
+         </DataTypeDomain>
+     </InputProperty>
+
+
+
+
+
+
+     <!-- Output data type: "vtkUnstructuredGrid" -->
+     <IntVectorProperty command="SetOutputDataSetType"
+                        default_values="4"
+                        name="OutputDataSetType"
+                        number_of_elements="1"
+                        panel_visibility="never">
+       <Documentation>The value of this property determines the dataset type
+       for the output of the programmable filter.</Documentation>
+     </IntVectorProperty>
+
+
+     <StringVectorProperty
+       name="Script"
+       command="SetScript"
+       number_of_elements="1"
+       default_values="import numpy as np&#xA;ipt = inputs[0].Points #coord&#xA;nel = inputs[0].CellTypes.shape[0]&#xA;p = inputs[0].PointData[&quot;pressure&quot;]&#xA;ind = np.argsort(p)&#xA;#tri = np.copy(inputs[0].Cells).reshape(nel,-1)&#xA;#nn = tri.shape[1]-1&#xA;pmax = p[ind[-1]]&#xA;pmin = p[ind[0]]&#xA;ymax = ipt[ind[-1],1]&#xA;ymin = ipt[ind[0],1]&#xA;ph = ((pmax-pmin)*ipt[:,1] + pmin*ymax-pmax*ymin)/(ymax-ymin)&#xA;p_dyn = p - ph&#xA;output.Points = ipt&#xA;output.PointData.append(p_dyn,&quot;dynamic_pressure&quot;)&#xA;"
+       panel_visibility="advanced">
+       <Hints>
+        <Widget type="multi_line"/>
+      </Hints>
+     <Documentation>This property contains the text of a python program that
+     the programmable source runs.</Documentation>
+     </StringVectorProperty>
+
+
   <Hints>
     <ShowInMenu category="MigFlow"/>
   </Hints>

--- a/paraview_utils/show_problem.py
+++ b/paraview_utils/show_problem.py
@@ -14,27 +14,34 @@ animationScene = GetAnimationScene()
 animationScene.PlayMode = 'Snap To TimeSteps'
 renderView = GetActiveViewOrCreate('RenderView')
 renderView.InteractionMode = '2D'
+renderView.Background = [1.0, 1.0, 1.0]


 # Show the fluid pressure
 def showFluid(fluid, renderView):
-    fluidDisplay = Show(fluid, renderView)
-    fluidDisplay.Representation = 'Surface'
-    ColorBy(fluidDisplay, ('POINTS', 'pressure'))
-    fluidDisplay.SetScalarBarVisibility(renderView, True)
-    RenameSource('Fluid', fluid)
    Hide(fluid, renderView)
-    return fluidDisplay
+    RenameSource('Fluid', fluid)
+    dynPressure = MigFlowdynamicpressure(Input=fluid)
+    RenameSource('Dynamic Pressure', dynPressure)
+    dynPressureDisplay = Show(dynPressure, renderView)
+    dynPressureDisplay.Representation = 'Surface'
+    ColorBy(dynPressureDisplay, ('POINTS', 'dynamic_pressure'))
+
+    # show color bar/color legend
+    dynPressureDisplay.SetScalarBarVisibility(renderView, True)
+    colorbar = GetScalarBar(GetColorTransferFunction('dynamic_pressure'), renderView)
+    colorbar.TitleColor = [0, 0, 0]
+    colorbar.LabelColor = [0, 0, 0]
+
+
+

 # Show the fluid velocity
 def showVelocity(fluid, renderView):
    try:
        velocity = MigFlowp2velocity(fluid)
-        velocityDisplay = Show(velocity, renderView)
-        ColorBy(velocityDisplay, ('POINTS', 'velocity', 'Magnitude'))
-        animationScene.GoToLast()
-        velocityDisplay.RescaleTransferFunctionToDataRange(True, False)
-        RenameSource('velocity', velocity)
+        Hide(velocity, renderView)
+        RenameSource('Velocity', velocity)

    except NameError:
        print("MigFlow p2 velocity filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.")
@@ -50,6 +57,9 @@ def showParticles(particleProblem, renderView):
    glyph.ScaleArray = ['POINTS', 'Radius']
    glyph.ScaleFactor = 2.0
    glyph.GlyphMode = 'All Points'
+    glyphDisplay = GetDisplayProperties(glyph, view=renderView)
+    glyphDisplay.AmbientColor = [0.0, 0.0, 0.0]
+    glyphDisplay.DiffuseColor = [0.0, 0.0, 0.0]
    Show(glyph, renderView)
    RenameSource('Particles', particles)
    RenameSource('Glyph', glyph)
@@ -90,17 +100,6 @@ def showContacts(particleProblem, renderView):
        print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.")
        return None

-def showFluidVelocity(fluid, renderView):
-    fluidVelocity = Calculator(Input = fluid)
-    fluidVelocity.ResultArrayName = 'FluidVelocity'
-    fluidVelocity.Function = 'velocity/porosity'
-    RenameSource('Fluid Velocity', fluidVelocity)
-    fluidVelocityDisplay = Show(fluidVelocity, renderView)
-    ColorBy(fluidVelocityDisplay, ('POINTS', 'FluidVelocity', 'Magnitude'))
-    fluidVelocityDisplay.SetScalarBarVisibility(renderView, True)
-    Hide(fluidVelocity)
-    return fluidVelocity
-



@@ -118,13 +117,10 @@ if isfile(particleProblemFile):

 if isfile(fluidFile):
    fluid = PVDReader(FileName=fluidFile)
-    fluidDisplay = showFluid(fluid, renderView)
+    showFluid(fluid, renderView)
    showVelocity(fluid, renderView)
    if isfile(particleProblemFile):
        Hide(fluid, renderView)
-    else:
-        animationScene.GoToLast()
-        fluidDisplay.RescaleTransferFunctionToDataRange(False, True)


 # Show the first time step

--- a/python/fluid.py
+++ b/python/fluid.py
@@ -195,14 +195,14 @@ class FluidProblem :
        """
        mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
        self._lib.fluid_problem_set_mesh(self._fp, mesh)
-        idx = self.get_numbering([1]*self.dimension())
-        self._lib.fluid_problem_mesh_local_map(self._fp, _np2c(idx, dtype=np.int32))
-        idx = self.get_numbering(self.order_fields)
-        idx_concentration = self.get_numbering(self.order_concentration)
-        idx_porosity = self.get_numbering(self.order_porosity)
-        self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
-        self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
-        self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
+        idx_mesh = self.get_numbering(self.elements(),[1]*self.dimension())
+        self._lib.fluid_problem_mesh_local_map(self._fp, _np2c(idx_mesh, dtype=np.int32), c_int(idx_mesh.size))
+        idx = self.get_numbering(self.elements(),self.order_fields)
+        idx_concentration = self.get_numbering(self.elements(),self.order_concentration)
+        idx_porosity = self.get_numbering(self.elements(),self.order_porosity)
+        self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32), c_int(idx.size))
+        self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32), c_int(idx_porosity.size))
+        self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32), c_int(idx_porosity.size))
        self._lib.fluid_problem_set_problem(self._fp)
        self.sys = None
        gmsh.model.remove()
@@ -335,7 +335,20 @@ class FluidProblem :
        """
        self._lib.fluid_problem_adapt_gen_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin))
        mesh = _load_msh("adapt/mesh.msh", self._lib, self.dimension())
-        self._lib.fluid_problem_adapt_mesh(self._fp, mesh, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
+        f = self._lib.mesh_n_elements
+        f.restype = c_int
+        n_elements = f(mesh)
+        f = getattr(self._lib,"mesh_elements")
+        f.restype = POINTER(c_int)
+        elements= np.ctypeslib.as_array(f(mesh),shape=(n_elements, self._dim+1))
+
+        idx = self.get_numbering(elements, self.order_fields)
+        idx_concentration = self.get_numbering(elements, self.order_concentration)
+        idx_porosity = self.get_numbering(elements, self.order_porosity)
+        idx_mesh = self.get_numbering(elements, [1]*self._dim)
+        self._lib.fluid_problem_adapt_mesh(self._fp, mesh, c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume),
+         _np2c(idx, np.int32), _np2c(idx_concentration, np.int32), _np2c(idx_porosity, np.int32), _np2c(idx_mesh, np.int32),
+         c_int(idx.size), c_int(idx_concentration.size), c_int(idx_porosity.size), c_int(idx_mesh.size))
        self.sys = None

    def mesh_boundaries(self) :
@@ -380,7 +393,7 @@ class FluidProblem :
        isp2p1 = np.max(self.order_fields) == 2
        if isp2p1 :
            field_data.append((b"velocity", self.velocity()))
-            idx = self.get_numbering([2])
+            idx = self.get_numbering(self.elements(),[2])
            cell_data.append(("velocity_map", idx.reshape(self.n_elements(),-1)))
        else:
            data.append(("velocity", self.velocity()))
@@ -409,13 +422,15 @@ class FluidProblem :
                _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
                _np2c(data["parent_node_id"],np.int32) if "parent_node_id"  in data else None))
        self._lib.fluid_problem_set_mesh(self._fp, _mesh)
-        ndof = np.max(self.get_numbering(self.order_fields))+1
+        ndof = np.max(self.get_numbering(self.elements(),self.order_fields))+1
        id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
            + self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
            np.arange(self.n_fields()*self.n_nodes(), ndof)]
        id_p = np.arange(self._dim, self.n_fields()*self.n_nodes(), self.n_fields())
        sol = self.solution()
-        # sol[id_v] = data["velocity"][:,:self._dim].reshape([-1])
+        isp2p1 = np.max(self.order_fields) == 2
+        velocity = fdata['velocity'] if isp2p1 else data["velocity"]
+        sol[id_v] = velocity.reshape([-1])
        sol[id_p] = data["pressure"].reshape([-1])
        if self._n_fluids == 2 :
            self.concentration_dg()[:] = cdata["concentration"]
@@ -445,7 +460,7 @@ class FluidProblem :
        return torques

    @timeit
-    def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
+    def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0., solution_old=None) :
        """Solves the fluid equations.

        Keyword arguments:
@@ -460,16 +475,12 @@ class FluidProblem :
        self._lib.fluid_problem_apply_boundary_conditions(self._fp)
        self._lib.fluid_problem_reset_boundary_force(self._fp)
        self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
-        idx = self.get_numbering(self.order_fields)
-        idx_concentration = self.get_numbering(self.order_concentration)
-        idx_porosity = self.get_numbering(self.order_porosity)
-        self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
-        self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
-        self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
+
+        idx = self.local_map()
        # print(time()-timer)
        idx = idx.reshape([self.n_elements(), -1])
        periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
-        solution_old = self.solution().copy()
+        solution_old = self.solution().copy() if solution_old is None else solution_old.copy()
        if self.sys is None :
            constraints = []
            if self.mean_pressure is not None:
@@ -498,7 +509,6 @@ class FluidProblem :
            print("norm",norm)
            if norm > check_residual_norm :
                raise ValueError("wrong derivative or linear system precision")
-        self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
        # print(time()-timer)

    def compute_cfl(self):
@@ -553,20 +563,14 @@ class FluidProblem :
        f.restype = POINTER(typec)
        return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))

-    # def solution(self) :
-        # """Gives access to the fluid field value at the mesh nodes."""
-        # sol = self.solution_flat()[:self.n_fields()*self.n_nodes()].reshape(-1, self.n_fields())
-        # return sol
-        # return self._get_matrix("solution", self.n_nodes(), self.n_fields())
-
    def solution(self) :
        """Gives access to the fluid field value at each degree of freedom"""
-        ndof = np.max(self.get_numbering(self.order_fields))+1
+        ndof = np.max(self.local_map())+1
        return self._get_matrix("solution", ndof, 1).reshape([-1])

    def velocity(self) :
        """Gives access to the fluid velocity value."""
-        ndof = np.max(self.get_numbering(self.order_fields))+1
+        ndof = np.max(self.local_map())+1
        id_v = np.r_[np.repeat(np.arange(self._dim)[None,:],self.n_nodes(),axis=0).reshape(-1)
            + self.n_fields()*np.repeat(np.arange(self.n_nodes()), self._dim),
            np.arange(self.n_fields()*self.n_nodes(), ndof)]
@@ -589,7 +593,7 @@ class FluidProblem :
        """Gives access to the coordinate of the mesh nodes."""
        x = self._get_matrix("coordinates",self.n_nodes(), 3)
        if order != 1:
-            connectivity = self.get_numbering([np.max(order)]).reshape(self.n_elements(),-1)
+            connectivity = self.get_numbering(self.elements(),[np.max(order)]).reshape(self.n_elements(),-1)
            x = np.zeros((connectivity.max()+1, self._dim))
            x[connectivity[:,:self._dim+1]] = self.coordinates()[connectivity[:,:self._dim+1], :self._dim]
            x[connectivity[:,self._dim+1]] = (x[connectivity[:,0]] + x[connectivity[:,1]])/2
@@ -694,12 +698,13 @@ class FluidProblem :
        self._lib.fluid_problem_local_size.restype = c_int
        return self._lib.fluid_problem_local_size(self._fp)

-    def local_map(self, idx):
-        self._lib.fluid_problem_local_map(self._fp, _np2c(idx, dtype=np.int32))
+    def local_map(self):
+        """Gives access to the map of all the degrees of freedom."""
+        return self._get_matrix("idx", self.n_elements(), self.local_size(), c_int)

    @timeit
-    def get_numbering(self, order=[]):
-        mesh_cl = MeshCl(self.elements())
+    def get_numbering(self, elements, order=[]):
+        mesh_cl = MeshCl(elements)
        maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
        ndofs = np.array([lmap.shape[1] for lmap in maps])
        dof_simplices = np.array([1, self._dim+1, 3*self._dim-3], dtype=np.int32)
@@ -709,7 +714,7 @@ class FluidProblem :
            [0, 1, 1]  #p2
        ], dtype=np.int32)
        localsize = np.sum(nloc*dof_simplices[None,:], axis=1)
-        idx_uvp = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
+        idx_uvp = np.full([len(elements), np.sum(ndofs)], -1, dtype=np.int32)
        offsets = np.zeros(ndofs.shape[0], dtype=np.int32)
        for s in range(0, dof_simplices.shape[0]):
            has_simplices = np.where(nloc[order, s] != 0)[0]

--- a/testcases/depot-2d/depot.py
+++ b/testcases/depot-2d/depot.py
@@ -35,48 +35,29 @@ import time
 import shutil
 import random

-def genInitialPosition(filename, r, H, ly, lx, rhop) :
-    """Set all the particles centre positions and create the particles objects to add in the computing structure
-
-    Keyword arguments:
-    filename -- name of the output file
-    r -- max radius of the particles
-    H -- domain height
-    ly - particles area height
-    lx -- particles area width
-    rhop -- particles density
-    """
-    # Particles structure builder
-    p = scontact.ParticleProblem(2)
-    # Load mesh.msh file specifying physical boundaries names
-    p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
-
-    #Definition of the points where the particles are located
-    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2*r)
-    y = np.arange(H/2-r, H/2-ly+r, -2*r)
-    x, y = np.meshgrid(x, y)
-    x = x.flat
-    y = y.flat
-    # Add a grain at each centre position
-    for i in range(len(x)) :
-        p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
-    p.write_vtk(filename,0,0)
+
+def gen_hexagonal_distribution(r,w,h,x0=np.array([0,0])):
+    step = 2*r*np.sqrt(3)
+    grid_0 = np.stack(np.meshgrid(np.arange(0, w, step), np.arange(0, h, 2*r)), axis=2).reshape(-1,2)
+    grid_1 = np.stack(np.meshgrid(np.arange(step/2, w, step), np.arange(r, h, 2*r)), axis=2).reshape(-1,2)
+    return np.r_[grid_0, grid_1] + x0[None,:]
+

 # Define output directory
-outputdir = "output"
+outputdir = "output_ratio5_subloop11"
 if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

 # Physical parameters
 g = np.array([0,-9.81])                         # gravity
-r = 1.5e-3                                        # particles radius
-rhop = 1500                                     # particles density
+r = 2e-3                                        # particles radius
+rhop = 5000                                     # particles density
 rho = 1000                                      # fluid density
-nu = 1e-6                                       # kinematic viscosity
+nu = 2e-6                                       # kinematic viscosity

 # Numerical parameters
-outf = 5                                        # number of iterations between output files
-dt = 1e-3                                     # time step
+outf = 1                                         # number of iterations between output files
+dt = 1e-3                                       # time step
 tEnd = 100                                      # final time

 # Geometrical parameters
@@ -88,11 +69,13 @@ H = 0.6                                         # domain height
 # PARTICLE PROBLEM
 #
 # Initialise particles
-genInitialPosition(outputdir, r, H, ly, lx, rhop)
 p = scontact.ParticleProblem(2)
-p.read_vtk(outputdir,0)
+x = gen_hexagonal_distribution(1*r, lx-2.1*r, ly, np.array([-lx/2+1.25*r, H/2-ly]))
+print(x.shape)
+p.add_particles(x, np.full(x.shape[0], r), np.full(x.shape[0], rhop*np.pi*r**2))
+p.load_msh_boundaries("mesh.msh", ["Bottom", "Lateral", "Top"])

-print("r = %g, m = %g\n" %  (p.r()[0], p.mass()[0]))
+print("r = %g, m = %g, n_particles = %d\n" %  (p.r()[0], p.mass()[0], p.n_particles()))
 print("RHOP = %g" % rhop)

 # Initial time and iteration
@@ -102,7 +85,7 @@ ii        =  0
 #
 # FLUID PROBLEM
 #
-fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu", stability=False)
+fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], order=[2,2,1], solver="petsc4py", solver_options="-pc_type lu")#, stability=True)
 # Set the mesh geometry for the fluid computation
 fluid.load_msh("mesh.msh")
 fluid.set_wall_boundary("Bottom")
@@ -118,9 +101,48 @@ tic = time.time()
 #
 # COMPUTATION LOOP
 #
+sub_iter_max = 20
+sub_iter = 0
+tol = 5e-2
+# stab_param = [1e-2]*sub_iter_max
+# stab_param = [10e-7, 5e-7, 1e-7, 1e-7]
+stab_param = 10*dt
 while t < tEnd :
-    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass(), check_residual_norm=1e-3)

+    initial_sol = fluid.solution().copy()
+    p.save_state()
+    sub_iter = 0
+    error = tol+1
+    pv = fluid.particle_velocity().copy()
+
+    while sub_iter < sub_iter_max :#and error > tol:
+        sol_star_old = fluid.solution().copy()
+        p.restore_state()
+        fluid.particle_velocity()[:] = pv
+        fext_old = fluid.compute_node_force(dt)+g*p.mass()
+        stab = stab_param/(5**min(sub_iter, 10))
+        fluid.implicit_euler(dt, check_residual_norm=1, solution_old=initial_sol, stab_param=stab)
+        fext = fluid.compute_node_force(dt)+g*p.mass()
+        time_integration._advance_particles(p, (fext+fext_old)/2, dt, min_nsub=5, contact_tol=1e-8)
+        pv = p.velocity().copy()
+        
+        error_max = np.max(np.abs(sol_star_old - fluid.solution()))
+        error_inc = np.linalg.norm(initial_sol - fluid.solution())
+        fmax = np.sqrt(np.max(np.sum(fext**2, axis=1)))
+        error_fext = np.linalg.norm(fext_old - fext)/fmax
+        sub_iter += 1
+        print("sub_iter : ", sub_iter, "/", sub_iter_max, "error u* max:", error_max, "u_{n+1} - u_n :", error_inc, "error f_ext : ", error_fext, "fmax : ", fmax, "stab :", stab)
+        error = error_fext
+