Commit cfe70a73 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

relaod bnd tags in fluids

parent ab04e7a2
......@@ -1187,31 +1187,32 @@ int particleProblemImportVtk(ParticleProblem *p, const char *filename)
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name,"Points")){
mb_read_vtk_data_array(r,&a);
size_t n_particles = a.n/3;
double *position = mb_vtk_data_array_to_double(&a,n_particles,3);
for (size_t i = 0; i < n_particles; ++i)
for (int j = 0; j < DIMENSION; ++j)
p->position[DIMENSION*i+j] = ((double*)a.data)[i*3+j];
free(a.name);
free(a.data);
p->position[DIMENSION*i+j] = position[i*3+j];
free(position);
mb_vtk_data_array_destroy(&a);
}
else if (!strcasecmp(e.name,"PointData")) {
while(mb_read_vtk_data_array(r,&a) != 0){
double *value = mb_vtk_data_array_to_double(&a,n_particles,-1);
if(!strcasecmp(a.name,"Velocity")) {
for (size_t i = 0; i < n_particles; ++i)
for (int j = 0; j < DIMENSION; ++j)
p->velocity[DIMENSION*i+j] = ((double*)a.data)[i*3+j];
p->velocity[DIMENSION*i+j] = value[i*3+j];
}
else if(!strcasecmp(a.name,"Mass")) {
for (size_t i = 0; i < n_particles; ++i)
p->particles[i].m = ((double*)a.data)[i];
p->particles[i].m = value[i];
}
else if(!strcasecmp(a.name,"Radius")) {
for (size_t i = 0; i < n_particles; ++i){
p->particles[i].r = ((double*)a.data)[i];
p->particles[i].r = value[i];
}
}
free(a.name);
free(a.data);
free(value);
mb_vtk_data_array_destroy(&a);
}
}
else {
......@@ -1238,11 +1239,17 @@ int particleProblemImportBoundaryVtk(ParticleProblem *p, const char *filename)
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name, "FieldData")) {
while(mb_read_vtk_data_array(r,&a)) {
free(a.name);
free(a.data);
mb_vtk_data_array_destroy(&a);
}
}
else if(!strcasecmp(e.name,"Piece")) {
int n_points, n_cells;
for(size_t i = 0; i < e.n_attr; ++i) {
if (!strcasecmp("NumberOfPoints",e.attr[i].name))
sscanf(e.attr[i].value,"%i",&n_points);
if (!strcasecmp("NumberOfCells",e.attr[i].name))
sscanf(e.attr[i].value,"%i",&n_cells);
}
MbXmlElement f;
for (size_t i = 0; i<e.n_attr; ++i) {
MBVtkDataArray points={0}, connectivity={0}, offsets={0}, types={0};
......@@ -1259,8 +1266,7 @@ int particleProblemImportBoundaryVtk(ParticleProblem *p, const char *filename)
offsets = a;
else {
printf("unexpected data array in \"%s\" : \"%s\"\n",filename, a.name);
if(a.name) free(a.name);
if(a.data) free(a.data);
mb_vtk_data_array_destroy(&a);
}
}
vectorClear(p->disks);
......@@ -1268,11 +1274,17 @@ int particleProblemImportBoundaryVtk(ParticleProblem *p, const char *filename)
#if DIMENSION == 3
vectorClear(p->triangles);
#endif
double *pts = points.data;
int *off = ((int*)offsets.data);
for(size_t i = 0; i < types.n; ++i) {
int t = ((int*)types.data)[i];
int *nodes = ((int*)connectivity.data)+(i==0 ? 0 : off[i-1]);
double *pts = mb_vtk_data_array_to_double(&points,n_points,3);
int *off = mb_vtk_data_array_to_int(&offsets,n_cells,1);
int *con = mb_vtk_data_array_to_int(&connectivity,off[n_cells-1],1);
int *typ = mb_vtk_data_array_to_int(&types,n_cells,1);
mb_vtk_data_array_destroy(&points);
mb_vtk_data_array_destroy(&offsets);
mb_vtk_data_array_destroy(&connectivity);
mb_vtk_data_array_destroy(&types);
for(size_t i = 0; i < n_cells; ++i) {
int t = typ[i];
int *nodes = con+(i==0 ? 0 : off[i-1]);
if (t == 1){
particleProblemAddBoundaryDisk(p,pts+nodes[0]*3,0.,"");
}
......@@ -1285,14 +1297,6 @@ int particleProblemImportBoundaryVtk(ParticleProblem *p, const char *filename)
else
printf("Unknown boundary object : %i\n", t);
}
free(points.name);
free(points.data);
free(connectivity.name);
free(connectivity.data);
free(offsets.name);
free(offsets.data);
free(types.name);
free(types.data);
}
else {
printf("unexpected element in \"%s\" : \"%s\"\n",filename, f.name);
......
......@@ -135,27 +135,53 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
return 0;
}
static void read_vtk_piece(MbXmlReader *r, FluidProblem *problem) {
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;
Mesh *m = problem->mesh;
while(mb_xml_read_element(r,&e)!=1) {
if (!strcasecmp(e.name,"Points")){
mb_read_vtk_data_array(r,&a);
if(a.name)free(a.name);
m->n_nodes = a.n/3;
if(m->x) free(m->x);
m->x = a.data;
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")) {
m->n_elements = a.n/(DIMENSION+1);
m->elements = a.data;
} else {
free(a.data);
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);
}
if(a.name)free(a.name);
}
}
else if (!strcasecmp(e.name,"PointData")) {
......@@ -166,20 +192,19 @@ static void read_vtk_piece(MbXmlReader *r, FluidProblem *problem) {
if(!strcasecmp(a.name,"Porosity")) {
if(problem->porosity)
free(problem->porosity);
problem->porosity = a.data;
problem->porosity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
}
else if(!strcasecmp(a.name,"Pressure")) {
memcpy(problem->solution+DIMENSION*m->n_nodes, a.data,sizeof(double)*m->n_nodes);
free(a.data);
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")) {
memcpy(problem->solution, a.data,sizeof(double)*m->n_nodes*DIMENSION);
free(a.data);
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);
}
else {
free(a.data);
}
if(a.name)free(a.name);
mb_vtk_data_array_destroy(&a);
}
}
else {
......@@ -188,17 +213,42 @@ static void read_vtk_piece(MbXmlReader *r, FluidProblem *problem) {
}
mb_xml_end_element(r,&e);
}
}
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);
read_vtk_piece(r,problem);
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);
......@@ -251,7 +301,6 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
#endif
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) {
......@@ -273,52 +322,17 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
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");
//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;
}
......
......@@ -26,26 +26,6 @@ char *mb_xml_start_element(MbXmlReader *reader)
return NULL;
}
char *mb_xml_read_content(MbXmlReader *reader, size_t *size) {
*size = 0;
for(int doc=getc(reader->f); doc!=EOF;doc=getc(reader->f)) {
yxml_ret_t r = yxml_parse(&reader->x,doc);
switch(r) {
case YXML_CONTENT:
if ((*size)+1 >= reader->content_buf_size){
reader->content_buf_size*=2;
reader->content_buf = realloc(reader->content_buf,sizeof(char)*reader->content_buf_size);
}
reader->content_buf[(*size)++] = *reader->x.data;
break;
case YXML_OK: case YXML_PISTART: case YXML_PIEND: case YXML_PICONTENT: break;
case YXML_ELEMEND: reader->content_buf[(*size)] = '\0';at_end = 1;return reader->content_buf;
default : printf("unknown xml in read content %i\n",r); exit(1); }
}
return NULL;
}
int mb_xml_read_attribute(MbXmlReader *reader, char *value,size_t bufsize)
{
size_t pos = 0;
......@@ -122,6 +102,24 @@ int mb_xml_read_element(MbXmlReader *reader, MbXmlElement *e) {
return 0;
}
static char *mb_xml_read_content(MbXmlReader *reader) {
size_t size = 0;
for(int doc=getc(reader->f); doc!=EOF;doc=getc(reader->f)) {
yxml_ret_t r = yxml_parse(&reader->x,doc);
switch(r) {
case YXML_CONTENT:
if (size+1 >= reader->content_buf_size){
reader->content_buf_size*=2;
reader->content_buf = realloc(reader->content_buf,sizeof(char)*reader->content_buf_size);
}
reader->content_buf[size++] = *reader->x.data;
break;
case YXML_OK: case YXML_PISTART: case YXML_PIEND: case YXML_PICONTENT: break;
case YXML_ELEMEND: reader->content_buf[size] = '\0';at_end = 1;return reader->content_buf;
default : printf("unknown xml in read content %i\n",r); exit(1); }
}
return NULL;
}
int mb_read_vtk_data_array(MbXmlReader *reader, MBVtkDataArray *a) {
MbXmlElement e;
......@@ -133,48 +131,25 @@ int mb_read_vtk_data_array(MbXmlReader *reader, MBVtkDataArray *a) {
exit(1);
}
a->name = NULL;
a->data = NULL;
a->type = 0;
int format = 0;
a->content = NULL;
a->type = NULL;
a->format = 0;
a->number_of_components = -1;
a->number_of_tuples = -1;
for (size_t i = 0; i < e.n_attr; ++i) {
if(!strcasecmp(e.attr[i].name,"type")) {
if(!strcasecmp(e.attr[i].value,"Float64"))
a->type = 1;
if(!strcasecmp(e.attr[i].value,"Float32"))
a->type = 1;
if(!strcasecmp(e.attr[i].value,"Int32"))
a->type = 2;
}
if(!strcasecmp(e.attr[i].name,"format") && !strcasecmp(e.attr[i].value,"binary"))
format = 1;
if(!strcasecmp(e.attr[i].name,"Name"))
if(!strcasecmp(e.attr[i].name,"type"))
a->type=strdup(e.attr[i].value);
else if(!strcasecmp(e.attr[i].name,"format") && !strcasecmp(e.attr[i].value,"binary"))
a->format = 1;
else if(!strcasecmp(e.attr[i].name,"Name"))
a->name = strdup(e.attr[i].value);
}
size_t size;
char *content = mb_xml_read_content(reader,&size);
if (a->type == 2 && format == 0) {
FILE *f = fmemopen(content,size,"r");
a->n = 0;
int buf;
while(fscanf(f,"%i",&buf) == 1) a->n++;
fseek(f,0,SEEK_SET);
a->data = malloc(sizeof(int)*a->n);
for (size_t p = 0;p < a->n;p++)
fscanf(f,"%i",&((int*)a->data)[p]);
fclose(f);
}
if (a->type == 1 && format == 0) {
FILE *f = fmemopen(content,size,"r");
a->n = 0;
double buf;
while(fscanf(f,"%le",&buf) == 1) a->n++;
fseek(f,0,SEEK_SET);
a->data = malloc(sizeof(double)*a->n);
for (size_t p = 0;p < a->n;p++) {
fscanf(f,"%le",&((double*)a->data)[p]);
else if(!strcasecmp(e.attr[i].name,"NumberOfTuples"))
sscanf(e.attr[i].value,"%i",&a->number_of_tuples);
else if(!strcasecmp(e.attr[i].name,"NumberOfComponents")){
sscanf(e.attr[i].value,"%i",&a->number_of_components);
}
fclose(f);
}
a->content = strdup(mb_xml_read_content(reader));
mb_xml_end_element(reader,&e);
return 1;
}
......@@ -202,3 +177,65 @@ void mb_xml_reader_destroy(MbXmlReader *xml){
fclose(xml->f);
free(xml);
}
static void validate_array(MBVtkDataArray *a, int number_of_tuples, int number_of_components){
if (number_of_tuples !=-1) {
if(number_of_tuples != a->number_of_tuples && a->number_of_tuples != -1){
printf("invalid number of tuples in vtk array \"%s\"\n", a->name);
exit(1);
}
a->number_of_tuples = number_of_tuples;
}
if (a->number_of_tuples == -1) {
printf("number of tuples not specified for vtk array \"%s\"\n", a->name);
}
if (number_of_components !=-1) {
if(number_of_components != a->number_of_components && a->number_of_components != -1){
printf("invalid number of components in vtk array \"%s\"\n", a->name);
exit(1);
}
a->number_of_components = number_of_components;
}
if (a->number_of_components == -1) {
printf("number of components not specified for vtk array \"%s\"", a->name);
}
}
double *mb_vtk_data_array_to_double(MBVtkDataArray *a, int number_of_tuples, int number_of_components)
{
validate_array(a, number_of_tuples, number_of_components);
size_t size = a->number_of_tuples*a->number_of_components;
double *v = malloc(size*sizeof(double));
FILE *f = fmemopen(a->content, strlen(a->content),"r");
for(size_t i = 0; i < size; ++i){
if(fscanf(f,"%le",&v[i]) != 1) {
printf("cannot parse double in vtk array \"%s\"", a->name);
exit(1);
}
}
fclose(f);
return v;
}
int *mb_vtk_data_array_to_int(MBVtkDataArray *a, int number_of_tuples, int number_of_components)
{
validate_array(a, number_of_tuples, number_of_components);
size_t size = a->number_of_tuples*a->number_of_components;
int *v = malloc(size*sizeof(int));
FILE *f = fmemopen(a->content, strlen(a->content),"r");
for(size_t i = 0; i < size; ++i){
if(fscanf(f,"%d",&v[i]) != 1) {
printf("cannot parse double in vtk array \"%s\"", a->name);
exit(1);
}
}
fclose(f);
return v;
}
void mb_vtk_data_array_destroy(MBVtkDataArray *a)
{
free(a->content);
free(a->name);
free(a->type);
}
......@@ -3,8 +3,11 @@
#include <stdlib.h>
typedef struct {
char *name;
void *data;
int type; //1 double, 2 int, 0 unknown
char *content;
char *type;
int format;
int number_of_tuples;
int number_of_components;
size_t n;
} MBVtkDataArray;
......@@ -26,4 +29,7 @@ void mb_xml_end_element(MbXmlReader *reader,MbXmlElement *e);
int mb_read_vtk_data_array(MbXmlReader *reader, MBVtkDataArray *a);
MbXmlReader *mb_xml_reader_create(const char *filename);
void mb_xml_reader_destroy(MbXmlReader *xml);
double *mb_vtk_data_array_to_double(MBVtkDataArray *a, int number_of_tuples, int number_of_components);
int *mb_vtk_data_array_to_int(MBVtkDataArray *a, int number_of_tuples, int number_of_components);
void mb_vtk_data_array_destroy(MBVtkDataArray *a);
#endif
......@@ -5,12 +5,12 @@ y = 0;
lc = 0.01;
Point(1) = {-L, H, -P, lc};
Point(2) = {-L, -H, -P, lc};
Point(3) = {L, -H, -P, lc};
Point(2) = {-L, 0, -P, lc};
Point(3) = {L, 0, -P, lc};
Point(4) = {L, H, -P, lc};
Point(5) = {-L, H, P, lc};
Point(6) = {-L, -H, P, lc};
Point(7) = {L, -H, P, lc};
Point(6) = {-L, 0, P, lc};
Point(7) = {L, 0, P, lc};
Point(8) = {L, H, P, lc};
Line(1) = {1, 2};
......
Supports Markdown
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