Commit 149bd673 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

new adapt2d

parent e2e2d4eb
......@@ -39,8 +39,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
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 adapt_mesh(self, lcmin, lcmax, n_el) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmin), c_double(lcmax), 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))
......
......@@ -778,19 +778,15 @@ void fluid_problem_free(FluidProblem *problem)
mesh_free(problem->mesh);
}
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmin, double lcmax, double n_el)
{
double ngradM = 0.;
Mesh *mesh = problem->mesh;
double *solution = problem->solution;
double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
double *num_lc = malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i){
new_mesh_size[i] = 0.;
num_lc[i] = 0.;
new_mesh_size[i] = 0;
}
gradmax = sqrt(gradmax);
gradmin = sqrt(gradmin);
double *porosity = problem->porosity;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
......@@ -809,58 +805,47 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
invDD(dxdxi, dxidx);
double dphi[N_SF][DIMENSION];
grad_shape_functions(dxidx, dphi);
double ngrad2 = 0;
double ngradU2 = 0;
for (int i=0; i<DIMENSION; ++i){
for (int j=0; j<DIMENSION; ++j){
double gradU = 0;
for (int k=0; k<N_SF; ++k){
gradU += dphi[k][j]*U[i][k]/C[k];
gradU += U[i][k]*dphi[k][j]/C[k];
}
ngrad2 += gradU*gradU;
ngradU2 += gradU*gradU;
}
}
/* double 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;*/
/* }*/
double ngrad = sqrt(ngrad2);
double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
double ngradU = sqrt(ngradU2);
double ngradP2 = 0;
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;
ngradP2 += gradP*gradP;
}
ngrad = sqrt(ngrad2);
lc = fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc);
double ngradP = sqrt(ngradP2);
double olc = ngradP;
for (int j = 0; j < N_N; ++j){
new_mesh_size[el[j]] += lc;
num_lc[el[j]] += 1.;
new_mesh_size[el[j]] = fmax(olc,new_mesh_size[el[j]]);
}
}
double nOfEl = 0.0;
#if DIMENSION==2
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));
double nOfElByNode = 4/sqrt(3)*problem->node_volume[i]*(new_mesh_size[i]*new_mesh_size[i]);
nOfEl += nOfElByNode;
}
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);
//new_mesh_size[i] = fmin(lcmax,fmax(sqrt(ratioEl)/(new_mesh_size[i]),lcmin));
new_mesh_size[i] = 1./fmin(1./lcmin,fmax(new_mesh_size[i]*sqrt(ratioEl),1./lcmax));//fmin(lcmax,fmax(sqrt(ratioEl)/(new_mesh_size[i]),lcmin));
}
printf("expected nEl = %i\n", nOfEl);
#else
for (int i = 0; i < mesh->n_nodes; ++i){
new_mesh_size[i] /= num_lc[i];
double nOfElByNode = 12*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*new_mesh_size[i]*sqrt(2));
nOfEl += nOfElByNode;
}
......@@ -884,14 +869,13 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double gradmin, double
fprintf(f,"};\n");
fclose(f);
free(new_mesh_size);
free(num_lc);
system("gmsh -"xstr(DIMENSION)" adapt/mesh.geo");
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double gradPmin, double gradPmax, double lcmin, double n_el)
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmin, double lcmax, double n_el)
{
fluid_problem_adapt_gen_mesh(problem, gradmin, gradmax, gradPmin, gradPmax, lcmin, n_el);
fluid_problem_adapt_gen_mesh(problem,lcmin, lcmax,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);
......
......@@ -43,7 +43,7 @@ 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);
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 gradPmin, double gradPmax, double lcmin, double n_el);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmin, double lcmax, double n_el);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
#endif
......@@ -48,7 +48,7 @@ tEnd = 100
lcmin = 0.0001 # approx r*100 but should match the mesh size
dt = 5e-4
alpha = 2.5e-2
epsilon = alpha*lcmin**2 /nu
epsilon = alpha*lcmin**2/nu/5*10
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
......@@ -74,8 +74,8 @@ ii = 0
tic = time.clock()
forces = g*p.mass()
while t < tEnd :
if (ii%20==0 and ii != 0):
fluid.adapt_mesh(0.01,200,1e-4)
if ((ii%20==0 and ii != 0) or ii == 2):
fluid.adapt_mesh(5e-5,2e-3,60000)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
......@@ -89,8 +89,6 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
p.write(outputdir+"/part-%05d" % ioutput)
p.write_boundary_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
L = 0.025;
H = 0.5;
y = 0;
lc = 0.001;
lc = 0.0001;
Point(1) = {-L, H, 0};
Point(2) = {-L, 0, 0};
......@@ -11,8 +11,8 @@ Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Point(5) = {0,0.4,0};
Point(6) = {0,H,0};
Point(5) = {0,0.448,0};
Point(6) = {0,0.452,0};
Line(5) = {5,6};
Line Loop(1) = {1:4};
......@@ -30,8 +30,8 @@ Field[1].NNodesByEdge = 200;
Field[2] = Threshold;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.008;
Field[2].LcMax = 0.05;
Field[2].LcMin = 0.001;
Field[2].LcMax = 0.01;
Field[2].LcMin = 0.0001;
Field[2].IField = 1;
Background Field = 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