Skip to content
GitLab
Commit 6b0ab1f9 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

numpy -> np

parent 21998c98
Pipeline #4633 passed with stage
in 21 seconds
Hide whitespace changes
Inline Side-by-side
python/fluid.py
View file @ 6b0ab1f9
......@@ -29,16 +29,15 @@
"""
from ctypes import *
import numpy
from numpy import ctypeslib
import numpy as np
import signal
import os
import sys
from . import VTK
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = ctypeslib.load_library("libmbfluid3",dir_path)
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
signal.signal(signal.SIGINT, signal.SIG_DFL)
......@@ -49,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)
......@@ -84,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
......@@ -95,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()))
......@@ -195,7 +194,7 @@ class FluidProblem :
"""
v = self.solution()[:,:self._dim]
if self._dim == 2 :
v = numpy.insert(v,self._dim,0,axis=1)
v = np.insert(v,self._dim,0,axis=1)
data = [
("pressure",self.solution()[:,[self._dim]]),
("velocity",v),
......@@ -210,8 +209,8 @@ class FluidProblem :
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))
connectivity = self.elements()
types = numpy.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
offsets = numpy.cumsum(numpy.repeat([self._dim+1],connectivity.shape[0]))
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) :
......@@ -221,15 +220,15 @@ 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,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.items() :
......@@ -242,8 +241,8 @@ 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()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
sol[:,[self._dim]] = data["pressure"]
......@@ -259,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
......@@ -285,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
......@@ -295,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."""
......@@ -341,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
......@@ -358,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
python/scontact.py
View file @ 6b0ab1f9
......@@ -32,23 +32,22 @@ 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 *
import signal
from numpy import ctypeslib
signal.signal(signal.SIGINT, signal.SIG_DFL)
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = ctypeslib.load_library("libscontact2",dir_path)
lib3 = ctypeslib.load_library("libscontact3",dir_path)
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=numpy.float64) :
tmp = numpy.require(a,dtype,"C")
def _np2c(a,dtype=np.float64) :
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
......@@ -62,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)
......@@ -70,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.
......@@ -104,9 +103,9 @@ 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."""
......@@ -141,13 +140,13 @@ class ParticleProblem :
if self._dim == 3 :
return self._get_matrix("Triangle",12)
else :
return numpy.ndarray([0,12])
return np.ndarray([0,12])
def triangle_tag(self):
if self._dim == 3 :
return self._get_idx("TriangleTag")
else :
return numpy.array([],numpy.int32)
return np.array([],np.int32)
def set_boundary_velocity(self, tag, v) :
self._lib.particleProblemGetTagId.restype = c_size_t
......@@ -172,7 +171,7 @@ 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."""
......@@ -180,8 +179,8 @@ class ParticleProblem :
v = self.velocity()
x = self.position()
if self._dim == 2 :
v = numpy.insert(v,self._dim,0,axis=1)
x = numpy.insert(x,self._dim,0,axis=1)
v = np.insert(v,self._dim,0,axis=1)
x = np.insert(x,self._dim,0,axis=1)
data = [("Mass",self.mass()), ("Radius",self.radius()),("Velocity",v)]
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp")
......@@ -191,25 +190,25 @@ class ParticleProblem :
ns = xsegment.shape[0]/2
xtriangles = self.triangles()[:,:self._dim*3].reshape([-1,self._dim])
nt = xtriangles.shape[0]/3
x = numpy.vstack([xdisk,xsegment,xtriangles])
x = np.vstack([xdisk,xsegment,xtriangles])
if self._dim == 2 :
x = numpy.insert(x,self._dim,0,axis=1)
connectivity = numpy.arange(0,x.shape[0],dtype=numpy.int32)
types = numpy.repeat(numpy.array([1,3,5],dtype=numpy.int32),[nd,ns,nt])
offsets = numpy.cumsum(numpy.repeat([1,2,3],[nd,ns,nt]),dtype=numpy.int32)
tags = numpy.array(numpy.hstack([self.disk_tag(),self.segment_tag(),self.triangle_tag()]),dtype=numpy.float64)
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 = numpy.ndarray((n,1),dtype=numpy.float64,order="C")
o = numpy.ndarray((n,2),dtype=numpy.uint64,order="C")
v = np.ndarray((n,1),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))
fdata.append((name+b"_idx",o))
x = numpy.array([[0.,0.,0.]])
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",["pyparticles_%05d.vtp"%i,"pybnd_%05d.vtu"%i,"pycontact_%05d.vtp"%i],i,t)
......@@ -223,19 +222,19 @@ class ParticleProblem :
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
self._lib.particleProblemClearBoundaries(self._p)
offsets = numpy.hstack([[0],cells["offsets"]])
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
tags = cdata["Tag"]
for idx in numpy.where(cells["types"] == 1)[0] :
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 numpy.where(cells["types"] == 3)[0] :
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 numpy.where(cells["types"] == 5)[0] :
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]
......@@ -243,10 +242,10 @@ class ParticleProblem :
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 = numpy.vstack([fdata[k] for k in ks])
o = numpy.vstack([fdata[k+"_idx"] for k in ks])
t = numpy.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,numpy.uint64),_np2c(v),_np2c(t,numpy.int32))
v = np.vstack([fdata[k] for k in ks])
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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment