Commit 160a62c2 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

Merge branch 'binary' into 'master'

Binary

See merge request !3
parents f02dcfab 6b0ab1f9
Pipeline #4635 passed with stage
in 21 seconds
......@@ -28,7 +28,6 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"
#define D DIMENSION
......
# 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)
import zlib
import struct
import os
import numpy as np
import re
from base64 import b64decode,b64encode
from xml.etree import ElementTree as ET
import os.path
import os
def _typename(a) :
if a.dtype == np.int32 :
return b"Int32"
elif a.dtype == np.int64 :
return b"Int64"
elif a.dtype == np.float64 :
return b"Float64"
elif a.dtype == np.float32 :
return b"Float32"
else :
raise(NameError("unown array type : " + str(a.dtype)))
def _write_array_ascii(f,a,attr) :
f.write(b'<DataArray %s type="%s" format="ascii">\n'%(attr,_typename(a)))
f.write((" ".join(map(str,a.flat))).encode("utf-8"))
f.write(b'\n</DataArray>\n')
def _write_array(f,appended,a,ds,attr) :
f.write(b'<DataArray %s type="%s" format="appended" offset="%i"/>\n'%(attr,_typename(a),len(appended)-1))
def _write_array(f,anc,attr) :
a = np.ascontiguousarray(anc)
tname = {
np.dtype('<i4'):b"Int32",
np.dtype('<i8'):b"Int64",
np.dtype('<f4'):b"Float32",
np.dtype('<f8'):b"Float64",
np.dtype('uint64'):b"UInt64"
}[a.dtype]
f.write(b'<DataArray %s type="%s" format="binary">\n'%(attr,tname))
data = zlib.compress(a,2)
# return appended + struct.pack("iiii",1,a.size*ds,0,len(data)) + data
return appended + b64encode(struct.pack("iiii",1,a.size*ds,0,len(data))) + b64encode(data)
f.write(b64encode(struct.pack("iiii",1,a.nbytes,0,len(data))) + b64encode(data))
f.write(b"\n</DataArray>\n")
def _read_array(a) :
raw = a.text.strip()
typename = a.attrib["type"]
dtype = {"Float64":np.float64,"Float32":np.float32,"Int64":np.int64,"Int32":np.int32,"UInt64":np.uint64}[typename]
_,n,_,s = struct.unpack("iiii",b64decode(raw[:24]))
data = zlib.decompress(b64decode(raw[24:24-(-s//3)*4]))
nc,nt = -1,-1
if "NumberOfComponents" in a.attrib :
nc = int(a.attrib["NumberOfComponents"])
if "NumberOfTuples" in a.attrib :
nt = int(a.attrib["NumberOfTuples"])
v = np.frombuffer(data,dtype)
if nc != -1 or nt != -1 :
v = v.reshape(nt,nc)
return v
def _write_index(idxname,filename,i,t,reset_index) :
if reset_index is None : reset_index = (i==0 or not os.path.isfile(idxname))
if reset_index :
f = open(idxname,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">\n')
f.write(b'<Collection>\n')
else:
f = open(idxname,"rb+")
f.seek(-25,os.SEEK_END)
f.write(b'<DataSet timestep="%.16g" group="" part="0" file="%s"/>\n'%(t,filename.encode()))
f.write(b'</Collection>\n')
f.write(b'</VTKFile>\n')
f.close()
def write(basename,i,t,elements,x,data,field_data=None,cell_data=None) :
appended = b"_"
filename = "%s_%05d.vtu"%(basename,i)
def write(basename,i,t,elements,x,data=None,field_data=None,cell_data=None,reset_index=None,vtktype="vtu") :
filename = "%s_%05d.%s"%(basename,i,vtktype)
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="UnstructuredGrid" format="0.1" compressor="vtkZLibDataCompressor">\n')
f.write(b'<UnstructuredGrid>\n')
nel = elements.shape[0]
nptbyel = elements.shape[1]
vtkdescr = b"UnstructuredGrid" if vtktype=="vtu" else b"PolyData"
f.write(b'<VTKFile type="%s" format="0.1" compressor="vtkZLibDataCompressor">\n' % vtkdescr)
f.write(b'<%s>\n'%vtkdescr)
nel = elements[0].shape[0] if elements is not None else 0
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
appended = _write_array(f,appended,x,8,b'NumberOfComponents="3"')
_write_array(f,x,b'NumberOfComponents="3"')
f.write(b'</Points>\n')
f.write(b'<Cells>\n')
appended = _write_array(f,appended,elements,4,b'Name="connectivity"')
offsets = np.array(range(nptbyel,(nel+1)*nptbyel,nptbyel),np.int32)
appended = _write_array(f,appended,offsets,4,b'Name="offsets"')
types = np.ones(nel,np.int32)*(5 if elements.shape[1] == 3 else 13)
appended = _write_array(f,appended,types,4,b'Name="types"')
f.write(b'</Cells>\n')
if elements is not None:
(types,connectivity,offsets) = elements
f.write(b'<Cells>\n')
_write_array(f,connectivity,b'Name="connectivity"')
_write_array(f,offsets,b'Name="offsets"')
_write_array(f,types,b'Name="types"')
f.write(b'</Cells>\n')
if data is not None :
f.write(b'<PointData>\n')
for name,v in data :
vc = np.ascontiguousarray(v)
appended = _write_array(f,appended,vc,8,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),vc.shape[-1]))
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</PointData>\n')
if cell_data is not None :
f.write(b'<CellData>\n')
for name,v in cell_data :
vc = np.ascontiguousarray(v)
appended = _write_array(f,appended,vc,8,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),vc.shape[-1]))
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</CellData>\n')
f.write(b'</Piece>\n')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
vc = np.ascontiguousarray(v)
appended = _write_array(f,appended,vc,4,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,vc.shape[0],vc.shape[1]))
_write_array(f,v,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,v.shape[0],v.shape[1]))
f.write(b'</FieldData>\n');
f.write(b'</UnstructuredGrid>\n')
f.write(b'<AppendedData encoding="base64">\n')
f.write(appended)
f.write(b'</AppendedData>\n')
f.write(b'</%s>\n' % vtkdescr)
f.write(b'</VTKFile>\n')
f.close()
def _get_array_info(a) :
ncomp = int(a["NumberOfComponents"]) if "NumberOfComponents" in a else None
ntuples = int(a["NumberOfTuples"]) if "NumberOfTuples" in a else None
return (a["Name"],a["type"],int(a["offset"]),(ntuples,ncomp))
def _get_binary_array(appended,pos,dtype,shape) :
data = appended[pos:pos+24]
_,n,_,s = struct.unpack("iiii",b64decode(data))
data = appended[pos+24:pos+24-(-s//3)*4]
data = b64decode(data)
data = zlib.decompress(data)
return np.frombuffer(data,dtype).reshape(shape)
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def read(fname) :
f = open(fname,"rb")
root = ET.parse(f).getroot()
grid = root.find("UnstructuredGrid")
root = ET.parse(open(fname,"rb")).getroot()
ftype = root.attrib["type"]
grid = root.find(ftype)
piece = grid.find("Piece")
npoints = int(piece.attrib["NumberOfPoints"])
ncells = int(piece.attrib["NumberOfCells"])
appended = root.find("AppendedData").text
shift = appended.find("_")+1
dim = 2
pointsoffset = int(piece.find("Points/DataArray").attrib["offset"])
x = _get_binary_array(appended,shift+pointsoffset,np.float64,(npoints,3))
for f in piece.find("Cells") :
if f.attrib["Name"] == "connectivity" :
connectivityoffset = int(f.attrib["offset"])
el = _get_binary_array(appended,shift+connectivityoffset,np.int32,(ncells,dim+1))
fdata = []
for f in grid.find("FieldData") :
name,ntype,offset,shape = _get_array_info(f.attrib)
v = _get_binary_array(appended,shift+offset,np.int32,shape)
fdata.append((name,v))
data = []
for f in piece.find("PointData") :
name,ntype,offset,shape = _get_array_info(f.attrib)
v = _get_binary_array(appended,shift+offset,np.float64,(npoints,shape[1]))
data.append((name,v))
return x,el,data,fdata
x = _read_array(piece.find("Points/DataArray"))
cells = piece.find("Cells")
fd = grid.find("FieldData")
cd = piece.find("CellData")
pd = piece.find("PointData")
c = dict((f.attrib["Name"],_read_array(f)) for f in cells) if cells else None
fdata = dict((f.attrib["Name"],_read_array(f)) for f in fd) if fd else None
cdata = dict((f.attrib["Name"],_read_array(f)) for f in cd) if cd else None
pdata = dict((f.attrib["Name"],_read_array(f)) for f in pd) if pd else None
return x,c,pdata,cdata,fdata
def write_multipart(basename,parts,i,t,reset_index=None):
filename = basename+"_%05d.vtm"%i
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="vtkMultiBlockDataSet" version="1.0" byte_order="LittleEndian">\n');
f.write(b'<vtkMultiBlockDataSet>\n');
for ip,n in enumerate(parts):
f.write(b'<DataSet index="%i" file="%s"/>\n'%(ip,n.encode()));
f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
......@@ -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),
......@@ -209,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.
......@@ -218,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)
......@@ -238,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"]
......@@ -256,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
......@@ -282,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
......@@ -292,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."""
......@@ -338,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
......@@ -355,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,22 +30,28 @@
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=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"""
......@@ -55,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)
......@@ -63,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.
......@@ -97,17 +103,21 @@ 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 radius(self):
"""Return the list of particle masses."""
return self._get_matrix("Particle", 2)[:, 0][:,None]
def position(self):
"""Return the list of particle centre position vectors."""
......@@ -117,24 +127,36 @@ class ParticleProblem :
"""Return the list of particle velocity vectors."""
return self._get_matrix("Velocity", self._dim)
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)
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")
......@@ -149,15 +171,81 @@ 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()
x = self.position()
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.radius()),("Velocity",v)]
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp")
xdisk = self._get_matrix("Disk",2*self._dim+1)[:,: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 :