Skip to content
GitLab
Commit 0e32d2ca authored by Matthieu Constant's avatar Matthieu Constant
Browse files

Merge branch 'master' into surfaceTension

parents f20d6e84 160a62c2
Expand all Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 0e32d2ca
......@@ -28,7 +28,6 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"
#define D DIMENSION
......@@ -328,7 +327,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
free(h0);
}
static void f_outflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid){
static void f_pressure(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
double p = s[P];
double u[D], du[D][D], dp[D];
for (int iD = 0; iD < D; ++iD) {
......@@ -350,6 +349,7 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
utau[i] += u[j]*tau[i][j];
}
}
f[P] = un;
if(problem->n_fluids == 1){
mu = problem->mu[0];
......@@ -360,14 +360,14 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a+problem->rho[1]*(1-a);
mu = problem->mu[0]*a+problem->mu[1]*(1-a);
f[A] = un*(un > 0 ? a : 1.);
f[A] = un*(un > 0 ? a : v[1]);
}
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
f[U+id] = v[0]*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
}
}
static void f_inflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid){
static void f_velocity(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
double p = s[P];
double u[D], du[D][D], dp[D];
for (int iD = 0; iD < D; ++iD) {
......@@ -379,11 +379,13 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
double divu = 0;
double tau[D][D];
double utau[D]={0.};
double rho, mu;
double rho, mu, rho_ext;
double unint = 0;
double un = 0;
for (int i = 0; i < D; ++i) {
un += u[i]*n[i];
un += c*v[i]*n[i];
unint += u[i]*n[i];
divu += du[i][i];
for (int j = 0; j < D; ++j) {
tau[i][j] = du[i][j]-u[i]*dc[j]/c;
......@@ -391,25 +393,27 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
}
}
f[P] = un;
if(problem->n_fluids == 1){
mu = problem->mu[0];
rho = problem->rho[0];
rho_ext = rho;
}
else{
double a = s[A];
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a+problem->rho[1]*(1-a);
rho_ext = problem->rho[0]*v[D]+problem->rho[1]*(1-v[D]);
mu = problem->mu[0]*a+problem->mu[1]*(1-a);
f[A] = un*(un > 0 ? a : 1.);
f[A] = un*(un > 0 ? a : v[D]);
//printf("%g\n",v[D]);
}
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id]+ (un> 0 ? rho*u[id]*un/c : 1);
f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*unint : rho_ext*un*un*n[id]);
}
}
static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc,const double dt,int eid)
static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc, const double dt, int eid, const double *v)
{
double p = s[P];
double u[D], du[D][D], dp[D];
......@@ -419,19 +423,30 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
}
double un = 0;
double a;
for (int i = 0; i < D; ++i) {
un += u[i]*n[i];
}
double rho;
if (problem->n_fluids == 1){
rho = problem->rho[0];
}
else{
a = s[A];
a = fmin(1.,fmax(0.,a));
rho = problem->rho[0]*a + problem->rho[1]*(1-a);
}
f[P] = 0;
for (int id = 0; id < D; ++id) {
f[U+id] = p*n[id];
}
if(problem->n_fluids == 2){
double a = s[A];
f[A] = 0;
}
for (int id = 0; id < D; ++id) {
f[U+id] = (v==NULL ? p : v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
}
}
static double pow2(double a) {return a*a;}
......@@ -453,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;
......@@ -568,12 +583,12 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
f1[P*D+i] = -u[i] + (stab_param==0 ? taup*R : stab_param*dp[i]);
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(extraa);
}
}
/*if(problem->n_fluids == 2){
if(problem->n_fluids == 2){
f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
}*/
}
}
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
......@@ -588,6 +603,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
double *f0 = malloc(sizeof(double)*n_fields);
double *g0 = malloc(sizeof(double)*n_fields);
double *h0 = malloc(sizeof(double)*n_fields*D);
int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
f_cb *f;
#if DIMENSION == 2
int elbnd[3][2] = {{0,1},{1,2},{2,0}};
......@@ -595,19 +611,35 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
for (int i =0; i < mesh->n_nodes; ++i)
bnd_node_local_id[i] = -1;
if(mesh->phys_dimension[iphys] != D-1) continue;
double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
int j = mesh->phys_nodes[iphys][i];
bnd_node_local_id[j] = i;
for (int k = 0; k < D; ++k)
x[i*D+k] = mesh->x[j*3+k];
}
f = NULL;
int n_value;
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary bnd = problem->weak_boundaries[ibnd];
if (strcmp(mesh->phys_name[iphys],bnd.tag) == 0){
f = bnd.cb;
WeakBoundary *bnd = problem->weak_boundaries + ibnd;
if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
f = bnd->cb;
n_value = bnd->n_value;
if (n_value > 0){
v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
bnd->apply(mesh->phys_n_nodes[iphys], x, v);
}
break;
}
}
if (f == NULL) {
f = f_wall_pressure;
//printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
//exit(1);
printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
exit(1);
}
for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
......@@ -675,9 +707,21 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
double c = 0;
double *vbnd = NULL;
if (n_value>0){
vbnd = malloc(sizeof(double)*n_value);
}
for (int i = 0; i < n_value; ++i){
vbnd[i] = 0;
}
double dc[D] = {0};
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
int node_local_id = bnd_node_local_id[el[i]];
if (node_local_id != -1)
for (int j = 0; j < n_value; j++){
vbnd[j] += v[node_local_id*n_value+j]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
s[j] += phi[i]*dof;
......@@ -689,7 +733,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
const double jw = LQW[iqp]*detbnd;
f(problem,n,f0,s,ds,c,dc,dt,iel);
f(problem,n,f0,s,ds,c,dc,dt,iel,vbnd);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -698,12 +742,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
f(problem,n,g0,s,ds,c,dc,dt,iel);
f(problem,n,g0,s,ds,c,dc,dt,iel,vbnd);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel);
f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel,vbnd);
ds[jfield*D+jd] = store;
}
int N2 = n_fields*N_SF;
......@@ -723,8 +767,11 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
}
}
free(vbnd);
}
}
free(x);
free(v);
}
free(s);
free(ds);
......@@ -1505,7 +1552,7 @@ int fluid_problem_n_nodes(const FluidProblem *p)
return p->mesh->n_nodes;
}
void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, int field, int row, StrongBoundaryCallback *apply)
void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, int field, int row, BoundaryCallback *apply)
{
problem->n_strong_boundaries++;
problem->strong_boundaries = realloc(problem->strong_boundaries,sizeof(StrongBoundary)*problem->n_strong_boundaries);
......@@ -1516,7 +1563,24 @@ void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, i
problem->strong_boundaries[i].row = row;
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype)
void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
printf("Weak boundary is already defined for tag \"%s\".\n", tag);
exit(1);
}
}
p->n_weak_boundaries++;
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
p->weak_boundaries[p->n_weak_boundaries-1].tag = strdup(tag);
p->weak_boundaries[p->n_weak_boundaries-1].n_value = (pressurecb ? 1 : 0);
p->weak_boundaries[p->n_weak_boundaries-1].cb = f_wall;
p->weak_boundaries[p->n_weak_boundaries-1].apply = pressurecb;
}
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
{
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
......@@ -1528,17 +1592,24 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const cha
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
p->weak_boundaries[p->n_weak_boundaries-1].tag = strdup(tag);
f_cb *f = NULL;
if (strcasecmp(bndtype,"wall") == 0)
f = f_wall_pressure;
else if (strcasecmp(bndtype,"inflow") == 0)
f = f_inflow;
else if (strcasecmp(bndtype,"outflow") == 0)
f = f_outflow;
if (strcasecmp(bndtype,"wall") == 0){
f = f_wall;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = 0;
}
else if (strcasecmp(bndtype,"velocity") == 0){
f = f_velocity;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = D+p->n_fluids-1;
}
else if (strcasecmp(bndtype,"pressure") == 0){
f = f_pressure;
p->weak_boundaries[p->n_weak_boundaries-1].n_value = p->n_fluids;
}
else {
printf("Unkown weak boundary type \"%s\".\n", bndtype);
exit(1);
}
p->weak_boundaries[p->n_weak_boundaries-1].cb = f;
p->weak_boundaries[p->n_weak_boundaries-1].apply = apply;
}
int fluid_problem_n_elements(const FluidProblem *p)
......
fluid/fluid_problem.h
View file @ 0e32d2ca
......@@ -34,19 +34,21 @@
#define N_N 4
#endif
typedef void StrongBoundaryCallback(int n, const double *x, double *bnd);
typedef void BoundaryCallback(int n, const double *x, double *bnd);
typedef struct {
char *tag;
int row;
int field;
StrongBoundaryCallback *apply;
BoundaryCallback *apply;
} StrongBoundary;
typedef struct FluidProblem FluidProblem;
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid);
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, const double *v);
typedef struct {
char *tag;
f_cb *cb;
int n_value;
BoundaryCallback *apply;
} WeakBoundary;
......@@ -98,12 +100,13 @@ int fluid_problem_n_fields(const FluidProblem *p);
double *fluid_problem_solution(const FluidProblem *p);
double *fluid_problem_coordinates(const FluidProblem *p);
int fluid_problem_n_nodes(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply);
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype);
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, StrongBoundaryCallback *apply);
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb);
int fluid_problem_private_n_physical(const FluidProblem *p);
void fluid_problem_private_physical(const FluidProblem *p, int i, char **phys_name, int *phys_dim, int *phys_id, int *phys_n_nodes, int **phys_nodes, int *phys_n_boundaries, int **phys_boundaries);
......
paraview_utils/refresh_mb_output.py deleted 100644 → 0
View file @ f20d6e84
# 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 @ 0e32d2ca
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") :
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)