Commit 0abeac5a authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove c msh reader

parent 96c21b9c
Pipeline #8710 failed with stages
in 1 minute and 1 second
......@@ -21,8 +21,6 @@
set(SRC
file_reader.c
gmsh_mesh.c
fluid_problem.c
fluid_problem_concentration.c
mesh.c
......
#include "file_reader.h"
#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <errno.h>
struct FileReaderStruct {
FILE *f;
char *file_name;
char line[4096];
char *pos;
char *prev;
long int iline;
};
void file_reader_delete(FileReader *r) {
fclose(r->f);
free(r->file_name);
free(r);
}
void file_reader_error(FileReader *r, const char *msg, ...) {
printf("Error in file '%s' : ", r->file_name);
va_list args;
va_start(args, msg);
vprintf(msg,args);
printf("\n");
file_reader_delete(r);
exit(1);
}
void file_reader_error_at(FileReader *r, const char *msg,...) {
printf("Error in file '%s' line %lu column %lu : ", r->file_name, r->iline+1, r->prev-r->line+1);
va_list args;
va_start(args, msg);
vprintf(msg,args);
printf("\n");
printf("%s", r->line);
while(isblank(*r->prev))r->prev++;
for (int i = 0; i < (r->prev-r->line);++i) printf(" ");
printf("^\n");
file_reader_delete(r);
exit(1);
}
int file_reader_next_line_or_eof(FileReader *r) {
r->prev = r->pos;
if (r->pos){
while(isblank(*r->pos))r->pos++;
if (*r->pos != 0 && *r->pos != '\n')
file_reader_error_at(r,"end of line expected");
}
r->iline += 1;
char *l = fgets(r->line, 4096, r->f);
if (l == NULL) {
return 0;
}
r->pos = r->line;
return 1;
}
void file_reader_next_line_skip(FileReader *r) {
r->prev = r->pos;
r->iline += 1;
char *l = fgets(r->line, 4096, r->f);
if (l == NULL) {
file_reader_error(r, "Unexpected end of file.");
}
r->pos = r->line;
}
void file_reader_next_line(FileReader *r) {
r->prev = r->pos;
if (r->pos){
while(isblank(*r->pos))r->pos++;
if (*r->pos != 0 && *r->pos != '\n')
file_reader_error_at(r,"end of line expected");
}
r->iline += 1;
char *l = fgets(r->line, 4096, r->f);
if (l == NULL) {
file_reader_error(r, "Unexpected end of file.");
}
r->pos = r->line;
}
double file_reader_get_double(FileReader *r) {
r->prev = r->pos;
errno = 0;
double d = strtod(r->pos,&r->pos);
if (errno != 0)
file_reader_error_at(r,"expecting a double precision number");
return d;
}
long file_reader_get_int(FileReader *r) {
r->prev = r->pos;
errno = 0;
long d = strtol(r->pos,&r->pos,0);
if (errno != 0)
file_reader_error_at(r,"expecting an integer");
return d;
}
unsigned long file_reader_get_unsigned_int(FileReader *r) {
r->prev = r->pos;
errno = 0;
unsigned long d = strtoul(r->pos,&r->pos,0);
if (errno != 0)
file_reader_error_at(r,"expecting an unsigned integer");
return d;
}
void file_reader_get_word(FileReader *r, char *txt, int maxl) {
r->prev = r->pos;
char *next = r->pos;
while (isblank(*next)) next += 1;
int p = 0;
if (isspace(*next))
file_reader_error_at(r,"expecting a string");
while (!isspace(*next)) {
if (*next == 0)break;
txt[p++] = *next;
if(p+1 >= maxl) file_reader_error_at(r,"string is longer than expected");
next++;
}
txt[p] = '\0';
r->pos = next;
}
void file_reader_get_quoted_string(FileReader *r, char *txt, int maxl) {
r->prev = r->pos;
char *next = r->pos;
while (isblank(*next)) next += 1;
if (*next != '"') file_reader_error_at(r,"expecting a quoted string");
int p = 0;
next++;
while (*next != '"') {
if (*next == '\n' || *next == 0) file_reader_error_at(r,"expecting a quoted string");
txt[p++] = *next;
if(p+1 >= maxl) file_reader_error_at(r,"quoted string is longer than expected");
next++;
}
next++;
txt[p] = '\0';
r->pos = next;
}
void file_reader_assert(FileReader *r, const char *txt) {
r->prev = r->pos;
size_t l = strlen(txt);
if (strncmp(txt, r->pos, l) == 0) {
r->pos += l;
return;
}
file_reader_error_at(r,"expecting '%s'",txt);
}
FileReader *file_reader_new(const char *file_name) {
FileReader *r = malloc(sizeof(FileReader));
r->file_name = strdup(file_name);
r->f = fopen(file_name,"r");
if (!r->f) file_reader_error(r,"Cannot open file.");
r->iline = -1;
r->pos = NULL;
file_reader_next_line(r);
return r;
}
#ifndef FILE_READER_H
#define FILE_READER_H
typedef struct FileReaderStruct FileReader;
void file_reader_delete(FileReader *r);
void file_reader_error(FileReader *r, const char *msg, ...);
void file_reader_error_at(FileReader *r, const char *msg,...);
int file_reader_next_line_or_eof(FileReader *r);
void file_reader_next_line(FileReader *r);
void file_reader_next_line_skip(FileReader *r);
double file_reader_get_double(FileReader *r);
long file_reader_get_int(FileReader *r);
unsigned long file_reader_get_unsigned_int(FileReader *r);
void file_reader_get_word(FileReader *r, char *txt, int maxl);
void file_reader_get_quoted_string(FileReader *r, char *txt, int maxl);
void file_reader_assert(FileReader *r, const char *txt);
FileReader *file_reader_new(const char *file_name);
#endif
......@@ -1600,7 +1600,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
}
// allocate all mesh-dependant structures
static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
if(problem->mesh)
mesh_free(problem->mesh);
problem->mesh = mesh;
......@@ -1651,11 +1651,8 @@ double * fluid_problem_bulk_force(FluidProblem *problem) {
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume)
void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume)
{
struct timespec tic, toc;
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin);
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *concentration_cg = (double*)malloc(sizeof(double)*new_mesh->n_elements*N_N);
......@@ -1724,14 +1721,6 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(old_particle_uvw);
}
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {
Mesh *mesh = mesh_load_msh(mesh_file_name);
if (!mesh){
printf("cannot import mesh file '%s'\n", mesh_file_name);
}
fluid_problem_set_mesh(problem,mesh);
}
void fluid_problem_after_import(FluidProblem *problem)
{
......@@ -1983,15 +1972,6 @@ int *fluid_problem_periodic_mapping(const FluidProblem *p)
return p->mesh->periodic_mapping;
}
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals)
{
Mesh *m = malloc(sizeof(Mesh));
mesh_set_elements(m,n_nodes,x,n_elements,elements,n_boundaries,boundaries,boundary_tags,n_physicals,physicals);
mesh_gen_edges(m,n_boundaries, boundaries, boundary_tags);
fluid_problem_set_mesh(p,m);
}
int fluid_problem_mesh_n_boundaries(const FluidProblem *p) {
return p->mesh->n_boundaries;
}
......
......@@ -110,8 +110,8 @@ void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor,int temporal, int advection);
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume);
void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume);
void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem);
......@@ -132,7 +132,7 @@ int fluid_problem_n_nodes(const FluidProblem *p);
double *fluid_problem_coordinates(const FluidProblem *p);
int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals);
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals, const int *periodic);
int fluid_problem_n_mesh_boundaries(const FluidProblem *p);
void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize);
void fluid_problem_mesh_boundary_interfaces(const FluidProblem *p, int bid, int *binterfaces);
......
This diff is collapsed.
#ifndef GMSH_MESH_H
#define GMSH_MESH_H
#include <stdlib.h>
typedef struct {
int tag;
int dim;
int parent_tag;
int parent_dim;
size_t n_physicals;
int *physicals;
size_t n_partitions;
int *partitions;
size_t n_nodes;
int element_type;
double *x;
size_t *nodes_tag;
size_t n_elements;
size_t *elements; // lines in 1d, triangles in 2d
size_t *elements_tag;
size_t n_nodes_by_element;
size_t n_bounds;
double bbmin[3],bbmax[3];
int *bounds;
} Entity;
typedef struct {
int n_affine;
double *affine;
int dim;
int tag;
int parent_tag;
size_t n_nodes;
size_t *nodes; // 2*n_nodes
}Periodic;
typedef struct {
int n_partitions;
size_t n_entities[4];
Entity **entities[4];
int n_physical_names;
char **physical_names;
int *physical_name_tag;
int *physical_name_dim;
int n_periodic;
Periodic **periodic;
} GmshMesh;
GmshMesh *gmsh_mesh_new();
void gmsh_mesh_read(GmshMesh *mesh, const char *file_name);
void gmsh_mesh_write(const GmshMesh *mesh, const char *filename);
void gmsh_mesh_free(GmshMesh *mesh);
void gmsh_mesh_partition_2d(GmshMesh *mesh, int npart);
Entity *gmsh_mesh_new_entity(GmshMesh *mesh, int dim, int tag,int parent_dim,int parent_tag);
Entity *gmsh_mesh_get_entity(const GmshMesh *mesh, int dim, int tag);
void gmsh_periodic_free(Periodic *p);
GmshMesh *gmsh_mesh_load_msh(const char *filename);
int gmsh_mesh_n_physical_names(GmshMesh *mesh);
char* gmsh_mesh_physical_name(GmshMesh *mesh, int i);
int* gmsh_mesh_physical_name_tags(GmshMesh *mesh);
int* gmsh_mesh_physical_name_dims(GmshMesh *mesh);
int gmsh_mesh_n_nodes(GmshMesh *mesh, const int dimension);
int* gmsh_mesh_nodes_map(GmshMesh *mesh, const int dimension);
double* gmsh_mesh_nodes(GmshMesh *mesh, const int dimension);
size_t* gmsh_mesh_nodes_tag(GmshMesh *mesh, const int dimension);
size_t gmsh_mesh_n_entities(GmshMesh *mesh, int dimension);
int gmsh_mesh_entity_tag(GmshMesh *mesh, int dimension, int ie);
int gmsh_mesh_entity_dim(GmshMesh *mesh, int dimension, int ie);
size_t gmsh_mesh_entity_n_physicals(GmshMesh *mesh, int dim, int ie);
int* gmsh_mesh_entity_physicals(GmshMesh *mesh, int dimension, int ie);
size_t gmsh_mesh_entity_n_elements(GmshMesh *mesh, int dimension, int ie);
size_t* gmsh_mesh_entity_elements(GmshMesh *mesh, int dim, int ie);
size_t gmsh_mesh_entity_n_nodes(GmshMesh *mesh, int dim, int ie);
size_t gmsh_mesh_entity_n_nodes_by_element(GmshMesh *mesh, int dim, int ie);
size_t* gmsh_mesh_entity_nodes_tag(GmshMesh *mesh, int dim, int ie);
double* gmsh_mesh_entity_nodes(GmshMesh *mesh, int dim, int ie);
int gmsh_mesh_element_type(GmshMesh *mesh, int dimension, int ie);
size_t gmsh_mesh_element_tag(GmshMesh *mesh, int dim, int ie, int id);
size_t gmsh_mesh_element_n_partitions(GmshMesh *mesh, int dim, int ie);
int* gmsh_mesh_element_partitions(GmshMesh *mesh, int dim, int ie);
size_t gmsh_mesh_element_n_nodes(GmshMesh *mesh, int dim, int ie);
size_t* gmsh_mesh_element_nodes_tag(GmshMesh *mesh, int dim, int ie, int ielt);
int gmsh_mesh_n_periodic(GmshMesh *mesh);
int gmsh_mesh_periodic_dim(GmshMesh *mesh, int ip);
int gmsh_mesh_periodic_tag(GmshMesh *mesh, int ip);
int gmsh_mesh_periodic_parent_tag(GmshMesh *mesh, int ip);
size_t gmsh_mesh_periodic_n_nodes(GmshMesh *mesh, int ip);
size_t gmsh_mesh_periodic_node(GmshMesh *mesh, int ip, int i);
#endif
......@@ -22,7 +22,6 @@
*/
#include "mesh.h"
#include "gmsh_mesh.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
......@@ -67,7 +66,7 @@ static void sort_edge_nodes(int *n) {
#endif
}
void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements;
HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges);
HalfEdge *he = halfedges; // current
......@@ -137,169 +136,6 @@ void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tag
}
}
int get_node_id(int id, int *nodes_map, int n_nodes) {
int *found = bsearch(&id, nodes_map, n_nodes, sizeof(int), intcmp);
return found ? found-nodes_map : -1;
}
Mesh *mesh_load_msh(const char *filename)
{
GmshMesh* gmsh_m = gmsh_mesh_new();
gmsh_mesh_read(gmsh_m, filename);
Mesh *m = malloc(sizeof(Mesh));
m->interfaces = NULL;
// elements
m->n_elements = 0;
for(int ie = 0; ie < gmsh_m->n_entities[DIMENSION]; ++ie){
const Entity *e = gmsh_m->entities[DIMENSION][ie];
if (e->n_nodes_by_element != DIMENSION+1) {
exit(1);
}
m->n_elements += e->n_elements;
}
m->elements = malloc(sizeof(int)*m->n_elements*(DIMENSION+1));
int i_element = 0;
for(int ie = 0; ie < gmsh_m->n_entities[DIMENSION]; ++ie){
const Entity *e = gmsh_m->entities[DIMENSION][ie];
for(size_t i = 0; i < e->n_elements; ++i) {
for(size_t j = 0; j < e->n_nodes_by_element; ++j) {
m->elements[i_element++] = e->elements[i*e->n_nodes_by_element+j];
}
}
}
// assign id to nodes touched by elements
int n_enodes = m->n_elements*(DIMENSION+1);
int *nodes_map = malloc(sizeof(int)*n_enodes);
memcpy(nodes_map, m->elements, sizeof(int)*n_enodes);
qsort(nodes_map, n_enodes, sizeof(int), intcmp);
m->n_nodes = 0;
for (int i = 0; i < n_enodes; ++i) {
if (i == 0 || nodes_map[i] != nodes_map[m->n_nodes-1]) {
nodes_map[m->n_nodes++] = nodes_map[i];
}
}
for (int i = 0; i < n_enodes; ++i) {
m->elements[i] = get_node_id(m->elements[i], nodes_map, m->n_nodes);
}
// Nodes
m->x = malloc(sizeof(double)*3*m->n_nodes);
for (int d = 0; d < 4; ++d){
for(int ie = 0; ie < gmsh_m->n_entities[d]; ++ie){
const Entity *e = gmsh_m->entities[d][ie];
for(size_t i = 0; i < e->n_nodes; ++i){
int nid = get_node_id(e->nodes_tag[i], nodes_map, m->n_nodes);
if (nid >= 0) {
for (int id = 0; id < 3; ++id) {
m->x[nid*3+id] = e->x[i*3+id];
}
}
}
}
}
// Boundary
m->n_boundaries = 0;
int *boundary_tag_map = malloc(sizeof(int)*gmsh_m->n_physical_names);
for (int i = 0; i < gmsh_m->n_physical_names; ++i) {
if (gmsh_m->physical_name_dim[i] == DIMENSION-1) {
boundary_tag_map[i] = m->n_boundaries++;
}
else {
boundary_tag_map[i] = -1;
}
}
m->boundary_names = malloc(sizeof(char*)*m->n_boundaries);
m->n_boundaries = 0;
for (int i = 0; i < gmsh_m->n_physical_names; ++i) {
if (gmsh_m->physical_name_dim[i] == DIMENSION-1) {
m->boundary_names[m->n_boundaries++] = strdup(gmsh_m->physical_names[i]);
}
}
int n_bnd_el = 0;
for (int i = 0; i < gmsh_m->n_entities[DIMENSION-1]; ++i) {
const Entity *e = gmsh_m->entities[DIMENSION-1][i];
for (int itag = 0; itag < e->n_physicals; ++itag) {
for (int ip = 0; ip < gmsh_m->n_physical_names; ++ip) {
if (gmsh_m->physical_name_dim[ip] == e->dim && gmsh_m->physical_name_tag[ip] == e->physicals[itag]) {
for (int iel = 0; iel < e->n_elements; ++iel) {
int keep = 1;
for (int inode = 0; inode < DIMENSION; ++inode) {
int nid = get_node_id(e->elements[iel*DIMENSION+inode], nodes_map, m->n_nodes);
keep &= (nid >= 0);
}
if (keep)
n_bnd_el +=1;
}
}
}
}
}
int *bnd_elements = malloc(sizeof(int)*n_bnd_el*DIMENSION);
int *bnd_tags = malloc(sizeof(int)*n_bnd_el);
n_bnd_el = 0;
for (int ient = 0; ient < gmsh_m->n_entities[DIMENSION-1]; ++ient) {
const Entity *e = gmsh_m->entities[DIMENSION-1][ient];
for (int itag = 0; itag < e->n_physicals; ++itag) {
for (int ip = 0; ip < gmsh_m->n_physical_names; ++ip) {
if (gmsh_m->physical_name_dim[ip] == e->dim && gmsh_m->physical_name_tag[ip] == e->physicals[itag]) {
for (int iel = 0; iel < e->n_elements; ++iel) {
int nodes[DIMENSION], keep = 1;
for (int inode = 0; inode < DIMENSION; ++inode) {
nodes[inode] = get_node_id(e->elements[iel*DIMENSION+inode], nodes_map, m->n_nodes);
keep &= nodes[inode] >= 0;
}
if (keep) {
for (int inode = 0; inode < DIMENSION; ++inode) {
bnd_elements[n_bnd_el*DIMENSION+inode] = nodes[inode];
}
sort_edge_nodes(&bnd_elements[n_bnd_el*DIMENSION]);
bnd_tags[n_bnd_el] = boundary_tag_map[ip];
n_bnd_el += 1;
}
}
}
}
}
}
// Periodic
m->n_periodic = gmsh_m->n_periodic;
m->periodic_mapping = malloc(sizeof(int)*m->n_nodes);
for(int i = 0; i < m->n_nodes; ++i){
m->periodic_mapping[i] = i;
}
for(int i = 0; i < m->n_periodic; ++i){
Periodic *p = gmsh_m->periodic[i];
for(int j = 0; j < p->n_nodes; ++j){
int nid = get_node_id(p->nodes[2*j+1], nodes_map, m->n_nodes);
int nidSlave = get_node_id(p->nodes[2*j+0], nodes_map, m->n_nodes);
m->periodic_mapping[nid] = m->periodic_mapping[nidSlave];
}
}
int *renumber = malloc(sizeof(int)*m->n_nodes);
int c = 0;
for (int i = 0; i < m->n_nodes; ++i) {
if (m->periodic_mapping[i] == i)
renumber[i] = c++;
}
for (int i = 0; i < m->n_nodes; ++i) {
m->periodic_mapping[i] = renumber[m->periodic_mapping[i]];
}
free(renumber);
free(boundary_tag_map);
mesh_gen_edges(m, n_bnd_el, bnd_elements, bnd_tags);
free(bnd_elements);
free(bnd_tags);
free(nodes_map);
gmsh_mesh_free(gmsh_m);
return m;
}
void mesh_free(Mesh *m)
{
free(m->elements);
......@@ -313,137 +149,29 @@ void mesh_free(Mesh *m)
free(m);
}
int mesh_write_msh(Mesh *mesh, FILE *f)
{
fprintf(f,"$MeshFormat\n2.2 0 0\n$EndMeshFormat\n");
fprintf(f,"$Nodes\n");
fprintf(f, "%d\n", mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%i %.16g %.16g %.16g\n", i+1, mesh->x[i*3], mesh->x[i*3+1], mesh->x[i*3+2]);
}
fprintf(f, "$EndNodes\n");
fprintf(f, "$Elements\n");
fprintf(f, "%i\n", mesh->n_elements);
for (int i= 0; i < mesh->n_elements; ++i) {
#if DIMENSION==2
int *tri = mesh->elements+i*3;
fprintf(f, "%i 2 2 1 1 %i %i %i\n", i+1, tri[0]+1, tri[1]+1, tri[2]+1);
#else
int *tet = mesh->elements+i*4;
fprintf(f, "%i 4 2 1 1 %i %i %i %i\n", i+1, tet[0]+1, tet[1]+1, tet[2]+1, tet[3]+1);
#endif
}
fprintf(f, "$EndElements\n");
return 0;
}
int mesh_write_msh_scalar(const Mesh *mesh, FILE *f, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id)
{
fprintf(f,"$NodeData\n");
fprintf(f, "1\n\"%s\"\n",field_name);