Skip to content

GitLab

Commit 2c0c3679 authored by Margaux Boxho's avatar Margaux Boxho
Browse files

Corrected local adapt

parent 95525b1f
Pipeline #5358 passed with stage
in 28 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 2c0c3679
......@@ -1863,11 +1863,11 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
fluid_problem_set_elements(problem, n_nodes_new, x_new, n_elements_new, elements_new, n_boundaries_new, boundaries_new, boundary_tags_new, mesh->n_boundaries, mesh->boundary_names,1);
/*char filename[128];
char filename[128];
sprintf(filename, "test_%d.msh",global_ind); global_ind++;
FILE *ftmp = fopen(filename,"w");
mesh_write_msh(problem->mesh,ftmp);
fclose(ftmp);*/
fclose(ftmp);
/*printf("HERE \n");*/
problem->map = map_new;
......
fluid/mesh.c
View file @ 2c0c3679
......@@ -618,7 +618,7 @@ void mesh_boundaries_free(MeshBoundary *bs, int n)
// Adding function to handle edges to create a gradient recovery estimator
void initialized_edges_array(Mesh *mesh, Edge *edges_vec, int n_edges){
int tag_edge = 0;
for (int iel=0; iel<mesh->n_elements; iel++) {
......@@ -639,8 +639,8 @@ void initialized_edges_array(Mesh *mesh, Edge *edges_vec, int n_edges){
bool_stop=1; // break the loop
}
ii++;
}
}
if (ii==tag_edge && bool_stop==0) {
edges_vec[ii].tag_Node[0]=el[jj];
edges_vec[ii].tag_Node[1]=el[(jj+1)%3];
......@@ -655,7 +655,7 @@ void initialized_edges_array(Mesh *mesh, Edge *edges_vec, int n_edges){
/*
This function aims to initialize the attribut min_h_target of the structure Mesh.
*/
void mesh_initialized_min_h_target(Mesh *mesh){
void mesh_initialized_min_h_target(Mesh *mesh){
double X1, Y1, X2, Y2, X3, Y3;
double a, b, c;
double J, h_k, vol;
......@@ -1282,7 +1282,9 @@ void deep_tree_count_boundaries(Tree_cell *current, int *ibd){
void deep_tree_create_the_data_structure(Mesh *mesh, Tree_cell *Tree, double *x_new, int n_nodes_new, int *elements_new, int n_elements_new, int *boundaries_new, int n_boundaries_new, int *boundary_tags_new, Tree_cell **map_new, int *map2Tree_new, int size){
int ibd = 0;
printf("Create the data structure size = %d\n", size);
for (int iel=0; iel<size; iel++) {
printf("---------------iel %d---------------\n", iel);
deep_tree_create_x_elt_bd(mesh, &Tree[iel], x_new, elements_new, boundaries_new, boundary_tags_new, &ibd);
}
// CHECK ORIENTATION
......@@ -1304,6 +1306,7 @@ void deep_tree_create_the_data_structure(Mesh *mesh, Tree_cell *Tree, double *x_
for (int iel=0; iel<size; iel++) {
deep_tree_get_map(&Tree[iel], map_new, map2Tree_new, &elt, iel);
}
printf("Printing the data structure \n");
printing_post_processing_files(mesh, elements_new, x_new, boundaries_new, n_elements_new, ibd);
}
......@@ -1313,12 +1316,17 @@ void deep_tree_create_x_elt_bd(Mesh *mesh, Tree_cell *current, double *x_new, in
if (current->parent->num_children==4) {
// ONLY CUTTING IN 4 CREATES NEW NODES
for (int ii=0; ii<3; ii++) {
//printf("HERE FOR CREATING NODES \n");
x_new[(current->nodes[ii])*3+0] = current->XY[ii*2+0];
x_new[(current->nodes[ii])*3+1] = current->XY[ii*2+1];
x_new[(current->nodes[ii])*3+2] = 0.;
if (current->tag==-1) {
printf("ERROR\n");
}
//printf("HERE FOR CREATING ELEMENTS (1) TAG %d\n", current->tag);
elements_new[current->tag*3+ii] = current->nodes[ii];
if (current->neighbours[ii]==-1) {
//printf("tag : %d => NODES (%d, %d)\n", current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
printf("ibd %d, tag : %d => NODES (%d, %d)\n", *ibd, current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
boundaries_new[*ibd*2+0] = current->nodes[ii];
boundaries_new[*ibd*2+1] = current->nodes[(ii+1)%3];
boundary_tags_new[*ibd] = current->phys[ii];
......@@ -1328,9 +1336,13 @@ void deep_tree_create_x_elt_bd(Mesh *mesh, Tree_cell *current, double *x_new, in
}
if (current->parent->num_children==2) {
for (int ii=0; ii<3; ii++) {
if (current->tag==-1) {
printf("ERROR\n");
}
//printf("HERE FOR CREATING ELEMENTS (2) TAG %d\n", current->tag);
elements_new[current->tag*3+ii] = current->nodes[ii];
if (current->neighbours[ii]==-1) {
//printf("tag : %d => NODES (%d, %d)\n", current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
printf("ibd %d, tag : %d => NODES (%d, %d)\n", *ibd, current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
boundaries_new[*ibd*2+0] = current->nodes[ii];
boundaries_new[*ibd*2+1] = current->nodes[(ii+1)%3];
boundary_tags_new[*ibd] = current->phys[ii];
......@@ -1340,9 +1352,13 @@ void deep_tree_create_x_elt_bd(Mesh *mesh, Tree_cell *current, double *x_new, in
}
}else{
for (int ii=0; ii<3; ii++) {
if (current->tag==-1) {
printf("ERROR\n");
}
//printf("HERE FOR CREATING ELEMENTS (3) %d\n", current->tag);
elements_new[current->tag*3+ii] = current->nodes[ii];
if (current->neighbours[ii]==-1) {
//printf("tag : %d => NODES (%d, %d)\n", current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
printf("ibd %d, tag : %d => NODES (%d, %d)\n", *ibd, current->tag, current->nodes[ii], current->nodes[(ii+1)%3]);
boundaries_new[*ibd*2+0] = current->nodes[ii];
boundaries_new[*ibd*2+1] = current->nodes[(ii+1)%3];
boundary_tags_new[*ibd] = current->phys[ii];
......@@ -1591,7 +1607,11 @@ void deep_tree_green_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *current,
int bool_1_node = ((current->sub_nodes[0]!=-1 && current->sub_nodes[1]==-1 && current->sub_nodes[2]==-1) || (current->sub_nodes[0]==-1 && current->sub_nodes[1]!=-1 && current->sub_nodes[2]==-1) || (current->sub_nodes[0]==-1 && current->sub_nodes[1]==-1 && current->sub_nodes[2]!=-1));
int bool_3_nodes = (current->sub_nodes[0]!=-1 && current->sub_nodes[1]!=-1 && current->sub_nodes[2]!=-1);
if (bool_1_node) {
if (current->parent!=NULL) {
if ((current->parent->num_children==2) && (bool_1_node==1)) {
deep_tree_transform_into_4(mesh, Tree, current, index, elt);
}
}else if (bool_1_node) {
// THE CURRENT ELEMENT NEEDS TO BE CUT INTO TWO SUB_ELEMENTS.
printf("CURRENT %d : CUT INTO 2\n", current->tag);
......@@ -1600,8 +1620,7 @@ void deep_tree_green_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *current,
init_children(current);
deep_tree_cut_into_2(current, index, elt);
}
if (bool_3_nodes) {
}else if (bool_3_nodes) {
// THE CURRENT ELEMENT NEEDS TO BE CUT INTO FOUR SUB_ELEMENTS.
printf("CURRENT %d : CUT INTO 4\n", current->tag);
......
testcases/drop-2d/SimpleDrop/drop.py
View file @ 2c0c3679
......@@ -85,7 +85,7 @@ init = True
# Physical parameters
g = -9.81 # gravity
rhop = 2450 # particles density
r = 2*154e-6 #154e-6 # particles radii
r = 2*154e-6 #154e-6 # particles radii
compacity = 0.50 # solid volume fraction in the drop
rho = 1030 # fluid density
nu = 1.17/rho # kinematic viscosity
......@@ -95,9 +95,9 @@ shift = 0.22 # vertical shift of the particle
print("RHOP = %g, r = %g, RHO = %g" % (rhop,r,rho))
# Numerical parameters
outf = 1 # number of iterations between output files
dt = 5e-2 # time step
tEnd = 20 #100 # final time
outf = 5 # number of iterations between output files
dt = 5e-3 # time step
tEnd = 100 # final time
#
# PARTICLE PROBLEM
......@@ -124,7 +124,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# 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.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
......@@ -135,7 +135,7 @@ if init:
#if jj!=0:
# Adaptation of the mesh.
# fluid.adapt_mesh(5e-3,8e-4,50000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
while t < tEnd1 :
# Fluid solver
......@@ -143,7 +143,7 @@ if init:
forces = fluid.compute_node_force(dt1)
# Changes the velocities of the particles without moving the particles
p.velocity()[:] += dt1*forces/p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
# Increase time step if convergence
if niter < 5 and dt1 < 20:
dt1 *= 3
......@@ -156,7 +156,7 @@ t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g
fluid.export_vtk(outputdir,0,0)
tic = time.time()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
#
# COMPUTATION LOOP
......@@ -166,8 +166,8 @@ while t < tEnd :
fluid.implicit_euler(dt,newton_max_it=100)
# Adaptation of the mesh.
if (ii%15==0 and ii != 0):
fluid.adapt_mesh(5e-3,8e-4,50000)
#if (ii%15==0 and ii != 0):
# fluid.adapt_mesh(5e-3,8e-4,50000)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
......@@ -181,7 +181,7 @@ while t < tEnd :
p.iterate(dt/nsub, forces)
t += dt
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
# Output files writting
if ii %outf == 0 :
......
Markdown is supported
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