Commit 4d2a7402 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

vtk output

parent 28e42b3e
......@@ -11,6 +11,7 @@ class FluidSolver :
def __init__(self, mesh, dt, rhop, rho, mu, g, imex=True) :
if isinstance(mesh, str):
mesh = Mesh(mesh)
self._mesh = mesh
self._jac = MeshJacobian(mesh)
self._pf = ParticleFluid(self._jac)
self._imex = imex
......@@ -84,16 +85,64 @@ class FluidSolver :
print("unknown export format : \"%s\"." % filetype)
def write_solution(self, odir, oiter, time, filetype="msh") :
dim = self._jac.dim
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(self._jac.dim, self._sol))
X = np.zeros(((self._jac.dim + 1) *3, self._jac.x.shape[0]))
for d in range (self._jac.dim) :
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(dim, self._sol))
X = np.zeros(((dim+1)*3, self._jac.x.shape[0]))
for d in range (dim) :
X[d::3, :] = self._dof.getField(d, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
elif filetype == "vtk" :
output = open("%s/f_%05i.vtu" %(odir, oiter), "w")
vertices = np.array(self._mesh.vertices)[:,:3]
elements = np.array([[v[3] for v in e.vertices] for e in self._dof.getElements(0)], np.int32) -1
v = np.zeros((vertices.shape[0], 3))
for d in range (dim) :
v[elements, d] = self._jac.x[:,:,d]
output.write("<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n")
output.write("<UnstructuredGrid>\n")
output.write("<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n" % (vertices.shape[0], elements.shape[0]))
output.write("<Points>\n")
output.write("<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n")
np.savetxt(output, v)
output.write("</DataArray>\n")
output.write("</Points>\n")
output.write("<Cells>\n")
output.write("<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n")
np.savetxt(output, elements, fmt="%i")
output.write("</DataArray>\n")
output.write("<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n")
np.savetxt(output, (np.array(range(elements.shape[0]), np.int32)+1)*(dim+1), fmt="%i", newline=" ")
output.write("</DataArray>\n")
output.write("<DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">\n")
t = np.ndarray((elements.shape[0]), np.int32)
t[:] = 5 if dim == 2 else 10
np.savetxt(output, t, fmt="%i", newline=" ")
output.write("</DataArray>\n")
output.write("</Cells>\n")
output.write("<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n")
output.write("<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
np.savetxt(output, self._pf.porosity)
output.write("</DataArray>\n")
output.write("<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
p = np.ndarray((vertices.shape[0]))
p[np.transpose(elements)] = self._dof.getField(dim, self._sol)
np.savetxt(output, p)
output.write("</DataArray>\n")
output.write("<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">")
v = np.zeros((vertices.shape[0], 3))
for d in range (self._jac.dim) :
v[np.transpose(elements), d] = self._dof.getField(d, self._sol)[:self._jac.dim+1,:]
np.savetxt(output, v)
output.write("</DataArray>\n")
output.write("</PointData>\n")
output.write("</Piece>\n")
output.write("</UnstructuredGrid>\n")
output.write("</VTKFile>\n")
else :
print("unknown export format : \"%s\"." % filetype)
print("unknown export format : \"%s\" (choose \"vtk\" or \"msh\")." % filetype)
def mesh_position(self) :
return self._jac.x
......
......@@ -473,7 +473,7 @@ void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, doub
}
}*/
if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN!!! %g %g\n", i, x[0], x[1]);
printf(" PARTICLE %lu OUTSIDE DOMAIN!!! %g %g\n", i, x[0], x[1]);
exit(1);
}
//}
......@@ -764,21 +764,97 @@ void particleProblemWrite(const ParticleProblem *p, const char *filename) {
}
fclose(output);
}
void particleProblemWriteVtu(ParticleProblem *p, const char *filename) {
FILE *f = fopen(filename, "w");
fprintf(f, "<VTKFile type=\"UnstructuredGrid\">\n<UnstructuredGrid>\n<Piece NumberOfPoints=\"%i\" NumberOfCells=\"0\">\n", vectorSize(p->particles));
fprintf(f, "<PointData Scalars=\"Radius\">\n<DataArray Name=\"Radius\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
fprintf(f, "<VTKFile type=\"PolyData\">\n");
fprintf(f, "<PolyData>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\">\n", vectorSize(p->particles));
fprintf(f, "<PointData Scalars=\"Radius\" Vectors=\"Velocity\">\n");
fprintf(f, "<DataArray Name=\"Radius\" NumberOfComponents = \"1\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->particles); ++i) {
fprintf(f, "%g ", p->particles[i].r);
fprintf(f, "%.16g ", p->particles[i].r);
}
fprintf(f, "</DataArray>\n</PointData>\n");
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->particles); ++i) {
fprintf(f, "%g %g 0. ", p->position[i * DIMENSION + 0], p->position[i * DIMENSION + 1]);
#if DIMENSION == 2
fprintf(f, "%.16g %.16g 0 ", p->velocity[i * 2], p->velocity[i * 2 + 1]);
#else
fprintf(f, "%.16g %.16g %.16g ", p->velocity[i * 3], p->velocity[i * 3 + 1], p->velocity[i * 3 + 2]);
#endif
}
fprintf(f, "</DataArray>\n</Points>\n");
fprintf(f, "<Cells>\n<DataArray type=\"Int32\" Name=\"connectivity\" />\n<DataArray type=\"Int32\" Name=\"offsets\" />\n<DataArray type=\"UInt8\" Name=\"types\" />\n</Cells>\n");
fprintf(f, "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n");
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</PointData>\n");
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->position); ++i) {
#if DIMENSION == 2
fprintf(f, "%.16g %.16g 0 ", p->position[i * 2], p->position[i * 2 + 1]);
#else
fprintf(f, "%.16g %.16g %.16g ", p->position[i * 3], p->position[i * 3 + 1], p->position[i * 3 + 2]);
#endif
}
fprintf(f, "\n</DataArray>\n</Points>\n");
//fprintf(f, "<Cells>\n<DataArray type=\"Int32\" Name=\"connectivity\" />\n<DataArray type=\"Int32\" Name=\"offsets\" />\n<DataArray type=\"UInt8\" Name=\"types\" />\n</Cells>\n");
//fprintf(f, "</Piece>\n");
fprintf(f, "</PolyData>\n");
fprintf(f, "</VTKFile>\n");
fclose(f);
}
void particleProblemWriteBoundaryVtu(ParticleProblem *p, const char *filename) {
FILE *f = fopen(filename, "w");
fprintf(f, "<VTKFile type=\"PolyData\">\n");
fprintf(f, "<PolyData>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfLines=\"%lu\">\n", vectorSize(p->segments)*2, vectorSize(p->segments));
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->segments); ++i) {
Segment *s = &p->segments[i];
fprintf(f, "%.16g %.16g %.16g ", s->p[0][0], s->p[0][1], DIMENSION == 2 ? 0. : s->p[0][2]);
fprintf(f, "%.16g %.16g %.16g ", s->p[1][0], s->p[1][1], DIMENSION == 2 ? 0. : s->p[1][2]);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "<Lines>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"connectivity\">\n");
for (int i = 0; i < vectorSize(p->segments); ++i) {
fprintf(f, "%i %i ", i*2, i*2+1);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"offsets\">\n");
for (int i = 0; i < vectorSize(p->segments); ++i) {
fprintf(f, "%i ", (i+1)*2);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Lines>\n");
fprintf(f, "</Piece>\n");
#if DIMENSION == 3
fprintf(f, "<Piece NumberOfPoints=\"%lu\" NumberOfPolys=\"%lu\">\n", vectorSize(p->triangles)*3, vectorSize(p->triangles));
fprintf(f, "<Points>\n<DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
for (int i = 0; i < vectorSize(p->triangles); ++i) {
Triangle *t = &p->triangles[i];
fprintf(f, "%.16g %.16g %.16g ", t->p[0][0], t->p[0][1], t->p[0][2]);
fprintf(f, "%.16g %.16g %.16g ", t->p[1][0], t->p[1][1], t->p[1][2]);
fprintf(f, "%.16g %.16g %.16g ", t->p[2][0], t->p[2][1], t->p[2][2]);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Points>\n");
fprintf(f, "<Polys>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"connectivity\">\n");
for (int i = 0; i < vectorSize(p->triangles)*3; ++i) {
fprintf(f, "%i ", i);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "<DataArray type=\"Int32\" Name=\"offsets\">\n");
for (int i = 0; i < vectorSize(p->triangles); ++i) {
fprintf(f, "%i ", (i+1)*3);
}
fprintf(f, "\n</DataArray>\n");
fprintf(f, "</Polys>\n");
fprintf(f, "</Piece>\n");
#endif
fprintf(f, "</PolyData>\n");
fprintf(f, "</VTKFile>\n");
fclose(f);
}
......
......@@ -28,6 +28,7 @@ double *particleProblemExternalForces(ParticleProblem *p);
unsigned long int particleProblemNParticle(ParticleProblem *p);
void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax);
void particleProblemWriteVtu(ParticleProblem *p, const char *filename);
void particleProblemWriteBoundaryVtu(ParticleProblem *p, const char *filename);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
int *particleProblemSegmentTag(ParticleProblem *p);
......
......@@ -143,6 +143,7 @@ struct ParticleProblem{};
PyObject *externalForces() {return coordVector2PyArray(particleProblemExternalForces($self->_this), particleProblemNParticle($self->_this), DIMENSION);}
void useJacobi(int jacobi, double relax = 0.5) {particleProblemUseJacobi($self->_this, jacobi, relax);}
void writeVtu(const char *filename) {particleProblemWriteVtu($self->_this, filename);}
void writeBoundaryVtu(const char *filename) {particleProblemWriteBoundaryVtu($self->_this, filename);}
};
PyObject *particleToMesh(PyObject *position, PyObject *elements);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment