Commit 3bd4504a authored by Margaux Boxho's avatar Margaux Boxho
Browse files

Final version of refinement strategy

parent 54b3fa70
Pipeline #5460 passed with stage
in 27 seconds
...@@ -1809,6 +1809,7 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double ...@@ -1809,6 +1809,7 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
printf("==> INITIALIZED THE TREE STRUCTURE \n"); printf("==> INITIALIZED THE TREE STRUCTURE \n");
problem->Tree = (Tree_cell*) malloc(mesh->n_elements*sizeof(Tree_cell)); problem->Tree = (Tree_cell*) malloc(mesh->n_elements*sizeof(Tree_cell));
problem->size_initial_mesh = mesh->n_elements; problem->size_initial_mesh = mesh->n_elements;
problem->size_initial_node = mesh->n_nodes;
deep_tree_initialized(mesh, problem->Tree, mesh->n_elements); deep_tree_initialized(mesh, problem->Tree, mesh->n_elements);
// INITIALISE THE MAPPING // INITIALISE THE MAPPING
...@@ -1820,23 +1821,89 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double ...@@ -1820,23 +1821,89 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
int n_nodes_new = 0; int n_nodes_new = 0;
int n_elements_new = 0; int n_elements_new = 0;
int n_boundaries_new = 0; int n_boundaries_new = 0;
deep_tree_refinement(mesh, problem->Tree, problem->map, problem->size_initial_mesh, lcmin, &n_nodes_new, &n_elements_new, &n_boundaries_new); int n_nodes_removed = 0;
int n_elements_removed = 0;
// Let's create the correspondence table for the removed node
int len_node = mesh->n_nodes - problem->size_initial_node;
int *tab_node = (int*) malloc(len_node*sizeof(int));
for (int ii=0; ii<len_node; ii++) {
tab_node[ii]=-1;
}
// Let's create the correpondence table for the removed element
int len_elt = mesh->n_elements - problem->size_initial_mesh;
int *tab_elt = (int*) malloc(len_elt*sizeof(int));
for (int ii=0; ii<len_elt; ii++) {
tab_elt[ii]=-1;
}
deep_tree_refinement(mesh, problem->Tree, problem->map, problem->size_initial_mesh, lcmin, &n_nodes_new, &n_elements_new, &n_boundaries_new, len_node, tab_node, problem->size_initial_node, len_elt, tab_elt, problem->size_initial_mesh, &n_nodes_removed, &n_elements_removed);
printf("There are %d nodes, %d elements in the initial mesh\n", problem->size_initial_node, problem->size_initial_mesh);
printf("There are %d nodes, %d elements in the old mesh\n", mesh->n_nodes, mesh->n_elements); printf("There are %d nodes, %d elements in the old mesh\n", mesh->n_nodes, mesh->n_elements);
printf("There are %d nodes, %d elements and %d boundaries in the mesh\n", n_nodes_new, n_elements_new, n_boundaries_new); printf("There are %d nodes, %d elements and %d boundaries in the mesh\n", n_nodes_new, n_elements_new, n_boundaries_new);
int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int)); int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int));
double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double)); double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double));
int *boundaries_new = (int*) calloc(2*n_boundaries_new, sizeof(int)); int *boundaries_new = (int*) calloc(2*n_boundaries_new, sizeof(int));
int *boundary_tags_new = (int*) malloc(n_boundaries_new*sizeof(int)); int *boundary_tags_new = (int*) malloc(n_boundaries_new*sizeof(int));
Tree_cell **map_new = malloc(n_elements_new*sizeof(Tree_cell*)); Tree_cell **map_new = malloc(n_elements_new*sizeof(Tree_cell*));
if (len_node==0 && len_elt==0) {
for (int n=0; n<mesh->n_nodes; n++) { for (int n=0; n<mesh->n_nodes; n++) {
for (int ii=0; ii<3; ii++) { for (int ii=0; ii<3; ii++) {
x_new[n*3+ii] = mesh->x[n*3+ii]; x_new[n*3+ii] = mesh->x[n*3+ii];
}
} }
}else{
for (int n=0; n<mesh->n_nodes; n++) {
int m = n;
if (n>=problem->size_initial_node && n<len_node+problem->size_initial_node) {
m = tab_node[n-problem->size_initial_node];
}else if(n>=len_node+problem->size_initial_node){
m = n-n_nodes_removed;
}
if (m>-1) {
for (int ii=0; ii<3; ii++) {
x_new[m*3+ii] = mesh->x[n*3+ii];
}
}
}
}
deep_tree_travel(mesh, problem->Tree, problem->map, map_new, x_new, elements_new, boundaries_new, boundary_tags_new, n_nodes_new, n_elements_new, n_boundaries_new, problem->size_initial_mesh, len_node, tab_node, problem->size_initial_node, len_elt, tab_elt, problem->size_initial_mesh, n_nodes_removed, n_elements_removed);
/*printf("Print in post_processing files \n");
printing_post_processing_files(mesh, elements_new, x_new, boundaries_new, n_elements_new, n_boundaries_new);
if (mesh->n_nodes>855) {
int tag1 = 855;
printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
}
if (mesh->n_nodes>916) {
int tag1 = 916;
printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
}
if (mesh->n_nodes>1158) {
int tag1 = 1158;
printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
}
if (mesh->n_nodes>1223) {
int tag1 = 1223;
printf("(%f, %f)\n", x_new[tag1*3+0], x_new[tag1*3+1]);
}
FILE *file2 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/element_new.dat", "w");
if(file2 == NULL){
printf("Error at opening file \n");
exit(1);
} }
deep_tree_travel(mesh, problem->Tree, problem->map, map_new, x_new, elements_new, boundaries_new, boundary_tags_new, n_nodes_new, n_elements_new, n_boundaries_new, problem->size_initial_mesh); for (int iel=0; iel<n_elements_new; iel++) {
fprintf(file2, "%4d, %4d, %4d, ", elements_new[iel*3+0], elements_new[iel*3+1], elements_new[iel*3+2]);
}
fclose(file2);*/
free(problem->map); free(problem->map);
problem->map = NULL; problem->map = NULL;
problem->map = map_new; problem->map = map_new;
...@@ -1863,8 +1930,8 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double ...@@ -1863,8 +1930,8 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
free(x_new); free(x_new);
free(boundaries_new); free(boundaries_new);
free(boundary_tags_new); free(boundary_tags_new);
//free(problem->map); free(tab_node);
//free(problem->Tree); free(tab_elt);
} }
} }
......
...@@ -66,7 +66,8 @@ struct FluidProblem { ...@@ -66,7 +66,8 @@ struct FluidProblem {
Tree_cell *Tree; Tree_cell *Tree;
Tree_cell **map; Tree_cell **map;
int *map2Tree; int *map2Tree;
int size_initial_mesh; int size_initial_mesh; // number of initial elements in the mesh
int size_initial_node; // number of initial nodes in the mesh
MeshBoundary *boundaries; MeshBoundary *boundaries;
double *porosity; double *porosity;
double *oldporosity; double *oldporosity;
......
This diff is collapsed.
...@@ -20,7 +20,7 @@ ...@@ -20,7 +20,7 @@
* along with this program (see COPYING and COPYING.LESSER files). If not, * along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>. * see <http://www.gnu.org/licenses/>.
*/ */
#ifndef MESH_H #ifndef MESH_H
#define MESH_H #define MESH_H
...@@ -131,6 +131,7 @@ typedef struct Tree_cell ...@@ -131,6 +131,7 @@ typedef struct Tree_cell
int tag; int tag;
int bool_boundary; int bool_boundary;
int bool_ref; int bool_ref;
int bool_def;
int num_children; int num_children;
int sub_nodes[3]; int sub_nodes[3];
...@@ -152,7 +153,18 @@ void deep_tree_get_mapping(Tree_cell *Tree, Tree_cell **map, int size); ...@@ -152,7 +153,18 @@ void deep_tree_get_mapping(Tree_cell *Tree, Tree_cell **map, int size);
void deep_tree_go_to_leaves(Tree_cell *current, Tree_cell **map); void deep_tree_go_to_leaves(Tree_cell *current, Tree_cell **map);
// REFINEMENT // REFINEMENT
void deep_tree_check_for_refinement(Mesh *mesh, Tree_cell *current, double h_target, double h_min); void deep_tree_check_for_refinement(Mesh *mesh, Tree_cell *current, double h_target, double h_min);
void deep_tree_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell **map, int size, int h_min, int *n_nodes_new, int *n_elements_new, int *n_boundaries_new); void deep_tree_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell **map, int size, int h_min, int *n_nodes_new, int *n_elements_new, int *n_boundaries_new, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt, int *n_nodes_removed, int *n_elements_removed);
// COARSENING
void deep_tree_correspondence(int *tab_node, int len_node, int size_init_node, int *tab_elt, int len_elt, int size_init_elt, int *index_remove, int *elt_remove);
void deep_tree_map_for_remove(Tree_cell *current, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt, int n_nodes_removed, int n_elements_removed, int bool_tag);
void deep_tree_check_for_coarsening(Tree_cell *current, Tree_cell **map, int bool_nghb[3], int bool_tag[3], int n_elements);
void deep_tree_remove_elt(Tree_cell *current, Tree_cell **map, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform4_into_3(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform4_into_2(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform4_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform3_into_2(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform3_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
void deep_tree_transform2_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt);
// CUTTING FUNCTION // CUTTING FUNCTION
void deep_tree_cut_into_4(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd); void deep_tree_cut_into_4(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd);
void deep_tree_cut_into_2(Tree_cell *current, int *elt); void deep_tree_cut_into_2(Tree_cell *current, int *elt);
...@@ -163,8 +175,9 @@ void deep_tree_transform3_into_4(Tree_cell *current, Tree_cell **map, int *index ...@@ -163,8 +175,9 @@ void deep_tree_transform3_into_4(Tree_cell *current, Tree_cell **map, int *index
void deep_tree_transform(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd); void deep_tree_transform(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd);
// CONSTRUCTION OF OUTPUTS DATA STRUCTURE // CONSTRUCTION OF OUTPUTS DATA STRUCTURE
void deep_tree_count_boundaries(Tree_cell *current, int *ibd); void deep_tree_count_boundaries(Tree_cell *current, int *ibd);
void deep_tree_create_data_structure(Mesh *mesh, Tree_cell *current, Tree_cell **map_new, double *x_new, int *elements_new, int *boundaries_new, int *boundary_tags_new, int n_nodes_new, int n_elements_new); void deep_tree_create_data_structure(Mesh *mesh, Tree_cell *current, Tree_cell **map_new, double *x_new, int *elements_new, int *boundaries_new, int *boundary_tags_new, int n_nodes_new, int n_elements_new, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt, int n_nodes_removed, int n_elements_removed);
void deep_tree_travel(Mesh *mesh, Tree_cell *Tree, Tree_cell **map, Tree_cell **map_new, double *x_new, int *elements_new, int *boundaries_new, int *boundary_tags_new, int n_nodes_new, int n_elements_new, int n_boundaries_new, int size); void deep_tree_travel(Mesh *mesh, Tree_cell *Tree, Tree_cell **map, Tree_cell **map_new, double *x_new, int *elements_new, int *boundaries_new, int *boundary_tags_new, int n_nodes_new, int n_elements_new, int n_boundaries_new, int size, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt, int n_nodes_removed, int n_elements_removed);
void deep_tree_travel_to_get_new_correspondence(Tree_cell *current, int len_node, int *tab_node, int size_init_node, int len_elt, int *tab_elt, int size_init_elt, int n_nodes_removed, int n_elements_removed);
// PRINTING TOOLS // PRINTING TOOLS
void printing_post_processing_files(Mesh *mesh, int *elements_new, double *x_new, int *boundaries_new, int n_elements_new, int n_boundaries_new); void printing_post_processing_files(Mesh *mesh, int *elements_new, double *x_new, int *boundaries_new, int n_elements_new, int n_boundaries_new);
void printing_current(Tree_cell *current); void printing_current(Tree_cell *current);
......
L = 0.05; L = 0.1;
H = 0.6; H = 0.6;
y = 0; y = 0;
...@@ -41,7 +41,7 @@ Physical Line("Lateral") = {1,2, 4, 5}; ...@@ -41,7 +41,7 @@ Physical Line("Lateral") = {1,2, 4, 5};
Physical Line("Top") = {6}; Physical Line("Top") = {6};
Physical Surface("Domain") = {1}; Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1}; Physical Point("PtFix") = {1};
Mesh.Algorithm = 8; Mesh.Algorithm = 5;
Merge "lc.pos"; Merge "lc.pos";
Field[1] = PostView; Field[1] = PostView;
......
...@@ -96,37 +96,37 @@ def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) : ...@@ -96,37 +96,37 @@ def genInitialPosition(filename, r, H, ly, lx, rhop, compacity) :
p.write_vtk(filename,0,0) p.write_vtk(filename,0,0)
# Define output directory filename = "/Volumes/Margaux/Memoire_Code/testcases/Boycott/"
outputdir = "output" outputdir = filename + "output"
if not os.path.isdir(outputdir) : if not os.path.isdir(outputdir) :
os.makedirs(outputdir) os.makedirs(outputdir)
# Physical parameters # Physical parameters
g = -9.81 # gravity g = -9.81 # gravity
r = 1e-3 # particles radius r = 1e-3 # particles radius os sand seed
rhop = 2200 # particles density as aluminium rhop = 1600 # particles density as sand
rho = 1000 # fluid density rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity nu = 1e-6 # kinematic viscosity
mu = rho*nu; mu = rho*nu;
compacity = 0.7 compacity = 0.2
# Numerical parameters # Numerical parameters
outf = 5 # number of iterations between output files outf = 5 # number of iterations between output files
dt = 0.1e-2 # time step dt = 0.5e-3 # time step
tEnd = 10 # final time tEnd = 3 # final time
# Geometrical parameters # Geometrical parameters
ly = 0.8 # particles area height ly = 0.8 # particles area height
lx = 0.05 # particles area widht lx = 0.1 # particles area widht
H = 0.6 # domain height H = 0.6 # domain height
# #
# PARTICLE PROBLEM # PARTICLE PROBLEM
# #
# Initialise particles # Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, rhop, compacity) #genInitialPosition(outputdir, r, H, ly, lx, rhop, compacity)
p = scontact.ParticleProblem(2) p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0) p.read_vtk("/Volumes/Margaux/Memoire_Code/testcases/Boycott/output_no_slip_wall_Ref",0)
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0])) print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop) print("RHOP = %g" % rhop)
...@@ -135,12 +135,12 @@ U_0 = 2*9.81*(rhop-rho)*r**2/(9*mu); ...@@ -135,12 +135,12 @@ U_0 = 2*9.81*(rhop-rho)*r**2/(9*mu);
phi = compacity phi = compacity
f_phi = ((1-phi)**2)/((1+phi**(1/3))*np.exp(5/3*phi/(1-phi))); f_phi = ((1-phi)**2)/((1+phi**(1/3))*np.exp(5/3*phi/(1-phi)));
U_0 = 0.065; U_0 = 0.8633;
P_0 = 130.755; P_0 = 255;
Re = rho*U_0*r/mu; Re = rho*U_0*r/mu;
N_0 = len(p.position()); N_0 = len(p.position());
e_target = 1e-7 e_target = 1e-7
print(" ======= The initial number of grains is %.5f [m/s]=======" %N_0) print(" ======= The initial number of grains is %.5f [-]=======" %N_0)
print(" ======= The initial velocity is %.5f [m/s]=======" %U_0) print(" ======= The initial velocity is %.5f [m/s]=======" %U_0)
print(" ======= The initial pressure is %.5f [Pa]=======" %P_0) print(" ======= The initial pressure is %.5f [Pa]=======" %P_0)
print(" ======= The Reynolds number is %.3f [Pa]=======" %Re) print(" ======= The Reynolds number is %.3f [Pa]=======" %Re)
...@@ -157,9 +157,9 @@ delta = 10 ...@@ -157,9 +157,9 @@ delta = 10
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], petsc_solver_type="-pc_type lu -ksp_max_it 50") fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], petsc_solver_type="-pc_type lu -ksp_max_it 50")
# Set the mesh geometry for the fluid computation # Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh") fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom") fluid.set_wall_boundary("Bottom", velocity=[0,0])
fluid.set_wall_boundary("Lateral") fluid.set_wall_boundary("Lateral", velocity=[0,0])
fluid.set_wall_boundary("Top",pressure=0) fluid.set_wall_boundary("Top",pressure=0, velocity=[0,0])
# Set location of the particles in the mesh and compute the porosity in each computation cell # Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces()) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.export_vtk(outputdir,0,0) fluid.export_vtk(outputdir,0,0)
...@@ -195,7 +195,8 @@ while t < tEnd : ...@@ -195,7 +195,8 @@ while t < tEnd :
jj += 1 jj += 1
print("ii=%i, t=%i, jj=%i" %(ii,t,jj)) print("ii=%i, t=%i, jj=%i" %(ii,t,jj))
# load the save past state # load the save past state
fluid.import_vtk("output/fluid_00000.vtu") folder = outputdir+"/fluid_00000.vtu"
fluid.import_vtk(folder)
p.read_vtk(outputdir,0) p.read_vtk(outputdir,0)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces()) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
# create the .msh on the adapt/lc.pos # create the .msh on the adapt/lc.pos
......
L = 0.05; L = 0.1;
H = 0.6; H = 0.6;
y = 0; y = 0;
...@@ -41,6 +41,6 @@ Physical Line("Lateral") = {1,2, 4, 5}; ...@@ -41,6 +41,6 @@ Physical Line("Lateral") = {1,2, 4, 5};
Physical Line("Top") = {6}; Physical Line("Top") = {6};
Physical Surface("Domain") = {1}; Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1}; Physical Point("PtFix") = {1};
Mesh.Algorithm = 8; Mesh.Algorithm = 5;
Mesh.MshFileVersion = 2; Mesh.MshFileVersion = 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