Commit ebb4bb39 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

test pour jonjonjon

parent 176e57c6
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;
}
static int fluid_problem_import_vtk(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, "r");
if (!f){
printf("Cannot open file \"%s\".\n", filename);
return 1;
}
check_word_in_file(f,"<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n");
check_word_in_file(f,"<UnstructuredGrid>\n");
if (fscanf(f, "<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n", &(mesh->n_nodes), &(mesh->n_elements)) != 2){
printf("Cannot read physical names\n");
return 1;
}
check_word_in_file(f,"<Points>\n");
check_word_in_file(f,"<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < 3*mesh->n_nodes; ++i) {
if (fscanf(f,"%g ", &(mesh->x[i])) != 1){
printf("Cannot read nodes\n");
return 1;
}
}
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"</Points>\n");
check_word_in_file(f,"<Cells>\n");
check_word_in_file(f,"<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_elements*(DIMENSION+1); ++i) {
if (fscanf(f,"%i ", &(mesh->elements[i])) != 1){
printf("Cannot read elements.\n");
return 1;
}
}
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n");
char word[256];
for(int i = 0; i < mesh->n_elements;++i) {
if (fscanf(f,"%i ", &word) != 1){
printf("Cannot read element first node id\n");
return 1;
}
}
#if DIMENSION == 2
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"<DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_elements; ++i) {
if (fscanf(f,"%i ", word) != 1){
printf("Cannot read element type.\n");
return 1;
}
}
#else
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"<DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_elements; ++i) {
if (fscanf(f,"%i ", word) != 1){
printf("Cannot read element type.\n");
return 1;
}
}
#endif
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"</Cells>\n");
check_word_in_file(f,"<PointData Scalars=\"Porosity Pressure\" Vectors=\"Velocity\">\n");
check_word_in_file(f,"<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
if (fscanf(f,"%g ", &(problem->porosity[i])) != 1){
printf("Cannot read porosity\n");
return 1;
}
}
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
if (fscanf(f,"%g ", &(problem->solution[i*(DIMENSION+1)+DIMENSION])) != 1){
printf("Cannot read pressure\n");
return 1;
}
}
check_word_in_file(f,"\n</DataArray>\n");
check_word_in_file(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)
if (fscanf(f,"%g ", &(problem->solution[i*(DIMENSION+1)+j])) != 1){
printf("Cannot read velocity\n");
return 1;
}
}
fclose(f);
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, "%i %i %*i\n", &nd, &ns)!=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");
for (size_t i = 0; i < nd; ++i) {
Disk *d = vectorPush(&p->disks);
if (fscanf(f, "%.16g %.16g %.16g ", &(d->x[0]), &(d->x[1]), DIMENSION == 2 ? &scratch : &(d->x[2]))!=3){
printf("Cannot read disks\n");
return 1;
}
}
for (size_t i = 0; i < ns; ++i) {
Segment *s = vectorPush(&p->segments);
if (fscanf(f, "%.16g %.16g %.16g ", &(s->p[0][0]), &(s->p[0][1]), DIMENSION == 2 ? &scratch : &(s->p[0][2]))!=3){
printf("Cannot read segments\n");
return 1;
}
if (fscanf(f, "%.16g %.16g %.16g ", &(s->p[1][0]), &(s->p[1][1]), DIMENSION == 2 ? &scratch : &(s->p[1][2]))!=3){
printf("Cannot read segments\n");
return 1;
}
}
#if DIMENSION == 3
for (int i = 0; i < nt; ++i) {
Triangle *t = vectoPush(&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, "%.16g ", &(particles.r))!=1){
printf("Cannot read radii\n");
return 1;
}
}
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, "%.16g %.16g 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, "%.16g %.16g %.16g ", &(p->velocity[i * 3]), &(p->velocity[i * 3 + 1]), &(p->velocity[i * 3 + 2]))!=3){
printf("Cannot read particle velocities\n");
return 1;
}
#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, "%.16g %.16g 0 ", &(p->position[i * 2]), &(p->position[i * 2 + 1]))!=2){
printf("Cannot read particles positions\n");
return 1;
}
#else
vectorPushN(&p->position, DIMENSION);
if (fscanf(f, "%.16g %.16g %.16g ", &(p->position[i * 3]), &(p->position[i * 3 + 1]), &(p->position[i * 3 + 2]))!=3){
printf("Cannot read particles positions\n");
return 1;
}
#endif
}
fclose(f);
return 0;
}
static int _readContactVtkSingle(ParticleProblem *p, FILE *f, int ctype)
{
size_t n = 0;
size_t 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, "%.2g ", &(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;
}
}
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,"particle_particle");
_readContactVtkSingle(p,f,PARTICLE_DISK,"particle_disk");
_readContactVtkSingle(p,f,PARTICLE_SEGMENT,"particle_segment");
#if DIMENSION == 3
_readContactVtkSingle(p,f,PARTICLE_TRIANGLE,"particle_triangle");
#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;
}
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