Commit 2d72f7a7 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove scontactplot and export vtk

parent bb5b56ec
......@@ -45,6 +45,9 @@ class fluid_problem :
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
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 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))
......
......@@ -839,8 +839,13 @@ void particleProblemWrite(const ParticleProblem *p, const char *filename) {
fclose(output);
}
void particleProblemWriteVtk(ParticleProblem *p, const char *filename) {
int particleProblemWriteVtk(ParticleProblem *p, const char *filename) {
FILE *f = fopen(filename, "w");
if (!f){
printf("cannot open file %s\n", filename);
return -1;
}
fprintf(f, "<VTKFile type=\"PolyData\">\n");
fprintf(f, "<PolyData>\n");
fprintf(f, "<Piece NumberOfPoints=\"%lu\">\n", vectorSize(p->particles));
......
......@@ -27,7 +27,7 @@ double *particleProblemPosition(ParticleProblem *p);
double *particleProblemExternalForces(ParticleProblem *p);
unsigned long int particleProblemNParticle(ParticleProblem *p);
void particleProblemUseJacobi(ParticleProblem *p, int jacobi, double relax);
void particleProblemWriteVtk(ParticleProblem *p, const char *filename);
int particleProblemWriteVtk(ParticleProblem *p, const char *filename);
void particleProblemWriteBoundaryVtk(ParticleProblem *p, const char *filename);
double *particleProblemRadius(ParticleProblem *p);
int *particleProblemDiskTag(ParticleProblem *p);
......
......@@ -123,6 +123,115 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
return 0;
}
int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
Mesh *mesh = problem->mesh;
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.vtu",output_dir, iter);
FILE *f = fopen(file_name, "w");
if (!f){
printf("Cannot open file \"%s\" for writing.\n", file_name);
return -1;
}
fprintf(f,"<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n");
fprintf(f,"<UnstructuredGrid>\n");
fprintf(f,"<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n", mesh->n_nodes, mesh->n_elements);
fprintf(f,"<Points>\n");
fprintf(f,"<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < 3*mesh->n_nodes; ++i) {
fprintf(f,"%g ", mesh->x[i]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Points>\n");
fprintf(f,"<Cells>\n");
fprintf(f,"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_elements*3; ++i) {
fprintf(f,"%i ", mesh->elements[i]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n");
for(int i = 0; i < mesh->n_elements;++i) {
fprintf(f,"%i ", (i+1)*3);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_elements; ++i) {
fprintf(f,"%i ", 5); // 10 for 3d
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Cells>\n");
//fprintf(f,"<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n");
fprintf(f,"<PointData Scalars=\"Porosity Pressure\" Vectors=\"Velocity\">\n");
fprintf(f,"<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->porosity[i]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+DIMENSION]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int j = 0; j < DIMENSION; ++j)
fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+j]);
if(DIMENSION==2)
fprintf(f, "0 ");
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</PointData>\n");
fprintf(f,"</Piece>\n");
fprintf(f,"</UnstructuredGrid>\n");
fprintf(f,"</VTKFile>\n");
//t[:] = 5 if dim == 2 else 10
/*
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("<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)
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\">")
p = np.ndarray((vertices.shape[0]))
p[np.transpose(elements)] = self._jac.dofManager.getField(0, self._pf.porosity)
np.savetxt(output, p)
output.write("</DataArray>\n")
output.write("<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
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")
*/
fclose(f);
return 0;
}
#if DIMENSION == 2
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop, double gradp[2])
{
......
......@@ -36,6 +36,7 @@ typedef struct {
} FluidProblem;
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter);
int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir, double t, int iter);
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
......
cmake_minimum_required(VERSION 2.6)
project(scontactplot)
find_package(FLTK)
option(ENABLE_LIBPNG "required to save screenshot" ON)
add_definitions(-std=c++11)
if (FLTK_LIBRARIES)
add_executable(scontactplot fltk.cc ObjectCollection.cc simpleMesh.cc)
find_package(OpenGL REQUIRED)
include_directories(${OPENGL_INCLUDE_DIR})
include_directories(${FLTK_INCLUDE_DIR})
set(EXT_LIB ${OPENGL_LIBRARIES} ${FLTK_LIBRARIES})
if (ENABLE_LIBPNG)
find_package(PNG)
if (PNG_FOUND)
include_directories(${PNG_INCLUDE_DIR})
set(EXT_LIB ${EXT_LIB} ${PNG_LIBRARIES})
list(APPEND COMP_FLAGS "-DHAVE_LIBPNG")
endif(PNG_FOUND)
endif(ENABLE_LIBPNG)
target_link_libraries(scontactplot ${EXT_LIB})
install(TARGETS scontactplot DESTINATION bin)
if(COMP_FLAGS)
set_target_properties(scontactplot PROPERTIES COMPILE_FLAGS ${COMP_FLAGS})
endif(COMP_FLAGS)
endif (FLTK_LIBRARIES)
// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _GMSH_DEFINES_H_
#define _GMSH_DEFINES_H_
// IO file formats (numbers should not be changed)
#define FORMAT_MSH 1
#define FORMAT_UNV 2
#define FORMAT_GREF 3
#define FORMAT_XPM 4
#define FORMAT_PS 5
#define FORMAT_BMP 6
#define FORMAT_GIF 7
#define FORMAT_GEO 8
#define FORMAT_JPEG 9
#define FORMAT_AUTO 10
#define FORMAT_PPM 11
#define FORMAT_YUV 12
#define FORMAT_DMG 13
#define FORMAT_SMS 14
#define FORMAT_OPT 15
#define FORMAT_VTK 16
#define FORMAT_MPEG 17
#define FORMAT_TEX 18
#define FORMAT_VRML 19
#define FORMAT_EPS 20
#define FORMAT_PNG 22
#define FORMAT_PDF 24
#define FORMAT_POS 26
#define FORMAT_STL 27
#define FORMAT_P3D 28
#define FORMAT_SVG 29
#define FORMAT_MESH 30
#define FORMAT_BDF 31
#define FORMAT_CGNS 32
#define FORMAT_MED 33
#define FORMAT_DIFF 34
#define FORMAT_BREP 35
#define FORMAT_STEP 36
#define FORMAT_IGES 37
#define FORMAT_IR3 38
#define FORMAT_INP 39
#define FORMAT_PLY2 40
// Element types
#define TYPE_PNT 1
#define TYPE_LIN 2
#define TYPE_TRI 3
#define TYPE_QUA 4
#define TYPE_TET 5
#define TYPE_PYR 6
#define TYPE_PRI 7
#define TYPE_HEX 8
#define TYPE_POLYG 9
#define TYPE_POLYH 10
// Element types in .msh file format (numbers should not be changed)
#define MSH_LIN_2 1
#define MSH_TRI_3 2
#define MSH_QUA_4 3
#define MSH_TET_4 4
#define MSH_HEX_8 5
#define MSH_PRI_6 6
#define MSH_PYR_5 7
#define MSH_LIN_3 8
#define MSH_TRI_6 9
#define MSH_QUA_9 10
#define MSH_TET_10 11
#define MSH_HEX_27 12
#define MSH_PRI_18 13
#define MSH_PYR_14 14
#define MSH_PNT 15
#define MSH_QUA_8 16
#define MSH_HEX_20 17
#define MSH_PRI_15 18
#define MSH_PYR_13 19
#define MSH_TRI_9 20
#define MSH_TRI_10 21
#define MSH_TRI_12 22
#define MSH_TRI_15 23
#define MSH_TRI_15I 24
#define MSH_TRI_21 25
#define MSH_LIN_4 26
#define MSH_LIN_5 27
#define MSH_LIN_6 28
#define MSH_TET_20 29
#define MSH_TET_35 30
#define MSH_TET_56 31
#define MSH_TET_34 32
#define MSH_TET_52 33
#define MSH_POLYG_ 34
#define MSH_POLYH_ 35
#define MSH_QUA_16 36
#define MSH_QUA_25 37
#define MSH_QUA_36 38
#define MSH_QUA_12 39
#define MSH_QUA_16I 40
#define MSH_QUA_20 41
#define MSH_TRI_28 42
#define MSH_TRI_36 43
#define MSH_TRI_45 44
#define MSH_TRI_55 45
#define MSH_TRI_66 46
#define MSH_QUA_49 47
#define MSH_QUA_64 48
#define MSH_QUA_81 49
#define MSH_QUA_100 50
#define MSH_QUA_121 51
#define MSH_TRI_18 52
#define MSH_TRI_21I 53
#define MSH_TRI_24 54
#define MSH_TRI_27 55
#define MSH_TRI_30 56
#define MSH_QUA_24 57
#define MSH_QUA_28 58
#define MSH_QUA_32 59
#define MSH_QUA_36I 60
#define MSH_QUA_40 61
#define MSH_LIN_7 62
#define MSH_LIN_8 63
#define MSH_LIN_9 64
#define MSH_LIN_10 65
#define MSH_LIN_11 66
#define MSH_LIN_B 67
#define MSH_TRI_B 68
#define MSH_POLYG_B 69
#define MSH_LIN_C 70
#define MSH_TET_84 71
#define MSH_TET_120 72
#define MSH_TET_165 73
#define MSH_TET_220 74
#define MSH_TET_286 75
#define MSH_HEX_64 76
#define MSH_HEX_125 77
#define MSH_HEX_196 78
#define MSH_TET_74 79
#define MSH_TET_100 80
#define MSH_TET_130 81
#define MSH_TET_164 82
#define MSH_TET_202 83
#define MSH_NUM_TYPE 84
// Geometric entities
#define ENT_NONE 0
#define ENT_POINT (1<<0)
#define ENT_LINE (1<<1)
#define ENT_SURFACE (1<<2)
#define ENT_VOLUME (1<<3)
#define ENT_ALL (ENT_POINT | ENT_LINE | ENT_SURFACE | ENT_VOLUME)
// 2D meshing algorithms (numbers should not be changed)
#define ALGO_2D_MESHADAPT 1
#define ALGO_2D_MESHADAPT_OLD 4
#define ALGO_2D_DELAUNAY 5
#define ALGO_2D_FRONTAL 6
#define ALGO_2D_BAMG 7
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
#define ALGO_3D_FRONTAL 4
// Meshing methods
#define MESH_NONE 0
#define MESH_TRANSFINITE 1
#define MESH_UNSTRUCTURED 2
#endif
#include "GmshDefines.h"
#include "ObjectCollection.h"
#include <algorithm>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <cfloat>
#ifdef HAVE_LIBPNG
#include <png.h>
#endif
#ifdef __APPLE__
#include <OpenGL/gl.h>
#else
#include <GL/gl.h>
#endif
#include <stdarg.h>
#include <sys/stat.h>
static void rainbow(float f, float &r, float &g, float &b) {
f = std::max(f, (float)0.);
f = std::min(f, (float)1.);
float a=(1-f)/0.25;
int x =floor(a);
float y = a - x;
switch(x)
{
case 0: r=1;g=y;b=0;break;
case 1: r=1 -y ;g=1;b=0;break;
case 2: r=0;g=1;b=y;break;
case 3: r=0;g=1-y;b=1;break;
case 4: r=0;g=0;b=1;break;
}
}
class SolidSphere
{
protected:
std::vector<GLfloat> vertices;
std::vector<GLfloat> normals;
std::vector<GLushort> indices;
public:
SolidSphere(float radius, unsigned int rings, unsigned int sectors)
{
float const R = 1./(float)(rings-1);
float const S = 1./(float)(sectors-1);
int r, s;
vertices.resize(rings * sectors * 3);
normals.resize(rings * sectors * 3);
std::vector<GLfloat>::iterator v = vertices.begin();
std::vector<GLfloat>::iterator n = normals.begin();
for(r = 0; r < rings; r++) for(s = 0; s < sectors; s++) {
float const y = sin( -M_PI_2 + M_PI * r * R );
float const x = cos(2*M_PI * s * S) * sin( M_PI * r * R );
float const z = sin(2*M_PI * s * S) * sin( M_PI * r * R );
*v++ = x * radius;
*v++ = y * radius;
*v++ = z * radius;
*n++ = x;
*n++ = y;
*n++ = z;
}
indices.resize(rings * sectors * 4);
std::vector<GLushort>::iterator i = indices.begin();
for(r = 0; r < rings-1; r++) for(s = 0; s < sectors-1; s++) {
*i++ = r * sectors + s;
*i++ = r * sectors + (s+1);
*i++ = (r+1) * sectors + (s+1);
*i++ = (r+1) * sectors + s;
}
}
void draw()
{
glEnableClientState(GL_VERTEX_ARRAY);
glEnableClientState(GL_NORMAL_ARRAY);
glVertexPointer(3, GL_FLOAT, 0, &vertices[0]);
glNormalPointer(GL_FLOAT, 0, &normals[0]);
glDrawElements(GL_QUADS, indices.size(), GL_UNSIGNED_SHORT, &indices[0]);
glDisableClientState(GL_VERTEX_ARRAY);
glDisableClientState(GL_NORMAL_ARRAY);
}
};
void ObjectCollection::display()
{
glMatrixMode(GL_MODELVIEW);
glScalef(1./_scale, 1./_scale, 1./_scale);
glTranslatef(_shift[0], _shift[1], 0.);
if (dim == 3){
glRotatef(_rotation[1], 1, 0, 0);
glRotatef(_rotation[0], 0, 1, 0);
}
glEnable(GL_COLOR_MATERIAL);
if (dim == 3) {
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 };
GLfloat mat_shininess[] = { 20 };
glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
glShadeModel (GL_SMOOTH);
glEnable(GL_DEPTH_TEST);
GLfloat light_position[] = { -10.0, -10.0, 10.0, 0.0 };
glLightfv(GL_LIGHT0, GL_POSITION, light_position);
glClearColor(1., 1., 1., 1.);
} else {
glClearColor(1., 1., 1., 1.);
glEnable(GL_BLEND);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
}
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
glColor4f(1, 0.5, 0, 1);
glPushMatrix();
double S = *std::max_element(radius.begin(), radius.end())/2;
size_t zidx = 0;
if (dim == 3) {
int nTotPart = 5;
SolidSphere sphere(S, 15, 15);
for(size_t i = 0; i < radius.size()/nTotPart; ++i) {
if (coord[i * 3 + 1] < 0 && _clip) continue;
glPushMatrix();
glTranslatef(coord[i * 3 * nTotPart], coord[i * 3 * nTotPart + 1] , coord[i * 3 *nTotPart + 2]);
if (fixed[i]) glColor4f(1, 0.5, 0, 1);
else {
float aa = (coordorig[i *3 * nTotPart +2 ]+ 1.5)/(-0.35 +1.5);//(hypot(coordorig[i * 3 + 1], coordorig[i*3]));
float r, g, b;
rainbow(aa, r, g, b);
glColor4f(1, 0, 1, .5);
}
double f = radius[i * nTotPart]/S;
glScalef(f, f, f);
sphere.draw();
glPopMatrix();
}
}else {
int nTotPart = 1;
for(size_t i = 0; i < radius.size()/nTotPart; ++i) {
glEnableClientState(GL_VERTEX_ARRAY);
glVertexPointer(2, GL_FLOAT, 0, &circleGeom[0]);
glPushMatrix();
glTranslatef(coord[i * 2 * nTotPart], coord[i * 2 * nTotPart + 1], 0);
glScalef(radius[nTotPart*i], radius[nTotPart*i], 1.);
if(!fixed[i]) {
glColor4f(0.2, 0.2, 0.6, 0.3);
glDrawArrays(GL_TRIANGLE_FAN, 0, circleGeom.size() / 2);
}
glColor4f(0.4, 0.4, 0.8, 1);
glDrawArrays(GL_LINE_STRIP, 1, circleGeom.size() / 2 -1);
glDisableClientState(GL_VERTEX_ARRAY);
glPopMatrix();
}
}
glEnableClientState(GL_VERTEX_ARRAY);
glEnable(GL_BLEND);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glColor4f(1., 1., 1., .1);
glPolygonMode( GL_FRONT_AND_BACK, GL_FILL ); // GL_FILL
glVertexPointer(dim, GL_FLOAT, 0, &tricoord[0]);
glDrawArrays(GL_TRIANGLES, 0, tricoord.size() / dim);
glPolygonMode( GL_FRONT_AND_BACK, GL_FILL );
glDisable(GL_LIGHTING);
glDisable(GL_LIGHT0);
glColor4f(0., 0., 0., 1);
glVertexPointer(dim, GL_FLOAT, 0, &segcoord[0]);
glDrawArrays(GL_LINES, 0, segcoord.size() / dim);
glDisableClientState(GL_VERTEX_ARRAY);
auto &triangles = _mesh.elements(MSH_TRI_3);
auto &nodes = _mesh.nodes();
glColor4f(1.0, 0.0, 0.0, 1);
for (auto triangle : triangles) {
glBegin(GL_LINE_STRIP);
for (int i=0; i < 3; i++ ) {
int id = triangle.node(i);
glVertex3f(nodes[id*3],nodes[id*3+1],nodes[id*3+2]);
}
int id = triangle.node(0);
glVertex3f(nodes[id*3],nodes[id*3+1],nodes[id*3+2]);
glEnd();
}
/*
for (int i=0; i < _pressure.nData(); ++i) {
int idElement;
std::vector<double> data;
_pressure.data(i,idElement,data);
const simpleElement &triangle = *_mesh.element(idElement);
GLfloat r,g,b;
glBegin(GL_POLYGON);
for (size_t j=0; j < data.size(); ++j) {
double f = (data[j] - _pressure.min()) / (_pressure.max() - _pressure.min());
rainbow(f,r,g,b);
glColor3f(r,g,b);
int id = triangle.node(j);
glVertex3f(nodes[id*3]+0.03,nodes[id*3+1],nodes[id*3+2]);
}
glEnd();
}*/
glPopMatrix();
if (_ffmpeg) {
int xywh[4];
glGetIntegerv(GL_VIEWPORT, xywh);
int s = xywh[2] * xywh[3];
int *data = new int[s];
glReadPixels(xywh[0], xywh[1], xywh[2], xywh[3], GL_RGBA, GL_UNSIGNED_BYTE, data);
fwrite(data, sizeof(int) * s, 1, _ffmpeg);
delete []data;
}
}
static void readCoord(std::istream &stream, std::vector<float> &coordVector, int size)
{
for (int i = 0; i < size; ++i) {
double x;
stream >> x;
coordVector.push_back(x);
}
}
const std::string ObjectCollection::getTagName(size_t id) const {
if (id >= _tagname.size())
return "";
return _tagname[id];
}
int ObjectCollection::read(int step) {
if (step == -1)