Commit ab04e7a2 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

import vtk

parent ebb4bb39
ALL: libmbfluid3.so libscontact2.so libscontact3.so
CFLAGS=-Wno-unused-function -O3 -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt
CFLAGS=-Wno-unused-function -O3 -g -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt -Isrc
LDFLAGS=-shared -lm
FLUID_C=src/fluid_problem.c src/mesh.c src/mesh_find.c scontact/quadtree.c hxt/hxt_linear_system.c hxt/hxt_linear_system_lu.c hxt/hxt_message.c
FLUID_C=src/fluid_problem.c src/mesh.c src/mesh_find.c scontact/quadtree.c hxt/hxt_linear_system.c hxt/hxt_linear_system_lu.c hxt/hxt_message.c src/mbxml.c src/yxml.c
FLUID_H=src/tools.h src/fluid_problem.h hxt/hxt_linear_system.h src/mesh_find.h src/mesh.h scontact/quadtree.h scontact/vector.h
USE_PETSC=1
......@@ -13,8 +13,8 @@ ifeq ($(USE_PETSC),1)
LDFLAGS += -L${PETSC_DIR}/${PETSC_ARCH}/lib -Wl,-rpath=${PETSC_DIR}/${PETSC_ARCH}/lib -lpetsc
endif
SCONTACT_C=scontact/quadtree.c scontact/scontact.c
SCONTACT_H=scontact/quadtree.h scontact/scontact.h
SCONTACT_C=scontact/quadtree.c scontact/scontact.c src/mbxml.c src/yxml.c
SCONTACT_H=scontact/quadtree.h scontact/scontact.h src/mbxml.h src/yxml.h
libmbfluid2.so : ${FLUID_C} ${FLUID_H}
${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=2
......
......@@ -48,6 +48,9 @@ class fluid_problem :
def export_vtk(self, output_dir, t, it) :
lib.fluid_problem_export_vtk(self._fp, output_dir.encode(), c_double(t), c_int(it))
def import_vtk(self, filename) :
lib.fluid_problem_import_vtk(self._fp, filename.encode())
def compute_node_force(self, dt) :
forces = numpy.ndarray([self.n_particles,2],"d",order="C")
lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
......
......@@ -48,6 +48,9 @@ class fluid_problem :
def export_vtk(self, output_dir, t, it) :
lib.fluid_problem_export_vtk(self._fp, output_dir.encode(), c_double(t), c_int(it))
def import_vtk(self, filename) :
lib.fluid_problem_import_vtk(self._fp, filename.encode())
def compute_node_force(self, dt) :
def np2c(a) :
......
......@@ -935,43 +935,6 @@ void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], d
}
}
void particleProblemWrite(const ParticleProblem *p, const char *filename) {
FILE *output = fopen(filename, "w");
fprintf(output, "DIMENSION %i\n", DIMENSION);
fprintf(output, "TAGS %lu", vectorSize(p->tagname));
for (size_t i = 0; i < vectorSize(p->tagname); ++i) {
fprintf(output, " %s", p->tagname[i]);
}
fprintf(output, "\n");
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
fprintf(output, "T");
coordWrite(output, p->triangles[i].p[0]);
coordWrite(output, p->triangles[i].p[1]);
coordWrite(output, p->triangles[i].p[2]);
fprintf(output, " %i\n", p->triangleTag[i]);
}
#endif
for (size_t i = 0; i < vectorSize(p->segments); ++i){
fprintf(output, "S");
coordWrite(output, p->segments[i].p[0]);
coordWrite(output, p->segments[i].p[1]);
fprintf(output, " %i\n", p->segmentTag[i]);
}
for (size_t i = 0; i < vectorSize(p->disks); ++i) {
fprintf(output, "D");
coordWrite(output, p->disks[i].x);
fprintf(output, " %.16g %i\n", p->disks[i].r, p->diskTag[i]);
}
for(size_t i = 0; i< vectorSize(p->particles); ++i) {
fprintf(output, "P");
coordWrite(output, &p->position[i * DIMENSION]);
coordWrite(output, &p->velocity[i * DIMENSION]);
fprintf(output, " %.16g %.16g\n", p->particles[i].r, p->particles[i].m);
}
fclose(output);
}
static void _writeContactVtkSingle(ParticleProblem *p, FILE *f, int ctype, const char *cname)
{
size_t n = 0;
......@@ -1134,7 +1097,7 @@ static int particleProblemWriteBndVtu(ParticleProblem *p, const char *dirname, i
}
#if DIMENSION == 3
for (size_t i = 0; i < nt; ++i) {
fprintf(f, "%i %i %i", (int)(nd+ns*2+i*3+0), (int)(nd+ns*2+i*3+1), (int)(nd+ns*2+i*3+2));
fprintf(f, "%i %i %i ", (int)(nd+ns*2+i*3+0), (int)(nd+ns*2+i*3+1), (int)(nd+ns*2+i*3+2));
}
#endif
fprintf(f, "\n</DataArray>\n");
......@@ -1197,299 +1160,172 @@ int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id) {
return 0;
}
static int particleProblemReadBndVtu(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/boundaries_%05d.vtu",dirname, id);
size_t nd;
size_t ns;
size_t nt = 0;
size_t scratch;
FILE *f = fopen(filename, "r");
check_word_in_file(f, "<VTKFile type=\"UnstructuredGrid\">\n");
check_word_in_file(f, "<UnstructuredGrid>\n");
check_word_in_file(f, "<FieldData>\n");
check_word_in_file(f, "<DataArray Name=\"sizes\" NumberOfTuples=\"3\" NumberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n");
if (fscanf(f, "%lu %lu %lu\n", &nd, &ns, &scratch)!=3){
printf("Cannot read number of disks and segments\n");
return 1;
}
check_word_in_file(f, "</DataArray>\n");
check_word_in_file(f, "</FieldData>\n");
if (fscanf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"%lu\">\n", &scratch,&scratch)!=2){
printf("Cannot read number of points/cells\n");
return 1;
}
nt = scratch-nd-ns;
check_word_in_file(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
double scratchd;
for (size_t i = 0; i < nd; ++i) {
Disk *d = vectorPush(&p->disks);
#if DIMENSION == 2
if (fscanf(f, "%lf %lf %lf ", &d->x[0], &d->x[1], &scratchd)!=3){
printf("Cannot read disks\n");
return 1;
}
#else
if (fscanf(f, "%lf %lf %lf ", &d->x[0], &d->x[1], &d->x[2])!=3){
printf("Cannot read disks\n");
return 1;
}
#endif
}
for (size_t i = 0; i < ns; ++i) {
Segment *s = vectorPush(&p->segments);
#if DIMENSION == 2
if (fscanf(f, "%lf %lf %lf ", &s->p[0][0], &s->p[0][1], &scratchd)!=3){
printf("Cannot read segments\n");
return 1;
}
#else
if (fscanf(f, "%lf %lf %lf ", &s->p[0][0], &s->p[0][1], &s->p[0][2])!=3){
printf("Cannot read segments\n");
return 1;
}
#endif
if (fscanf(f, "%lf %lf %lf ", &s->p[1][0], &s->p[1][1], DIMENSION == 2 ? &scratchd : &s->p[1][2])!=3){
printf("Cannot read segments\n");
return 1;
}
}
#if DIMENSION == 3
for (int i = 0; i < nt; ++i) {
Triangle *t = vectorPush(&p->triangles);
if (fscanf(f, "%.16g %.16g %.16g ", &(t->p[0][0]), &(t->p[0][1]), &(t->p[0][2]))!=3||
fscanf(f, "%.16g %.16g %.16g ", &(t->p[1][0]), &(t->p[1][1]), &(t->p[1][2]))!=3||
fscanf(f, "%.16g %.16g %.16g ", &(t->p[2][0]), &(t->p[2][1]), &(t->p[2][2]))!=3){
printf("Cannot read triangles\n");
return 1;
}
}
#endif
fclose(f);
free(filename);
return 0;
}
static int particleProblemReadParticleVtp(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/particles_%05d.vtp",dirname, id);
FILE *f = fopen(filename, "r");
if (!f){
printf("cannot open file %s\n", filename);
return -1;
}
size_t n_pp_contacts = 0;
size_t n_p;
check_word_in_file(f, "<VTKFile type=\"PolyData\">\n");
check_word_in_file(f, "<PolyData>\n");
if (fscanf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfCells=\"0\">\n", &n_p)!=1){
printf("Cannot read number of particles\n");
return 1;
}
check_word_in_file(f, "<PointData Scalars=\"Radius\" Vectors=\"Velocity\">\n");
check_word_in_file(f, "<DataArray Name=\"Radius\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < n_p; ++i) {
Particle *particle = vectorPush(&p->particles);
if (fscanf(f, "%lf ", &particle->r)!=1){
printf("Cannot read radii\n");
return 1;
}
}
check_word_in_file(f, "\n</DataArray>\n");
check_word_in_file(f, "<DataArray Name=\"Mass\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < n_p; ++i) {
if (fscanf(f, "%lf ", &p->particles[i].m)!=1){
printf("Cannot read radii\n");
return 1;
#include "mbxml.h"
int particleProblemImportVtk(ParticleProblem *p, const char *filename)
{
MbXmlReader *r = mb_xml_reader_create(filename);
MbXmlElement fileelement,gridelement,pieceelement;
mb_xml_read_element(r,&fileelement);
mb_xml_read_element(r,&gridelement);
mb_xml_read_element(r,&pieceelement);
MbXmlElement e;
MBVtkDataArray a;
size_t n_particles;
for (size_t i = 0; i<pieceelement.n_attr; ++i) {
if(!strcasecmp(pieceelement.attr[i].name,"NumberOfPoints")){
sscanf(pieceelement.attr[i].value,"%lu", &n_particles);
vectorClear(p->position);
vectorClear(p->particles);
vectorClear(p->externalForces);
vectorClear(p->velocity);
vectorPushN(&p->position,n_particles*DIMENSION);
vectorPushN(&p->particles,n_particles);
vectorPushN(&p->externalForces,n_particles*DIMENSION);
vectorPushN(&p->velocity, n_particles*DIMENSION);
}
}
check_word_in_file(f, "\n</DataArray>\n");
check_word_in_file(f, "<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < n_p; ++i) {
#if DIMENSION == 2
vectorPushN(&p->velocity, DIMENSION);
if (fscanf(f, "%lf %lf 0 ", &p->velocity[i * 2], &p->velocity[i * 2 + 1])!=2){
printf("Cannot read particles velocities\n");
return 1;
}
#else
vectorPushN(&p->velocity, DIMENSION);
if (fscanf(f, "%lf %lf %lf ", &(p->velocity[i * 3]), &(p->velocity[i * 3 + 1]), &(p->velocity[i * 3 + 2]))!=3){
printf("Cannot read particle velocities\n");
return 1;
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name,"Points")){
mb_read_vtk_data_array(r,&a);
size_t n_particles = a.n/3;
for (size_t i = 0; i < n_particles; ++i)
for (int j = 0; j < DIMENSION; ++j)
p->position[DIMENSION*i+j] = ((double*)a.data)[i*3+j];
free(a.name);
free(a.data);
}
#endif
}
check_word_in_file(f, "\n</DataArray>\n");
check_word_in_file(f, "</PointData>\n");
check_word_in_file(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (size_t i = 0; i < n_p; ++i) {
#if DIMENSION == 2
vectorPushN(&p->position, DIMENSION);
if (fscanf(f, "%lf %lf 0 ", &(p->position[i * 2]), &(p->position[i * 2 + 1]))!=2){
printf("Cannot read particles positions\n");
return 1;
else if (!strcasecmp(e.name,"PointData")) {
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Velocity")) {
for (size_t i = 0; i < n_particles; ++i)
for (int j = 0; j < DIMENSION; ++j)
p->velocity[DIMENSION*i+j] = ((double*)a.data)[i*3+j];
}
else if(!strcasecmp(a.name,"Mass")) {
for (size_t i = 0; i < n_particles; ++i)
p->particles[i].m = ((double*)a.data)[i];
}
else if(!strcasecmp(a.name,"Radius")) {
for (size_t i = 0; i < n_particles; ++i){
p->particles[i].r = ((double*)a.data)[i];
}
}
free(a.name);
free(a.data);
}
}
#else
vectorPushN(&p->position, DIMENSION);
if (fscanf(f, "%lf %lf %lf ", &(p->position[i * 3]), &(p->position[i * 3 + 1]), &(p->position[i * 3 + 2]))!=3){
printf("Cannot read particles positions\n");
return 1;
else {
printf("unknown type in piece : \"%s\"\n", e.name);
exit(1);
}
#endif
mb_xml_end_element(r,&e);
}
fclose(f);
mb_xml_end_element(r,&pieceelement);
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
mb_xml_reader_destroy(r);
return 0;
}
static int _readContactVtkSingle(ParticleProblem *p, FILE *f, int ctype)
int particleProblemImportBoundaryVtk(ParticleProblem *p, const char *filename)
{
size_t n = 0;
char scratch;
if (fscanf(f, "<DataArray Name=\"%s\" NumberOfTuples=\"%zu\" NumberOfComponents=\"1\" type=\"Float32\" format=\"ascii\">\n",&scratch, &n)!=2){
printf("Cannot read number of contacts for impulsion\n");
return 1;
}
for (size_t i = 0; i < n; ++i) {
Contact *c = vectorPush(&p->contacts);
c[i].type = ctype;
if (fscanf(f, "%lf ", &(c[i].dv))!=1){
printf("Cannot read impulsion\n");
return 1;
}
}
check_word_in_file(f, "\n</DataArray>\n");
if (fscanf(f, "<DataArray Name=\"%s_idx\" NumberOfTuples=\"%zu\" NumberOfComponents=\"2\" type=\"UInt64\" format=\"ascii\">\n",&scratch,&n)!=2){
printf("Cannot read number of contacts for objects\n");
return 1;
}
for (size_t i = 0; i < n; ++i) {
if (fscanf(f, "%zu %zu ", &(p->contacts[i].o0), &(p->contacts[i].o1))!=2){
printf("Cannot read objects in contacts\n");
return 1;
MbXmlReader *r = mb_xml_reader_create(filename);
MbXmlElement fileelement,gridelement;
mb_xml_read_element(r,&fileelement);
mb_xml_read_element(r,&gridelement);
MbXmlElement e;
MBVtkDataArray a;
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name, "FieldData")) {
while(mb_read_vtk_data_array(r,&a)) {
free(a.name);
free(a.data);
}
}
}
check_word_in_file(f, "\n</DataArray>\n");
return 0;
}
static int particleProblemReadContactVtp(ParticleProblem *p, const char *dirname, int id) {
char *filename = malloc((strlen(dirname)+50)*sizeof(char));
sprintf(filename,"%s/contacts_%05d.vtp",dirname, id);
FILE *f = fopen(filename, "r");
if (!f){
printf("cannot open file %s\n", filename);
return -1;
}
check_word_in_file(f, "<VTKFile type=\"PolyData\">\n");
check_word_in_file(f, "<PolyData>\n");
check_word_in_file(f, "<FieldData>\n");
_readContactVtkSingle(p,f,PARTICLE_PARTICLE);
_readContactVtkSingle(p,f,PARTICLE_DISK);
_readContactVtkSingle(p,f,PARTICLE_SEGMENT);
else if(!strcasecmp(e.name,"Piece")) {
MbXmlElement f;
for (size_t i = 0; i<e.n_attr; ++i) {
MBVtkDataArray points={0}, connectivity={0}, offsets={0}, types={0};
while(mb_xml_read_element(r,&f)!=1) {
if (!strcasecmp(f.name,"Points"))
mb_read_vtk_data_array(r,&points);
else if (!strcasecmp(f.name,"Cells")) {
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"connectivity"))
connectivity = a;
else if (!strcasecmp(a.name,"types"))
types = a;
else if (!strcasecmp(a.name,"offsets"))
offsets = a;
else {
printf("unexpected data array in \"%s\" : \"%s\"\n",filename, a.name);
if(a.name) free(a.name);
if(a.data) free(a.data);
}
}
vectorClear(p->disks);
vectorClear(p->segments);
#if DIMENSION == 3
_readContactVtkSingle(p,f,PARTICLE_TRIANGLE);
vectorClear(p->triangles);
#endif
fclose(f);
free(filename);
return 0;
}
int particleProblemReadVtk(ParticleProblem *p, const char *dirname, int id) {
if (particleProblemReadParticleVtp(p,dirname,id) != 0)
return -1;
if (particleProblemReadBndVtu(p,dirname,id) != 0)
return -1;
if (particleProblemReadContactVtp(p,dirname,id) != 0)
return -1;
return 0;
}
void particleProblemLoad(ParticleProblem *p, const char *filename)
{
FILE *input = fopen(filename, "r");
char type[128];
int tag;
char stag[256];
while (fscanf(input, "%127s", type) == 1) {
if (!strcmp(type, "DIMENSION")) {
int dim;
fscanf(input, "%d", &dim);
if (dim != DIMENSION){
printf("model in file \"%s\" is of dimension %i\n", filename, dim);
exit(1);
}
}
else if (!strcmp(type, "TAGS")) {
for (size_t i = 0; i < vectorSize(p->tagname); ++i)
free(p->tagname[i]);
vectorFree(p->tagname);
p->tagname = NULL;
int n;
fscanf(input, "%d", &n);
for (int i = 0; i < n; ++i) {
fscanf(input, "%255s", stag);
if (strlen(stag) == 255) {
printf("tag name cannot exceed 255 characters\n");
exit(1);
double *pts = points.data;
int *off = ((int*)offsets.data);
for(size_t i = 0; i < types.n; ++i) {
int t = ((int*)types.data)[i];
int *nodes = ((int*)connectivity.data)+(i==0 ? 0 : off[i-1]);
if (t == 1){
particleProblemAddBoundaryDisk(p,pts+nodes[0]*3,0.,"");
}
else if (t == 3)
particleProblemAddBoundarySegment(p, pts+nodes[0]*3,pts+nodes[1]*3,"");
#if DIMENSION == 3
else if (t == 5)
particleProblemAddBoundaryTriangle(p, pts+nodes[0]*3,pts+nodes[1]*3,pts+nodes[2]*3,"");
#endif
else
printf("Unknown boundary object : %i\n", t);
}
free(points.name);
free(points.data);
free(connectivity.name);
free(connectivity.data);
free(offsets.name);
free(offsets.data);
free(types.name);
free(types.data);
}
else {
printf("unexpected element in \"%s\" : \"%s\"\n",filename, f.name);
exit(1);
}
mb_xml_end_element(r,&f);
}
particleProblemGetTagId(p, stag);
}
}
else if (!strcmp(type, "P")) {
double point[DIMENSION], v[DIMENSION];
double r, m;
coordRead(input, point);
coordRead(input, v);
fscanf(input, "%le %le", &r, &m);
particleProblemAddParticle(p, point, r, m);
for (int i = 0; i < DIMENSION; ++i)
p->velocity[(vectorSize(p->particles) - 1) * DIMENSION + i] = v[i];
}
else if (!strcmp(type, "D")) {
double point[DIMENSION];
double r;
coordRead(input, point);
fscanf(input, "%le %d", &r, &tag);
particleProblemAddBoundaryDiskTagId(p, point, r, tag);
}
else if (!strcmp(type, "S")) {
double x0[DIMENSION], x1[DIMENSION];
coordRead(input, x0);
coordRead(input, x1);
fscanf(input, "%d", &tag);
particleProblemAddBoundarySegmentTagId(p, x0, x1, tag);
}
#if DIMENSION == 3
else if (!strcmp(type, "T")) {
double x0[DIMENSION], x1[DIMENSION], x2[DIMENSION];
coordRead(input, x0);
coordRead(input, x1);
coordRead(input, x2);
fscanf(input, "%d", &tag);
particleProblemAddBoundaryTriangleTagId(p, x0, x1, x2, tag);
}
#endif
else {
printf("unkown type \"%s\" in file \"%s\".\n", type, filename);
exit(0);
printf("unexpected element in \"%s\" : \"%s\"\n",filename, e.name);
exit(1);
}
mb_xml_end_element(r,&e);
}
fclose(input);
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
mb_xml_reader_destroy(r);
return 0;
}
int particleProblemReadVtk(ParticleProblem *p, const char *dirname, int iter)
{
char *f0 = malloc(sizeof(char)*(strlen(dirname)+30));
char *f1 = malloc(sizeof(char)*(strlen(dirname)+30));
sprintf(f0, "%s/particles_%05d.vtp",dirname,iter);
sprintf(f1, "%s/boundaries_%05d.vtu",dirname,iter);
particleProblemImportVtk(p,f0);
particleProblemImportBoundaryVtk(p,f1);
free(f0);
free(f1);
return 0;
}
int *particleProblemDiskTag(ParticleProblem *p) {return p->diskTag;}
int *particleProblemSegmentTag(ParticleProblem *p) {return p->segmentTag;};
......
......@@ -28,6 +28,7 @@ double *particleProblemExternalForces(ParticleProblem *p);
unsigned long int particleProblemNParticle(ParticleProblem *p);
void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax);
int particleProblemWriteVtk(ParticleProblem *p, const char *dirname, int id);
int particleProblemReadVtk(ParticleProblem *p, const char *dirname, int iter);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
int *particleProblemSegmentTag(ParticleProblem *p);
......
......@@ -33,13 +33,11 @@ class ParticleProblem :
buf = (size*c_int).from_address(ptr)
return numpy.ctypeslib.as_array(buf)
def __init__(self, inputfile=None) :
def __init__(self) :
lib.particleProblemNew.restype = c_void_p
self._p = c_void_p(lib.particleProblemNew())
if self._p == None :
raise NameError("cannot create particle problem")
if inputfile :
lib.particleProblemLoad(self._p, inputfile.encode())
def volume(self):
return numpy.pi * self._get_matrix("Particle", 2)[:, [0]]**2
......@@ -76,10 +74,6 @@ class ParticleProblem :
self._get_matrix("ExternalForces",2)[:] = forces
lib.particleProblemIterate(self._p, c_double(numpy.min(self.r())), c_double(dt), c_double(1e-8), c_int(-1))