VTK

fluid/fluid_problem.h

@@ -138,7 +138,6 @@ 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 *fluid_problem_parent_tag(const FluidProblem *p);
 
 
 void fluid_problem_set_stab_param(FluidProblem *p, double  stab_param);

python/fluid.py

@@ -345,10 +345,11 @@ class FluidProblem :
             ("velocity",v),
             ("porosity",self.porosity()),
             ("old_porosity",self.old_porosity()),
-            ("grad",da)
+            ("grad",da),
+            ("parent_node_id", self.parent_nodes())
             ]
-        cell_data.append(("curvature",self.curvature()))
         if self._n_fluids == 2 :
+            cell_data.append(("curvature",self.curvature()))
             cell_data.append(("concentration",self.concentration_dg()))
             cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
         field_data = [(b"Boundary %s"%(name), nodes) for name,nodes in self._mesh_boundaries().items()]
@@ -371,11 +372,31 @@ class FluidProblem :
         bnds = np.vstack(list(mesh_boundaries.values()))
         bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
         bnd_tags = np.require(bnd_tags,np.int32,"C")
-        self._lib.fluid_problem_set_elements(self._fp,
+        # self._lib.fluid_problem_set_elements(self._fp,
+        #         c_int(x.shape[0]),_np2c(x,np.float64),
+        #         c_int(el.shape[0]),_np2c(el,np.int32),
+        #         c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
+        #         )
+
+        # lib.mesh_new_from_elements.restype = c_void_p
+        # return c_void_p(lib.mesh_new_from_elements(
+                # c_int(x.shape[0]),_np2c(x,np.float64),
+                # c_int(el.shape[0]),_np2c(el,np.int32),
+                # c_int(bnd.shape[0]),_np2c(bnd,np.int32),
+                # _np2c(btag,np.int32),c_int(len(cbname)),cbname,
+                # _np2c(periodic_parent,np.int32)))
+        self.parent_nodes()[:] = np.require(data["parent_node_id"],np.int32,"C")
+        self._lib.mesh_new_from_elements.restype = c_void_p
+        _mesh = c_void_p(self._lib.mesh_new_from_elements(
                 c_int(x.shape[0]),_np2c(x,np.float64),
                 c_int(el.shape[0]),_np2c(el,np.int32),
-                c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
-                )
+                c_int(bnds.shape[0]),_np2c(bnds,np.int32),
+                _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
+                _np2c(data["parent_node_id"],np.int32)))
+
+        self._lib.fluid_problem_set_mesh.restype = None
+        self._lib.fluid_problem_set_mesh(self._fp, _mesh)
+
         sol = self.solution()
         sol[:,:self._dim] = data["velocity"][:,:self._dim]
         sol[:,[self._dim]] = data["pressure"]
@@ -523,6 +544,10 @@ class FluidProblem :
         """Gives access to the coordinate of the mesh nodes."""
         return self._get_matrix("coordinates",self.n_nodes(), 3)
 
+    def parent_nodes(self):
+        """Gives access to the parent nodes of each node."""
+        return self._get_matrix("periodic_mapping", self.n_nodes(),1, typec=c_int32)
+
     def n_fluids(self) :
         """Returns the number of fluids."""
         return self._n_fluids
@@ -534,17 +559,17 @@ class FluidProblem :
     def n_elements(self):
         """Returns the number of mesh nodes."""
         self._lib.fluid_problem_n_elements.restype = c_int
-        return self._lib.fluid_problem_n_elements(self._fp);
+        return self._lib.fluid_problem_n_elements(self._fp)
 
     def n_fields(self):
         """Returns the number of fluid fields."""
         self._lib.fluid_problem_n_fields.restype = c_int
-        return self._lib.fluid_problem_n_fields(self._fp);
+        return self._lib.fluid_problem_n_fields(self._fp)
 
     def n_nodes(self):
         """Returns the number of mesh nodes."""
         self._lib.fluid_problem_n_nodes.restype = c_int
-        return self._lib.fluid_problem_n_nodes(self._fp);
+        return self._lib.fluid_problem_n_nodes(self._fp)
 
     def porosity(self):
         """Returns the porosity at nodes"""

testcases/periodicTest/2D/poiseuille.py

@@ -51,7 +51,7 @@ outputdir = "outputPoiseuille"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
-subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","2"])
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","0.5"])
 t = 0
 ii = 0
 
@@ -64,40 +64,29 @@ mu = nu*rho                                     # dynamic viscosity
 tEnd = 1000                                     # final time
 
 #numerical parameters
-lcmin = .1                                  # mesh size
-dt = 10                                       # time step
-alpha = 1e-4                                    # stabilization coefficient
-epsilon = alpha*lcmin**2 /nu                    # stabilization parametre
-print('epsilon',epsilon)
+dt = 10                                         # time step
+outf = 10                                       # number of iterations between output files
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
 
-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")
 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)
+# # 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
+t = ii*dt
 
 #set initial_condition
+fluid.export_vtk(outputdir,t,ii)
 
-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)

testcases/periodicTest/2D/testVTK.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/>.
+
+#!/usr/bin/env python
+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
+import unittest
+
+def get_nodes(fluid,tag):
+    edges = get_boundary(fluid,tag)
+    nodes = edges.reshape(np.size(edges))
+    nodes = np.unique(nodes)
+    nodes = sorted(nodes,key=lambda node: fluid.coordinates()[node,0])
+    return nodes
+
+# Return the edges which belong to the boundary identified by tag
+def get_boundary(fluid,tag):
+    bnds = fluid._mesh_boundaries()
+    edges = bnds.get(tag.encode())
+    return(edges)
+
+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 = "outputPoiseuilleVTK"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+# subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","0.5"])
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  np.array([0.01,0])                                      # gravity
+rho = 1000                                      # fluid density
+nu = 1e-3                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 1000                                     # final time
+
+#numerical parameters
+dt = 10                                         # time step
+outf = 10                                       # number of iterations between output files
+
+# shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#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") # a retirer
+fluid.import_vtk("outputPoiseuille/fluid_%05d.vtu"%ii)
+print("Import done !")
+
+
+
+ii = 0
+t = ii*dt
+
+#set initial_condition
+
+fluid.export_vtk(outputdir,t,ii)
+tic = time.time()
+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))
+
+
+nodesRight = get_nodes(fluid, "Right")
+nodesLeft = get_nodes(fluid, "Left")
+s = fluid.solution()
+error = np.abs(np.max(s[nodesRight,:],axis = 0) - np.max(s[nodesLeft,:], axis = 0))
+eps = 1e-6
+print(f"Periodic success ? {error} <? {eps} => {np.all(error < eps)}")