Commit 06d91fdd authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

split fluid_problem_io.c

parent cfe70a73
......@@ -3,7 +3,7 @@ ALL: libmbfluid3.so libscontact2.so libscontact3.so
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 src/mbxml.c src/yxml.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 src/fluid_problem_io.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
......
......@@ -6,7 +6,7 @@
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include <mbxml.h>
#include "mbxml.h"
#define DIMENSION 3
#define n_fields (DIMENSION+1)
......@@ -15,7 +15,7 @@
#define newton_atol 1e-6
#define newton_rtol 1e-5
#if DIMENSION == 2
#if DIMENSION==2
#define N_SF 3
#define N_N 3
......@@ -100,243 +100,6 @@ static double invDD(const double m[3][3], double inv[3][3]) {
#endif
static int check_word_in_file(FILE *f, const char *w) {
char word[256];
if(fscanf(f,"%255s", word) != 1){
printf("Invalid mesh file (\"%s\" expected).\n", w);
return -1;
}
if(strcmp(word, w) != 0){
printf("Invalid mesh file (\"%s\" expected).\n", w);
return -1;
}
return 0;
}
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.msh",output_dir, iter);
FILE *f = fopen(file_name, "w");
if (!f){
printf("Cannot open file \"%s\" for writing.\n", file_name);
return -1;
}
if(mesh_write_msh(problem->mesh, f))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
return -1;
int comp[]={0,1,2};
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, DIMENSION))
return -1;
fclose(f);
return 0;
}
int fluid_problem_import_vtk(FluidProblem *problem, 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);
Mesh *m = problem->mesh;
m->n_nodes = -1, m->n_elements = -1;
for(size_t i = 0; i < pieceelement.n_attr; ++i) {
if (!strcasecmp("NumberOfPoints",pieceelement.attr[i].name))
sscanf(pieceelement.attr[i].value,"%i",&m->n_nodes);
if (!strcasecmp("NumberOfCells",pieceelement.attr[i].name))
sscanf(pieceelement.attr[i].value,"%i",&m->n_elements);
}
if (m->n_nodes<0 || m->n_elements < 0){
printf("NumberOfPoints or NumberOfCells not specified in vtk piece");
}
MBVtkDataArray a;
MbXmlElement e;
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name,"Points")){
mb_read_vtk_data_array(r,&a);
if(m->x) free(m->x);
m->x = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
mb_vtk_data_array_destroy(&a);
}
else if (!strcasecmp(e.name,"Cells")) {
MBVtkDataArray connectivity;
int ncon = -1;
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Connectivity")) {
connectivity = a;
}
else if(!strcasecmp(a.name,"Offsets")){
int *offsets = mb_vtk_data_array_to_int(&a,m->n_elements,1);
ncon = offsets[m->n_elements-1];
free(offsets);
mb_vtk_data_array_destroy(&a);
}
else if(!strcasecmp(a.name,"Types")){
mb_vtk_data_array_destroy(&a);
}
else {
printf("unknown data array \"%s\" in vtk piece\n", a.name);
exit(1);
}
}
}
else if (!strcasecmp(e.name,"PointData")) {
if (problem->solution)
free(problem->solution);
problem->solution = malloc(sizeof(double)*(DIMENSION+1)*m->n_nodes);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Porosity")) {
if(problem->porosity)
free(problem->porosity);
problem->porosity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
}
else if(!strcasecmp(a.name,"Pressure")) {
double *p = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
memcpy(problem->solution+DIMENSION*m->n_nodes,p,sizeof(double)*m->n_nodes);
free(p);
}
else if(!strcasecmp(a.name,"Velocity")) {
double *v = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
memcpy(problem->solution,v,sizeof(double)*m->n_nodes*DIMENSION);
free(v);
}
mb_vtk_data_array_destroy(&a);
}
}
else {
printf("unknown type in piece : \"%s\"\n", e.name);
exit(1);
}
mb_xml_end_element(r,&e);
}
mb_xml_end_element(r,&pieceelement);
MbXmlElement fielddataelement;
mb_xml_read_element(r,&fielddataelement);
int n_phys=0;
int n_phys_alloc = 10;
MBVtkDataArray *phys = malloc(sizeof(MBVtkDataArray*)*10);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strncasecmp(a.name,"Physical ",9)) {
if (n_phys == n_phys_alloc) {
n_phys_alloc *= 2;
phys = realloc(phys,n_phys_alloc*sizeof(MBVtkDataArray*));
}
}
else {
printf("unknown data array \"%s\" in vtk piece\n", a.name);
exit(1);
}
}
for (int i = 0; i < m->n_phys; ++i){
free(m->phys_name[i]);
free(m->phys_nodes[i]);
}
m->phys_id=realloc(m->phys_id,n_phys*sizeof(int));
m->phys_dimension=realloc(m->phys_dimension,n_phys*sizeof(int));
m->phys_n_nodes=realloc(m->phys_n_nodes,n_phys*sizeof(int));
m->phys_name = realloc(m->phys_name,n_phys*sizeof(char*));
m->phys_nodes = realloc(m->phys_name,n_phys*sizeof(int*));
for(int i = 0; i < n_phys; ++i) {
int p;
sscanf(a.name+9,"%i %i%n",&m->phys_dimension[i],&m->phys_id[i],&p);
m->phys_nodes[i] = mb_vtk_data_array_to_int(&phys[i],-1,-1);
m->phys_n_nodes[i] = phys[i].number_of_tuples;
mb_vtk_data_array_destroy(&phys[i]);
}
free(phys);
mb_xml_end_element(r,&fielddataelement);
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
mb_xml_reader_destroy(r);
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*(DIMENSION+1); ++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)*(DIMENSION+1));
}
#if DIMENSION == 2
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
}
#else
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 ", 10); // 10 for 3d
}
#endif
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Cells>\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,"<FieldData>\n");
for(int i = 0; i < mesh->n_phys; ++i) {
fprintf(f,"<DataArray Name=\"Physical %i %i %s\" NumberOfTuples=\"%i\" NomberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n",mesh->phys_dimension[i],mesh->phys_id[i],mesh->phys_name[i],mesh->phys_n_nodes[i]);
for (int j = 0; j < mesh->phys_n_nodes[i]; ++j) {
fprintf(f, "%i ", mesh->phys_nodes[i][j]);
}
fprintf(f,"\n</DataArray>\n");
}
fprintf(f,"</FieldData>\n");
fprintf(f,"</UnstructuredGrid>\n");
fprintf(f,"</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])
{
......
#include "fluid_problem.h"
#include "mbxml.h"
#include <string.h>
#define n_fields (DIMENSION+1)
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.msh",output_dir, iter);
FILE *f = fopen(file_name, "w");
if (!f){
printf("Cannot open file \"%s\" for writing.\n", file_name);
return -1;
}
if(mesh_write_msh(problem->mesh, f))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
return -1;
int comp[]={0,1,2};
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, DIMENSION))
return -1;
fclose(f);
return 0;
}
int fluid_problem_import_vtk(FluidProblem *problem, 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);
Mesh *m = problem->mesh;
m->n_nodes = -1, m->n_elements = -1;
for(size_t i = 0; i < pieceelement.n_attr; ++i) {
if (!strcasecmp("NumberOfPoints",pieceelement.attr[i].name))
sscanf(pieceelement.attr[i].value,"%i",&m->n_nodes);
if (!strcasecmp("NumberOfCells",pieceelement.attr[i].name))
sscanf(pieceelement.attr[i].value,"%i",&m->n_elements);
}
if (m->n_nodes<0 || m->n_elements < 0){
printf("NumberOfPoints or NumberOfCells not specified in vtk piece");
}
MBVtkDataArray a;
MbXmlElement e;
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name,"Points")){
mb_read_vtk_data_array(r,&a);
if(m->x) free(m->x);
m->x = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
mb_vtk_data_array_destroy(&a);
}
else if (!strcasecmp(e.name,"Cells")) {
MBVtkDataArray connectivity;
int ncon = -1;
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Connectivity")) {
connectivity = a;
}
else if(!strcasecmp(a.name,"Offsets")){
int *offsets = mb_vtk_data_array_to_int(&a,m->n_elements,1);
ncon = offsets[m->n_elements-1];
free(offsets);
mb_vtk_data_array_destroy(&a);
}
else if(!strcasecmp(a.name,"Types")){
mb_vtk_data_array_destroy(&a);
}
else {
printf("unknown data array \"%s\" in vtk piece\n", a.name);
exit(1);
}
}
}
else if (!strcasecmp(e.name,"PointData")) {
if (problem->solution)
free(problem->solution);
problem->solution = malloc(sizeof(double)*(DIMENSION+1)*m->n_nodes);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Porosity")) {
if(problem->porosity)
free(problem->porosity);
problem->porosity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
}
else if(!strcasecmp(a.name,"Pressure")) {
double *p = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
memcpy(problem->solution+DIMENSION*m->n_nodes,p,sizeof(double)*m->n_nodes);
free(p);
}
else if(!strcasecmp(a.name,"Velocity")) {
double *v = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
memcpy(problem->solution,v,sizeof(double)*m->n_nodes*DIMENSION);
free(v);
}
mb_vtk_data_array_destroy(&a);
}
}
else {
printf("unknown type in piece : \"%s\"\n", e.name);
exit(1);
}
mb_xml_end_element(r,&e);
}
mb_xml_end_element(r,&pieceelement);
MbXmlElement fielddataelement;
mb_xml_read_element(r,&fielddataelement);
int n_phys=0;
int n_phys_alloc = 10;
MBVtkDataArray *phys = malloc(sizeof(MBVtkDataArray*)*10);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strncasecmp(a.name,"Physical ",9)) {
if (n_phys == n_phys_alloc) {
n_phys_alloc *= 2;
phys = realloc(phys,n_phys_alloc*sizeof(MBVtkDataArray*));
}
}
else {
printf("unknown data array \"%s\" in vtk piece\n", a.name);
exit(1);
}
}
for (int i = 0; i < m->n_phys; ++i){
free(m->phys_name[i]);
free(m->phys_nodes[i]);
}
m->phys_id=realloc(m->phys_id,n_phys*sizeof(int));
m->phys_dimension=realloc(m->phys_dimension,n_phys*sizeof(int));
m->phys_n_nodes=realloc(m->phys_n_nodes,n_phys*sizeof(int));
m->phys_name = realloc(m->phys_name,n_phys*sizeof(char*));
m->phys_nodes = realloc(m->phys_name,n_phys*sizeof(int*));
for(int i = 0; i < n_phys; ++i) {
int p;
sscanf(a.name+9,"%i %i%n",&m->phys_dimension[i],&m->phys_id[i],&p);
m->phys_nodes[i] = mb_vtk_data_array_to_int(&phys[i],-1,-1);
m->phys_n_nodes[i] = phys[i].number_of_tuples;
mb_vtk_data_array_destroy(&phys[i]);
}
free(phys);
mb_xml_end_element(r,&fielddataelement);
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
mb_xml_reader_destroy(r);
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*(DIMENSION+1); ++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)*(DIMENSION+1));
}
#if DIMENSION == 2
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
}
#else
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 ", 10); // 10 for 3d
}
#endif
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Cells>\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,"<FieldData>\n");
for(int i = 0; i < mesh->n_phys; ++i) {
fprintf(f,"<DataArray Name=\"Physical %i %i %s\" NumberOfTuples=\"%i\" NomberOfComponents=\"1\" type=\"Int32\" format=\"ascii\">\n",mesh->phys_dimension[i],mesh->phys_id[i],mesh->phys_name[i],mesh->phys_n_nodes[i]);
for (int j = 0; j < mesh->phys_n_nodes[i]; ++j) {
fprintf(f, "%i ", mesh->phys_nodes[i][j]);
}
fprintf(f,"\n</DataArray>\n");
}
fprintf(f,"</FieldData>\n");
fprintf(f,"</UnstructuredGrid>\n");
fprintf(f,"</VTKFile>\n");
fclose(f);
return 0;
}
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