Commit 8c77ef87 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

really merge master

parent 2b84fd8a
Pipeline #4707 passed with stage
in 22 seconds
......@@ -28,7 +28,6 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"
#define D DIMENSION
......@@ -469,8 +468,8 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
double normu = {0.};
const uint32_t *el = &mesh->elements[iel*N_N];
double a = 0;
double amax = 0;
double amin = 0;
double amax = 0.5;
double amin = 0.5;
double c = 0;
for (int i=0; i< N_N; ++i) {
double normup = 0;
......
# 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/>.
import glob
import os
from paraview.simple import *
for sname,sid in GetSources() :
s = FindSource(sname)
if not hasattr(s,'FileName') :
continue
l = s.FileName.GetData()
for ext in [".vtu",".vtp",".vtm"] :
if l and l[0][-4:]==ext:
p = l[0].rfind("_")
basename = l[0][:p]
idname = l[0][p+1:-4]
idnamelength = len(idname)
fmt = basename+"_%0"+str(idnamelength)+"i"+ext
t0 = os.path.getmtime(l[0])
allfiles = []
i = 1
while True :
fname = fmt%i
try :
mt = os.path.getmtime(fname)
if mt < t0 :
break
except :
break
allfiles += [fname]
i += 1
if allfiles :
s.FileName.SetData(allfiles)
......@@ -29,15 +29,15 @@
"""
from ctypes import *
import numpy
import numpy as np
import signal
import os
import sys
from . import VTK
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = CDLL(os.path.join(dir_path,"libmbfluid2.so"))
lib3 = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
signal.signal(signal.SIGINT, signal.SIG_DFL)
......@@ -48,8 +48,8 @@ class _Bnd :
self._dim = dim
def apply(self, n, xp, vp) :
nv = len(self._b)
x = numpy.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,-1])
v = numpy.frombuffer(cast(vp, POINTER(int(n)*nv*c_double)).contents).reshape([n,nv])
x = np.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,-1])
v = np.frombuffer(cast(vp, POINTER(int(n)*nv*c_double)).contents).reshape([n,nv])
for i in range(nv):
if callable(self._b[i]) :
v[:,i] = self._b[i](x)
......@@ -83,7 +83,7 @@ class FluidProblem :
self.strong_cb_ref = []
self.weak_cb_ref = []
def np2c(a) :
tmp = numpy.require(a,"float","C")
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
......@@ -94,7 +94,7 @@ class FluidProblem :
else :
raise ValueError("dim should be 2 or 3.")
self._lib.fluid_problem_new.restype = c_void_p
n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
self._dim = dim
self._fp = c_void_p(self._lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), c_double(coeff_stab), petsc_solver_type.encode()))
......@@ -147,7 +147,7 @@ class FluidProblem :
raise NameError("Pressure or Velocity (but not both) should be specified at open boundaries")
if velocity is not None :
n_value = self._dim+self._n_fluids-1 # u, v, w, a
cb_or_value = [*velocity,concentration] if self._n_fluids == 2 else velocity
cb_or_value = velocity+[concentration] if self._n_fluids == 2 else velocity
bnd_type = "velocity"
else :
n_value = self._n_fluids # p, a
......@@ -192,9 +192,12 @@ class FluidProblem :
t -- computation time
it -- computation iteration
"""
v = self.solution()[:,:self._dim]
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
data = [
("pressure",self.solution()[:,[self._dim]]),
("velocity",self.solution()[:,:self._dim]),
("velocity",v),
("porosity",self.porosity()),
("old_porosity",self.old_porosity())
]
......@@ -205,7 +208,10 @@ class FluidProblem :
name,dim,pid,nodes,bnd = self._physical(i)
field_data.append((b"Physical %i %i %s" % (dim,pid,name),nodes[:,None]))
field_data.append((b"Boundary %i %i %s" % (dim,pid,name),bnd))
VTK.write(output_dir+"/fluid",it,t,self.elements(),self.coordinates(),data,field_data)
connectivity = self.elements()
types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0]))
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data)
def import_vtk(self, filename) :
"""Read output file to reload computation.
......@@ -214,17 +220,18 @@ class FluidProblem :
filename -- name of the file to read
"""
def np2c(a,dtype) :
tmp = numpy.require(a,dtype,"C")
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
x,el,data,fdata = VTK.read(filename)
x,cells,data,cdata,fdata = VTK.read(filename)
el = cells["connectivity"].reshape([-1,self._dim+1])
self._lib.fluid_problem_private_set_elements(self._fp,
c_int(x.shape[0]),np2c(x,numpy.float64),
c_int(el.shape[0]),np2c(el,numpy.int32))
c_int(x.shape[0]),np2c(x,np.float64),
c_int(el.shape[0]),np2c(el,np.int32))
phys = {}
bnds = {}
for n,d in fdata :
for n,d in fdata.items() :
k,dim,pid,name = n.split(None,4)
if k == "Physical" :
phys[(dim,pid)] = (name,d)
......@@ -234,11 +241,10 @@ class FluidProblem :
bnd = bnds[(dim,pid)]
self._lib.fluid_problem_private_add_physical(self._fp,
c_int(int(pid)), c_int(int(dim)), c_char_p(name.encode("utf-8")),
c_int(nodes.shape[0]), np2c(nodes,numpy.int32),
c_int(bnd.shape[0]), np2c(bnd,numpy.int32))
c_int(nodes.shape[0]), np2c(nodes,np.int32),
c_int(bnd.shape[0]), np2c(bnd,np.int32))
sol = self.solution()
data = dict(data)
sol[:,:self._dim] = data["velocity"]
sol[:,:self._dim] = data["velocity"][:,:self._dim]
sol[:,[self._dim]] = data["pressure"]
if self._n_fluids == 2 :
sol[:,[self._dim+1]] = data["concentration"]
......@@ -252,7 +258,7 @@ class FluidProblem :
Keyword arguments:
dt -- computation time step
"""
forces = numpy.ndarray([self.n_particles,self._dim],"d",order="C")
forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
self._lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
return forces
......@@ -278,7 +284,7 @@ class FluidProblem :
reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
def np2c(a) :
tmp = numpy.require(a,"float","C")
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
......@@ -288,7 +294,7 @@ class FluidProblem :
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
f = getattr(self._lib,"fluid_problem_"+f_name)
f.restype = POINTER(typec)
return numpy.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
def solution(self) :
"""Give access to the fluid field value at the mesh nodes."""
......@@ -334,7 +340,7 @@ class FluidProblem :
ptr = f(self._fp)
size = fs(self._fp)
buf = (size*c_int).from_address(ptr)
return numpy.ctypeslib.as_array(buf)
return np.ctypeslib.as_array(buf)
def _n_physicals(self):
f = self._lib.fluid_problem_private_n_physical
......@@ -351,7 +357,7 @@ class FluidProblem :
phys_n_boundaries = c_int()
phys_boundaries = POINTER(c_int)()
f(self._fp,c_int(ibnd),byref(phys_name),byref(phys_dim),byref(phys_id),byref(phys_n_nodes),byref(phys_nodes),byref(phys_n_boundaries),byref(phys_boundaries))
nodes = numpy.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else numpy.ndarray([0],numpy.int32)
bnd = numpy.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else numpy.ndarray([2,0],numpy.int32)
nodes = np.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else np.ndarray([0],np.int32)
bnd = np.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else np.ndarray([2,0],np.int32)
return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd
......@@ -30,8 +30,9 @@
from __future__ import division
from . import mesh
from . import VTK
import shutil
import numpy
import numpy as np
import os
import sys
from ctypes import *
......@@ -40,11 +41,17 @@ import signal
signal.signal(signal.SIGINT, signal.SIG_DFL)
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = CDLL(os.path.join(dir_path,"libscontact2.so"))
lib3 = CDLL(os.path.join(dir_path,"libscontact3.so"))
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
is_64bits = sys.maxsize > 2**32
def _np2c(a,dtype=np.float64) :
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
class ParticleProblem :
"""Create the numerical structure containing all the physical particles that appear in the problem"""
......@@ -54,7 +61,7 @@ class ParticleProblem :
ptr = f(self._p)
size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//8
buf = (size*c_double).from_address(ptr)
return numpy.ctypeslib.as_array(buf).reshape([-1,ncol])
return np.ctypeslib.as_array(buf).reshape([-1,ncol])
def _get_idx(self, fName) :
f = getattr(self._lib,"particleProblem"+fName)
......@@ -62,7 +69,7 @@ class ParticleProblem :
ptr = f(self._p)
size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//4
buf = (size*c_int).from_address(ptr)
return numpy.ctypeslib.as_array(buf)
return np.ctypeslib.as_array(buf)
def __init__(self, dim) :
"""Build the particles contanier structure.
......@@ -96,17 +103,17 @@ class ParticleProblem :
def volume(self):
"""Return the list of particle volumes."""
if self._dim == 2 :
return numpy.pi * self._get_matrix("Particle", 2)[:, [0]]**2
return np.pi * self._get_matrix("Particle", 2)[:, [0]]**2
else :
return 4./3.*numpy.pi * self._get_matrix("Particle", 2)[:, [0]]**3
return 4./3.*np.pi * self._get_matrix("Particle", 2)[:, [0]]**3
def r(self) :
"""Return the list of particle radii."""
return self._get_matrix("Particle", 2)[:, [0]]
return self._get_matrix("Particle", 2)[:, 0][:,None]
def mass(self):
"""Return the list of particle masses."""
return self._get_matrix("Particle", 2)[:, [1]]
return self._get_matrix("Particle", 2)[:, 1][:,None]
def position(self):
"""Return the list of particle centre position vectors."""
......@@ -116,24 +123,40 @@ class ParticleProblem :
"""Return the list of particle velocity vectors."""
return self._get_matrix("Velocity", self._dim)
def omega(self):
"""Return the list of particle angular velocity."""
return self._get_matrix("Omega", 1)
def disk_tag(self):
return self._get_idx("DiskTag")
def segments(self):
return self._get_matrix("Segment",6)
return self._get_matrix("Segment",3*self._dim+1)
def segment_tag(self):
return self._get_idx("SegmentTag")
def triangles(self):
return self._get_matrix("Triangle",12)
if self._dim == 3 :
return self._get_matrix("Triangle",12)
else :
return np.ndarray([0,12])
def triangle_tag(self):
if self._dim == 3 :
return self._get_idx("TriangleTag")
else :
return np.array([],np.int32)
def set_boundary_velocity(self, tag, v) :
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
disks = self._get_matrix("Disk",5)
disks = self._get_matrix("Disk",2*self._dim+1)
diskTag = self._get_idx("DiskTag")
disks[diskTag == tag, 2:4] = v
segments = self._get_matrix("Segment",6)
disks[diskTag == tag, self._dim:2*self._dim] = v
segments = self._get_matrix("Segment",3*self._dim)
segmentTag = self._get_idx("SegmentTag")
segments[segmentTag == tag, 4:6] = v
segments[segmentTag == tag, 2*self._dim:3*self._dim] = v
if self._dim == 3 :
triangles = self._get_matrix("Triangle",12)
triangleTag = self._get_idx("TriangleTag")
......@@ -158,15 +181,87 @@ class ParticleProblem :
tol -- optional argument defining the interpenetration tolerance to stop the NLGS iterations of the NSCD
"""
self._get_matrix("ExternalForces",self._dim)[:] = forces
self._lib.particleProblemIterate(self._p, c_double(numpy.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
self._lib.particleProblemIterate(self._p, c_double(np.min(self.r())), c_double(dt), c_double(tol), c_int(-1))
def write_vtk(self, odir, i, t) :
"""Write output files for post-visualization."""
self._lib.particleProblemWriteVtk(self._p, odir.encode(), i)
v = self.velocity()
omega = self.omega()
x = self.position()
tag = self._get_idx("ParticleTag").reshape(-1,1)
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
x = np.insert(x,self._dim,0,axis=1)
data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Tag",tag)]
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp")
xdisk = self._get_matrix("Disk",2*self._dim+2)[:,:self._dim]
nd = xdisk.shape[0]
xsegment = self.segments()[:,:self._dim*2].reshape([-1,self._dim])
ns = xsegment.shape[0]/2
xtriangles = self.triangles()[:,:self._dim*3].reshape([-1,self._dim])
nt = xtriangles.shape[0]/3
x = np.vstack([xdisk,xsegment,xtriangles])
if self._dim == 2 :
x = np.insert(x,self._dim,0,axis=1)
connectivity = np.arange(0,x.shape[0],dtype=np.int32)
types = np.repeat(np.array([1,3,5],dtype=np.int32),[nd,ns,nt])
offsets = np.cumsum(np.repeat([1,2,3],[nd,ns,nt]),dtype=np.int32)
tags = np.array(np.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]),dtype=np.float64)
VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1]))])
self._lib.particleProblemContactN.restype = c_size_t
fdata = []
for (ctype,name) in ((0,b"particle_particle"),(1,b"particle_disk"),(2,b"particle_segment"),(3,b"particle_triangle")) :
n = self._lib.particleProblemContactN(self._p,c_int(ctype))
v = np.ndarray((n,2),dtype=np.float64,order="C")
o = np.ndarray((n,2),dtype=np.uint64,order="C")
self._lib.particleProblemContact(self._p,c_int(ctype),c_void_p(o.ctypes.data),c_void_p(v.ctypes.data))
fdata.append((name,v[:,[0]]))
fdata.append((name+b"_t",v[:,[1]]))
fdata.append((name+b"_idx",o))
x = np.array([[0.,0.,0.]])
VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
def read_vtk(self,dirname,i) :
"""Read an existing output file to set the particle data."""
self._lib.particleProblemReadVtk(self._p, dirname.encode(),i)
x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
print(d["Omega"])
self._lib.particleProblemClearParticles(self._p)
self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
_np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(d["Tag"]))
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
self._lib.particleProblemClearBoundaries(self._p)
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
tags = cdata["Tag"]
for idx in np.where(cells["types"] == 1)[0] :
x0 = x[el[offsets[idx]],:self._dim]
t = int(tags[idx,0])
self._lib.particleProblemAddBoundaryDiskTagId(self._p, self._coord_type(*x0), c_double(0.), c_int(t))
for idx in np.where(cells["types"] == 3)[0] :
x0 = x[el[offsets[idx]],:self._dim]
x1 = x[el[offsets[idx]+1],:self._dim]
t = int(tags[idx,0])
self._lib.particleProblemAddBoundarySegmentTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), c_int(t))
for idx in np.where(cells["types"] == 5)[0] :
x0 = x[el[offsets[idx]],:self._dim]
x1 = x[el[offsets[idx]+1],:self._dim]
x2 = x[el[offsets[idx]+2],:self._dim]
t = int(tags[idx,0])
self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(t))
_,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
v = np.vstack([fdata[k] for k in ks])
vt = np.vstack([fdata[k+"_t"] for k in ks])
v = np.hstack([v,vt])
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32))
def add_boundary_disk(self, x0, r, tag) :
self._lib.particleProblemAddBoundaryDisk.restype = c_size_t
......
This diff is collapsed.
......@@ -24,7 +24,6 @@
#ifndef _SCONTACT_H_
#define _SCONTACT_H_
#include <stdlib.h>
#include "quadtree.h"
typedef struct _ParticleProblem ParticleProblem;
......@@ -54,8 +53,6 @@ double *particleProblemExternalForces(ParticleProblem *p);
void particleProblemSetUseQueue(ParticleProblem *p, int use_queue);
void particleProblemSetFrictionRelaxation(ParticleProblem *p, double friction_relaxation);
unsigned long int particleProblemNParticle(ParticleProblem *p);
int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id);
int particleProblemReadVtk(ParticleProblem *p, const char *dirname, int iter);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
int *particleProblemSegmentTag(ParticleProblem *p);
......
......@@ -74,15 +74,10 @@ ii = 0
fluid = fluid.FluidProblem(2,g,[nuf*rhof],[rhof])
# Set the mesh geometry for the fluid computation.
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Right",0,0)
fluid.set_strong_boundary("Left",0,0)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Left","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Right","Wall")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set locations of the grains in the mesh and compute the porosity in each computation cell.
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
# Write output files for post-visualization.
......
......@@ -116,15 +116,10 @@ ii = 0
fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Right",0,0)
fluid.set_strong_boundary("Left",0,0)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Left","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Right","Wall")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
#Access to fluid fields and node coordiantes
......
......@@ -78,15 +78,10 @@ ii = 0
fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof],coeff_stab=0.001)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",2,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Right",0,0)
fluid.set_strong_boundary("Left",0,0)
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Left","Wall")
fluid.set_weak_boundary("Top","Wall")
fluid.set_weak_boundary("Right","Wall")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Left")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_wall_boundary("Right")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
# Access to fluid fields and node coordiantes
......
......@@ -26,7 +26,7 @@
# This example shows how to set boundary conditions as a function of some parameters.
# A boolean parameter gives the possibility to use lmgc90 to solve contacts instead of scontact.
use_lmgc90 = 0
use_lmgc90 = True
from migflow import fluid
from migflow import scontact
......@@ -52,6 +52,7 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
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, y = np.meshgrid(x, x)
......@@ -69,12 +70,11 @@ def genInitialPosition(filename, r, rout, rin, rhop) :
for i in range(x.shape[0]) :
if y[i]<0:
rhop1 = random.choice([rhop*.9,1.1*rhop,rhop])
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop1,0)
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) :
......@@ -100,7 +100,6 @@ r = 397e-6/2 # grains radius
# PARTICLE PROBLEM
#
# Object particles creation
genInitialPosition(outputdir, r, rout, rin, rhop)
if use_lmgc90 :
friction = 0.1 # friction coefficient
......@@ -125,8 +124,8 @@ fluid.set_strong_boundary("Outer",1,lambda x : x[:, 0])
fluid.set_strong_boundary("Inner",0,0)
fluid.set_strong_boundary("Inner",1,0)
fluid.set_strong_boundary("PtFix",2,0)
fluid.set_weak_boundary("Inner","Wall")