Skip to content
GitLab
Commit ce5c4863 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Revert "Merging master into friction validation" 

This reverts commit be161af7, reversing
changes made to 3e4dbe64.
parent 4fde75e0
Pipeline #4662 passed with stage
in 37 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ ce5c4863
......@@ -28,6 +28,7 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"
#define D DIMENSION
......@@ -468,8 +469,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.5;
double amin = 0.5;
double amax = 0;
double amin = 0;
double c = 0;
for (int i=0; i< N_N; ++i) {
double normup = 0;
......
paraview_utils/refresh_mb_output.py 0 → 100644
View file @ ce5c4863
# 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)
python/VTK.py
View file @ ce5c4863
# 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 _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))
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))
data = zlib.compress(a,2)
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()
# 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)
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)
def write(basename,i,t,elements,x,data,field_data) :
appended = b"_"
filename = "%s_%05d.vtu"%(basename,i)
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
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
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]
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
_write_array(f,x,b'NumberOfComponents="3"')
appended = _write_array(f,appended,x,8,b'NumberOfComponents="3"')
f.write(b'</Points>\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'<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')
f.write(b'</Piece>\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'<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')
f.write(b'</VTKFile>\n')
f.close()
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
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)
def read(fname) :
"""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)
f = open(fname,"rb")
root = ET.parse(f).getroot()
grid = root.find("UnstructuredGrid")
piece = grid.find("Piece")
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)
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
python/fluid.py
View file @ ce5c4863
......@@ -29,15 +29,15 @@
"""
from ctypes import *
import numpy as np
import numpy
import signal
import os
import sys
from . import VTK
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
lib2 = CDLL(os.path.join(dir_path,"libmbfluid2.so"))
lib3 = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
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 = 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])
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])
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 = np.require(a,"float","C")
tmp = numpy.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 = np.require(rho,"float","C").reshape([-1]).shape[0]
n_fluids = numpy.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,12 +192,9 @@ 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",v),
("velocity",self.solution()[:,:self._dim]),
("porosity",self.porosity()),
("old_porosity",self.old_porosity())
]
......@@ -208,10 +205,7 @@ 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))
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)
VTK.write(output_dir+"/fluid",it,t,self.elements(),self.coordinates(),data,field_data)
def import_vtk(self, filename) :
"""Read output file to reload computation.
......@@ -220,18 +214,17 @@ class FluidProblem :
filename -- name of the file to read
"""
def np2c(a,dtype) :
tmp = np.require(a,dtype,"C")
tmp = numpy.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])
x,el,data,fdata = VTK.read(filename)
self._lib.fluid_problem_private_set_elements(self._fp,
c_int(x.shape[0]),np2c(x,np.float64),
c_int(el.shape[0]),np2c(el,np.int32))
c_int(x.shape[0]),np2c(x,numpy.float64),
c_int(el.shape[0]),np2c(el,numpy.int32))
phys = {}
bnds = {}
for n,d in fdata.items() :
for n,d in fdata :
k,dim,pid,name = n.split(None,4)
if k == "Physical" :
phys[(dim,pid)] = (name,d)
......@@ -241,10 +234,11 @@ 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,np.int32),
c_int(bnd.shape[0]), np2c(bnd,np.int32))
c_int(nodes.shape[0]), np2c(nodes,numpy.int32),
c_int(bnd.shape[0]), np2c(bnd,numpy.int32))
sol = self.solution()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
data = dict(data)
sol[:,:self._dim] = data["velocity"]
sol[:,[self._dim]] = data["pressure"]
if self._n_fluids == 2 :
sol[:,[self._dim+1]] = data["concentration"]
......@@ -258,7 +252,7 @@ class FluidProblem :
Keyword arguments:
dt -- computation time step
"""
forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
forces = numpy.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
......@@ -284,7 +278,7 @@ class FluidProblem :
reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
def np2c(a) :
tmp = np.require(a,"float","C")
tmp = numpy.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
......@@ -294,7 +288,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 np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
return numpy.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
def solution(self) :
"""Give access to the fluid field value at the mesh nodes."""
......@@ -340,7 +334,7 @@ class FluidProblem :
ptr = f(self._fp)
size = fs(self._fp)
buf = (size*c_int).from_address(ptr)
return np.ctypeslib.as_array(buf)
return numpy.ctypeslib.as_array(buf)
def _n_physicals(self):
f = self._lib.fluid_problem_private_n_physical
......@@ -357,7 +351,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 = 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)
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)
return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd
python/lmgc90Interface.py
View file @ ce5c4863
......@@ -21,7 +21,6 @@
from pylmgc90 import chipy
import numpy as np
from . import VTK
import os
import math
import sys
......@@ -95,7 +94,7 @@ class ParticleProblem(object):
chipy.OpenPostproFiles()
chipy.OpenDisplayFiles(restart+1)
#chipy.WriteDisplayFiles()
chipy.WriteDisplayFiles()
chipy.ComputeMass()
......@@ -103,6 +102,7 @@ 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,41 +160,16 @@ class ParticleProblem(object):
'gsit2' : 20
}
def set_boundary_velocity(self, tag, v) :
def setBoundaryVelocity(self, tag, v) :
self._tag2bnd[tag] = v
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)
n = np.hstack(CHIPY(l,"GetNbPointOutlines")[1:]-CHIPY(l,"GetNbPointOutlines")[:-1] for l in lmgcids)
offsets = np.cumsum(n)