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

bug in 2 fluids, reload weak boundary vtk, fluid constructor

parent df587e97
Pipeline #4065 passed with stage
in 1 minute and 28 seconds
......@@ -36,7 +36,7 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, g, mu, rho, strong_boundaries, tf=0):
def __init__(self, g, mu, rho, strong_boundaries):
nsb = len(strong_boundaries)
class Bnd :
def __init__(self, b) :
......@@ -56,7 +56,8 @@ class fluid_problem :
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2],BNDCB(Bnd(i[3]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), np2c(mu), np2c(rho), nsb, asb, c_int(tf)))
n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
self._fp = c_void_p(lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......@@ -64,6 +65,9 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def load_msh(self, mesh_file_name) :
lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
def set_weak_boundary(self, tag, bndtype) :
lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bndtype.encode())
......
......@@ -131,7 +131,7 @@ static double invDD(const double m[3][3], double inv[3][3]) {
#define U 0
#define A (D+1)
int fluid_problem_n_fields(const FluidProblem *p) {if(!p->tf) return (D+1);else return (D+2);}
int fluid_problem_n_fields(const FluidProblem *p) {return D + p->n_fluids;}
#if D==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
......@@ -175,7 +175,7 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
double Ra, a, mu, rho;
//printf("c=%e, dc[0]=%e, dc[1]=%e, cold=%e, dt=%e\n",c,dc[0],dc[1],cold,dt);
if (!problem->tf){
if (problem->n_fluids == 1){
rho = problem->rho[0];
mu = problem->mu[0];
}
......@@ -195,7 +195,7 @@ static void node_force_f(FluidProblem *problem, double *f0, double *sol, double
double vol = problem->particle_volume[ip];
double g = problem->g;
double gamma;
if(problem->tf){
if(problem->n_fluids == 2){
double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
gamma = gamma0*a + gamma1*(1-a);
......@@ -346,7 +346,7 @@ static void f_outflow(FluidProblem *problem,const double *n, double *f,const dou
}
}
f[P] = un;
if(!problem->tf){
if(problem->n_fluids == 1){
mu = problem->mu[0];
rho = problem->rho[0];
}
......@@ -387,7 +387,7 @@ static void f_inflow(FluidProblem *problem,const double *n, double *f,const doub
}
f[P] = un;
if(!problem->tf){
if(problem->n_fluids == 1){
mu = problem->mu[0];
rho = problem->rho[0];
}
......@@ -425,7 +425,7 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
f[U+id] = p*n[id];
}
if(problem->tf){
if(problem->n_fluids == 2){
double a = s[A];
f[A] = 0;
}
......@@ -434,6 +434,12 @@ static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,con
static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
const Mesh *mesh = problem->mesh;
problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
problem->taua = realloc(problem->taua,sizeof(double)*mesh->n_elements);
problem->extraa = realloc(problem->extraa,sizeof(double)*mesh->n_elements);
double maxtaup = 0;
double mintaup = DBL_MAX;
const int n_fields = fluid_problem_n_fields(problem);
......@@ -448,14 +454,15 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
for (int i=0; i< N_N; ++i) {
double normup = 0;
c += problem->porosity[el[i]]/N_N;
if(problem->tf){
if(problem->n_fluids == 2){
a += problem->solution[el[i]*n_fields+A];
amax = fmax(amax,problem->solution[el[i]*n_fields+A]);
amin = fmin(amin,problem->solution[el[i]*n_fields+A]);
}
for(int j = 0; j < DIMENSION; ++j)
for(int j = 0; j < DIMENSION; ++j){
normup += pow2(problem->solution[el[i]*n_fields+j]);
}
normu = fmax(normu,sqrt(normup));
}
for (int i=0; i<DIMENSION; ++i){
......@@ -468,7 +475,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
const double h = 2*cbrt(detDD(dxdxi)/(8*M_PI));
#endif
double rho, mu;
if(!problem->tf){
if(problem->n_fluids == 1){
rho = problem->rho[0];
mu = problem->mu[0];
}else{
......@@ -499,7 +506,7 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
double Ra, a, aold, mu, rho, rhoold, drhodt;
if (!problem->tf){
if (problem->n_fluids == 1){
rho = problem->rho[0];
mu = problem->mu[0];
}
......@@ -525,9 +532,15 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
tau[i][j] = du[i][j]-u[i]*dc[j]/c;
utau[i] += u[j]*tau[i][j];
}
if(problem->tf){
if(problem->n_fluids == 2){
da[i] = dsol[A*D+i];
Ra += u[i]*da[i];
}
}
if(problem->n_fluids == 2){
Ra += a*divu;
drhodt = (rho-rhoold)/dt;
for (int i = 0; i < D; ++i) {
for (int j=0;j<D;++j){
uugradrho[i] += u[i]*u[j]/c *(problem->rho[0]*da[j] - problem->rho[1]*da[j]);
}
......@@ -536,10 +549,8 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
f0[P] = (c-cold)/dt;
if(problem->tf){
Ra += a*divu;
if(problem->n_fluids == 2){
f0[A] = (a*c-aold*cold)/dt;
drhodt = (rho-rhoold)/dt;
}
for (int i = 0; i < D; ++i) {
......@@ -547,16 +558,16 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==(D-1) ? -c*problem->g : 0);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==(D-1) ? -c*rho*problem->g : 0);
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup;
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup;
}
f1[(U+i)*D+i] += divu*tauc*rho;
f1[P*D+i] = -u[i] + taup*R;
if(problem->tf){
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
}
}
if(problem->tf){
if(problem->n_fluids == 2){
f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
}
}
......@@ -833,7 +844,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
break;
}
if (iphys == mesh->n_phys)
printf("Boundary tag \"%s\" not found.", bnd->tag);
printf("Boundary tag \"%s\" not found.\n", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
}
......@@ -860,12 +871,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
/*for (int i = 0; i < local_size; ++i)
printf("%8.2g ", local_vector[i]);
printf("\n");*/
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
//exit(0);
free(forced_row);
free(all_local_matrix);
free(all_local_vector);
......@@ -881,7 +888,7 @@ static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_
break;
}
if (iphys == mesh->n_phys){
printf("Boundary tag \"%s\" not found.", bnd->tag);
printf("Boundary tag \"%s\" not found.\n", bnd->tag);
}
double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
......@@ -996,8 +1003,11 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries, int two_fluids)
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
if (n_fluids != 1 && n_fluids != 2) {
printf("error : only 1 or 2 fluids are supported\n");
}
static int initialized = 0;
if (!initialized) {
hxtInitializeLinearSystems(0, NULL);
......@@ -1007,29 +1017,17 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
#endif
}
FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
Mesh *mesh = mesh_load_msh(mesh_file_name);
if (!mesh)
return NULL;
problem->mesh_tree = mesh_tree_create(mesh);
problem->tf = two_fluids;
problem->n_fluids = n_fluids;
int n_fields = fluid_problem_n_fields(problem);
problem->weak_boundaries = NULL;
problem->n_weak_boundaries = 0;
if(!problem->tf){
problem->mu = (double*)malloc(sizeof(double));
problem->rho = (double*)malloc(sizeof(double));
problem->mu[0] = mu[0];
problem->rho[0] = rho[0];
}else{
problem->mu = (double*)malloc(sizeof(double)*2);
problem->rho = (double*)malloc(sizeof(double)*2);
for (int i = 0; i < 2; ++i){
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
}
problem->mu = (double*)malloc(sizeof(double)*problem->n_fluids);
problem->rho = (double*)malloc(sizeof(double)*problem->n_fluids);
for (int i = 0; i < problem->n_fluids; ++i){
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
}
problem->g = g;
problem->mesh = mesh;
problem->strong_boundaries = (StrongBoundary*)malloc(sizeof(StrongBoundary)*n_strong_boundaries);
problem->n_strong_boundaries = n_strong_boundaries;
for (int i = 0; i < n_strong_boundaries; ++i) {
......@@ -1038,28 +1036,17 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->strong_boundaries[i].field = strong_boundaries[i].field;
problem->strong_boundaries[i].row = strong_boundaries[i].row;
}
problem->taup = (double*)malloc(sizeof(double)*mesh->n_elements);
problem->tauc = (double*)malloc(sizeof(double)*mesh->n_elements);
problem->taua = (double*)malloc(sizeof(double)*mesh->n_elements);
problem->extraa = (double*)malloc(sizeof(double)*mesh->n_elements);
problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
for(int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
problem->node_volume = NULL;
fluid_problem_compute_node_volume(problem);
problem->mesh_tree = NULL;
problem->mesh = NULL;
problem->taup = NULL;
problem->tauc = NULL;
problem->taua = NULL;
problem->extraa = NULL;
problem->porosity = NULL;
problem->oldporosity = NULL;
problem->solution = NULL;
problem->linear_system = NULL;
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
problem->n_particles = 0;
problem->particle_uvw = NULL;
problem->particle_element_id = NULL;
......@@ -1125,7 +1112,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double C[N_N],Un[D][N_N],An[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
if(problem->tf) An[i] = solution[el[i]*n_fields+D+1];
if(problem->n_fluids ==2) An[i] = solution[el[i]*n_fields+D+1];
for (int j=0; j<D; ++j) {
Un[j][i] = solution[el[i]*n_fields+j];
}
......@@ -1156,7 +1143,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
volT += vol;
gradmax = fmax(gradmax,ngrad2);
gradmin = fmin(gradmin,ngrad2);
if(problem->tf){
if(problem->n_fluids == 2){
double ngradA = 0;
for (int j=0; j<D; ++j){
double gradA = 0;
......@@ -1213,7 +1200,7 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double ngrad = sqrt(ngrad2)/gradM;
double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
if(problem->tf){
if(problem->n_fluids == 2){
double ngradA2 = 0;
//printf("2fluids\n");
for (int j=0; j<D; ++j){
......@@ -1279,6 +1266,38 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
system("gmsh -" xstr(D) " adapt/mesh.geo");
}
// allocate all mesh-dependant structures
static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
if(problem->mesh)
mesh_free(problem->mesh);
problem->mesh = mesh;
if(problem->mesh_tree)
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(mesh);
free(problem->porosity);
problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
for(int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
free(problem->solution);
int n_fields = fluid_problem_n_fields(problem);
problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
if (problem->linear_system)
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
fluid_problem_compute_node_volume(problem);
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
{
struct timespec tic, toc;
......@@ -1324,45 +1343,28 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(new_eid);
free(new_xi);
free(newx);
fluid_problem_set_mesh(problem,new_mesh);
free(problem->solution);
free(problem->taup);
free(problem->tauc);
free(problem->taua);
free(problem->extraa);
problem->taup = malloc(sizeof(double)*new_mesh->n_elements);
problem->tauc = malloc(sizeof(double)*new_mesh->n_elements);
problem->taua = malloc(sizeof(double)*new_mesh->n_elements);
problem->extraa = malloc(sizeof(double)*new_mesh->n_elements);
problem->solution = new_solution;
free(problem->porosity);
problem->porosity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(sizeof(double)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
{
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
mesh_free(problem->mesh);
problem->mesh = new_mesh;
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
for (int i = 0; i < problem->n_particles; ++i)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
#endif
}
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)
{
mesh_tree_free(problem->mesh_tree);
if(problem->mesh_tree)
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
int n_fields = fluid_problem_n_fields(problem);
for (int i = 0; i < problem->n_particles; ++i)
......@@ -1370,7 +1372,8 @@ void fluid_problem_after_import(FluidProblem *problem)
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
fluid_problem_compute_porosity(problem);
hxtLinearSystemDelete(&problem->linear_system);
if(problem->linear_system)
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
#else
......
......@@ -72,7 +72,7 @@ struct FluidProblem {
double *particle_velocity;
double *particle_mass;
double *particle_force;
int tf;
int n_fluids;
//stabilisation coefficients
double *taup;
double *tauc;
......@@ -89,7 +89,8 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries, int n_fluids);
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries);
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);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
......
......@@ -29,13 +29,14 @@
int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
{
const int n_fields = fluid_problem_n_fields(problem);
printf("n_fields = %d\n",n_fields);
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);
if (problem->mesh)
mesh_free(problem->mesh);
problem->mesh = malloc(sizeof(Mesh));
Mesh *m = problem->mesh;
m->n_nodes = -1, m->n_elements = -1;
for(size_t i = 0; i < pieceelement.n_attr; ++i) {
......@@ -52,7 +53,6 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
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);
}
......@@ -60,7 +60,6 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
int ncon = -1;
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Connectivity")) {
free(m->elements);
m->elements = mb_vtk_data_array_to_int(&a,m->n_elements,DIMENSION+1);
mb_vtk_data_array_destroy(&a);
}
......@@ -82,8 +81,8 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
else if (!strcasecmp(e.name,"PointData")) {
if (problem->solution){
free(problem->solution);
problem->solution = (double*)malloc(sizeof(double)*n_fields*m->n_nodes);
}
}
problem->solution = (double*)malloc(sizeof(double)*n_fields*m->n_nodes);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Porosity")) {
......@@ -126,6 +125,9 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
int n_phys=0;
int n_phys_alloc = 10;
MBVtkDataArray *phys =(MBVtkDataArray*) malloc(sizeof(MBVtkDataArray)*10);
int n_boundaries = 0;
int n_boundaries_alloc = 10;
MBVtkDataArray *boundaries =(MBVtkDataArray*) 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) {
......@@ -134,20 +136,29 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
}
phys[n_phys++] = a;
}
else if(!strncasecmp(a.name,"Boundary ",9)) {
if (n_boundaries == n_boundaries_alloc) {
n_boundaries_alloc *= 2;
boundaries = (MBVtkDataArray*)realloc(boundaries,n_boundaries_alloc*sizeof(MBVtkDataArray));
}
boundaries[n_boundaries++] = a;
}
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 = malloc(n_phys*sizeof(int));
m->phys_dimension = (int*)malloc(n_phys*sizeof(int));
m->phys_n_nodes = (int*)malloc(n_phys*sizeof(int));
m->phys_name = (char**)malloc(n_phys*sizeof(char*));
m->phys_nodes = (int**)malloc(n_phys*sizeof(int*));
m->phys_n_boundaries = (int*)malloc(n_phys*sizeof(int));
m->phys_boundaries = (int**)malloc(n_phys*sizeof(int*));
if(n_phys != n_boundaries) {
printf("import vtk errors : number of physical and boundaries data does not match");
}
m->phys_id = (int*)realloc(m->phys_id,n_phys*sizeof(int));
m->phys_dimension = (int*)realloc(m->phys_dimension,n_phys*sizeof(int));
m->phys_n_nodes = (int*)realloc(m->phys_n_nodes,n_phys*sizeof(int));
m->phys_name = (char**)realloc(m->phys_name,n_phys*sizeof(char*));
m->phys_nodes = (int**)realloc(m->phys_nodes,n_phys*sizeof(int*));
m->n_phys = n_phys;
for(int i = 0; i < n_phys; ++i) {
MBVtkDataArray a = phys[i];
int p, n;
......@@ -158,6 +169,15 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
mb_vtk_data_array_destroy(&phys[i]);
}
free(phys);
for(int i = 0; i < n_boundaries; ++i) {
MBVtkDataArray a = boundaries[i];
int p, n;
sscanf(a.name+9,"%i %i%n",&m->phys_dimension[i],&m->phys_id[i],&p);
m->phys_boundaries[i] = mb_vtk_data_array_to_int(&boundaries[i],-1,-1);
m->phys_n_boundaries[i] = boundaries[i].number_of_tuples;
mb_vtk_data_array_destroy(&boundaries[i]);
}
free(boundaries);
mb_xml_end_element(r,&fielddataelement);
mb_xml_end_element(r,&gridelement);
mb_xml_end_element(r,&fileelement);
......@@ -228,7 +248,7 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
fprintf(f, "%g ", problem->solution[i*n_fields+(DIMENSION)]);
}
fprintf(f,"\n</DataArray>\n");
if(problem->tf){
if(problem->n_fluids == 2){
fprintf(f,"<DataArray Name=\"Concentration\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*n_fields+(DIMENSION+1)]);
......@@ -252,6 +272,11 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
fprintf(f, "%i ", mesh->phys_nodes[i][j]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Boundary %i %i %s\" NumberOfTuples=\"%i\" NumberOfComponents=\"2\" type=\"Int32\" format=\"ascii\">\n",mesh->phys_dimension[i],mesh->phys_id[i],mesh->phys_name[i],mesh->phys_n_boundaries[i]);
for (int j = 0; j < mesh->phys_n_boundaries[i]; ++j) {
fprintf(f, "%i %i ", mesh->phys_boundaries[i][j*2+0], mesh->phys_boundaries[i][j*2+1]);
}
fprintf(f,"\n</DataArray>\n");
}
fprintf(f,"</FieldData>\n");
fprintf(f,"</UnstructuredGrid>\n");
......
......@@ -86,7 +86,8 @@ def outerBndV(x) :
strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
# ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)
]#,("Lateral",0,0,0),("Bottom",1,1,0)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries,1)
fluid = fluid.fluid_problem(g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries)
fluid.load_msh("mesh.msh")
fluid.set_weak_boundary("Bottom","Wall")
fluid.set_weak_boundary("Lateral","Wall")
fluid.set_weak_boundary("Top","Outflow")
......