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

two criteria for adapt

parent 173a8e2b
......@@ -39,8 +39,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, gradmin, gradmax, lcmin) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(lcmin))
def adapt_mesh(self, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(gradPmin), c_double(gradPmax), c_double(lcmin), c_double(n_el))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
......
......@@ -647,7 +647,7 @@ void fluid_problem_free(FluidProblem *problem)
mesh_free(problem->mesh);
}
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
{
double ngradM = 0.;
Mesh *mesh = problem->mesh;
......@@ -697,19 +697,42 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
/* ngrad2 += gradP*gradP;*/
/* }*/
double ngrad = sqrt(ngrad2);
#if DIMENSION == 2
double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
#else
double lc = lcmin /fmin(1,fmax(ngrad*ngrad/(gradmax*gradmax),gradmin*gradmin/(gradmax*gradmax)));
#endif
ngrad = pow(ngrad2, 0.25);
for (int j=0; j<DIMENSION; ++j){
double gradP = 0;
for (int k=0; k<N_SF; ++k){
gradP += dphi[k][j]*P[k];
}
ngrad2 += gradP*gradP;
}
ngrad = sqrt(ngrad2);
lc = fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc);
for (int j = 0; j < N_N; ++j){
new_mesh_size[el[j]] += lc;
num_lc[el[j]] += 1.;
}
}
double nOfEl = 0.0;
for (int i = 0; i < mesh->n_nodes; ++i){
new_mesh_size[i] /= num_lc[i];
double nOfElByNode = 4*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*sqrt(3));
nOfEl += nOfElByNode;
}
printf("Total number of element before rescale = %g\n",nOfEl);
double ratioEl = n_el/nOfEl;
for (int i = 0; i < mesh->n_nodes; ++i){
new_mesh_size[i] = fmax(new_mesh_size[i]/sqrt(ratioEl),lcmin);
}
nOfEl = 0.;
for (int i = 0; i < mesh->n_nodes; ++i){
double nOfElByNode = 4*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*sqrt(3));
nOfEl += nOfElByNode;
}
printf("Total number of element after rescale = %g\n",nOfEl);
FILE *f = fopen("adapt/lc.pos", "w");
if(!f)
printf("cannot open file adapt/lc.pos for writing\n");
......@@ -730,9 +753,9 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
{
fluid_problem_adapt_gen_mesh(problem, gradmin, gradmax, lcmin);
fluid_problem_adapt_gen_mesh(problem, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el);
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
......
......@@ -42,7 +42,7 @@ void 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);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
#endif
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