Commit be161af7 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Merging master into friction validation

parents 3e4dbe64 d0a23bcf
Pipeline #4657 failed with stage
in 19 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)
# 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 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) :
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") :
"""Write VTK UnstructuredMesh (vtu) or PolyData (vtp) files.
The arrays are written in base64-encoded gzipped binary format.
Keyword arguments:
basename -- file path without the extension and the output number
e.g. "dir/file" becomes "dir/file_00012.vtp" and a collection vtk file "dir/file.pvd" will be created
i -- output number (will be appended to the filename and used as index in the collection)
t -- output time
elements -- list of elements. It can be None, otherwise it should be a tupple (types,connectivity,offsets)
where:
types -- is a vector of int with the type of each element
3:line,4:polyline,5:triangle,7:polygon,9:quad,10:tetra,12:hexahedron,...
see https://www.vtk.org/VTK/img/file-formats.pdf for a complete list
connectivity -- vector of int containing the indices of all elements nodes. Indices starts at 0.
offsets -- vector of int with the position of the *END* of the list of each element nodes inside the
connectivity vector
x -- two-dimensionnal array of float with the nodes coordinates, [[x0,y0,z0],[x1,y1,z1],...]].
it should have 3 columns, even for 2-dimensional meshes
node_data -- vector of pairs [(name0,data0),(name1,data1),...] with the name (as a string) and the fields attached
to the nodes. The field data should be two-dimensional arrays even if they have only one column.
cell_data -- like node_data, for the fields attached to the cells
field_data -- like node_data for the fields attached neither to cells or nodes.
reset_index -- boolean flag determining if the collection ".pvd" file index will be overwritten or appended.
if not specified, the .pvd is overwritten if i == 0.
vtktype -- "vtu" for UnstructuredMesh or "vtp" for PolyData. (The only difference seems to be that PolyData cannot
have 3-dimensional elements.)
"""
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')
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]))
f.write(b'</PointData>\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 :
_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 :
_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')
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]))
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')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
_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'</%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")
"""Read VTK UnstructuredMesh (vtu) or PolyData (vtp) files written by the write function
all arrays are returned as numpy array
Keyword arguments :
fname -- complete file path
Return (points,cells,node_data,cell_data,field_data) :
points -- nodes coordinates
cells -- dictionary with three entries : "connectivity", "offsets" and "types"
see the "write" function for a description of those arrays
node_data,cell_data,field_data -- dictionaries whose keys are the name of the node,cell and field data.
"""
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):
"""Write VTK MultiBlockDataSet (vtm) files.
Keyword arguments:
basename -- file path without the extension and the output number
e.g. "dir/file" becomes "dir/file_00012.vtm" and a collection vtk file "dir/file.pvd" will be created
i -- output number (will be appended to the filename and used as index in the collection)
t -- output time
parts -- list of data file paths
reset_index -- boolean flag determining if the collection ".pvd" file index will be overwritten or appended.
if not specified, the .pvd is overwritten if i == 0.
"""
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,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
......@@ -21,6 +21,7 @@
from pylmgc90 import chipy
import numpy as np
from . import VTK
import os
import math
import sys
......@@ -94,7 +95,7 @@ class ParticleProblem(object):
chipy.OpenPostproFiles()
chipy.OpenDisplayFiles(restart+1)
chipy.WriteDisplayFiles()
#chipy.WriteDisplayFiles()
chipy.ComputeMass()
......@@ -102,7 +103,6 @@ class ParticleProblem(object):
self._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._nbPoly = chipy.POLYG_GetNbPOLYG()
self._p2bMap = chipy.POLYG_GetPtrPOLYG2RBDY2()
self._position = np.zeros([self._nbDisk,3], 'd')
self._velocity = np.zeros([self._nbDisk,3], 'd')
......@@ -160,16 +160,41 @@ class ParticleProblem(object):
'gsit2' : 20
}
def setBoundaryVelocity(self, tag, v) :
def set_boundary_velocity(self, tag, v) :
self._tag2bnd[tag] = v
def write_vtk(self, filemame, t, i) :
chipy.WriteOutDof(1)
chipy.WriteOutVlocRloc(1)
def write_vtk(self, odir, iiter, t) :
v = self.velocity()
x = self.position()
if self._dim == 2 :
v = np.insert(v,2,0,axis=1)
x = np.insert(x,2,0,axis=1)
data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v)]
VTK.write(odir+"/particles",iiter,t,None,x,data,vtktype="vtp")
lmgcids = ("POLYG","PLANx","POLYR")
vtkids = (7,12,13)
def CHIPY(lid,fct) :
return getattr(chipy,lid+"_"+fct)()
def getx(l) :
x = CHIPY(l,"InitOutlines")
CHIPY(l,"UpdatePostdata")
if x.shape[1] == 2 : x = np.insert(x,2,0,axis=1)
return x.reshape((-1,3))
x = np.vstack(getx(l) for l in lmgcids)
connectivity = np.arange(0,x.shape[0],dtype=np.int32)