Local adpat with recursive procedure

fluid/fluid_problem.c

@@ -1712,11 +1712,9 @@ void fluid_problem_set_mesh_and_transfer_solution(FluidProblem *problem, Mesh *n
   free(new_eid);
   free(new_xi);
   free(newx);
-    printf("9 before set_mesh\n");
   fluid_problem_set_mesh(problem,new_mesh);
   free(problem->solution);
   problem->solution = new_solution;
-    printf("10\n");
   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);
@@ -1779,9 +1777,109 @@ void fluid_problem_estimator_evaluation(FluidProblem *problem, double lcmax, dou
     }
 }
 
-void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin){
-    
+
+/*
+ This function aims to create the new mesh according to the estimator. If bool_init_tree is set to 1, it means that it is the first time the mesh needs to be refined thus the 
+ Tree_cell *Tree structure has to be initialized and to be maintained all along the simulation. 
+ 
+ Otherwise, the Tree_cell *Tree structure already exists and a mapping between children tree and mesh elements was established. Using this Tree and the mapping, we are able to
+ refine the current mesh.
+ */
+void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, int bool_init_tree){
     Mesh *mesh = problem->mesh;
+    
+    int n_elements_new = 0;
+    int n_nodes_new = 0;
+    int n_boundaries_new = 0, delta_boundaries_new = 0;
+    
+    if (bool_init_tree) {
+        Tree_cell *Tree = (Tree_cell*) malloc(mesh->n_elements*sizeof(Tree_cell));
+        
+        problem->size_initial_mesh = mesh->n_elements;
+        printf("size = %d\n", problem->size_initial_mesh);
+        
+        create_the_tree(mesh, Tree, lcmin, &n_nodes_new, &n_elements_new, &delta_boundaries_new);
+        
+        for (int kk=0; kk<mesh->n_interfaces; kk++) {
+            if (mesh->interfaces[kk*4+2]==-1) {
+                n_boundaries_new++;
+            }
+        }
+        int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int));
+        double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double));
+        int *boundaries_new = (int*) calloc(2*(n_boundaries_new+delta_boundaries_new), sizeof(int));
+        int *boundary_tags_new = (int*) malloc((n_boundaries_new+delta_boundaries_new)*sizeof(int));
+        
+        // initialize to zero
+        for (int iel=0; iel<n_elements_new; iel++) {
+            for (int kk=0; kk<3; kk++) {
+                elements_new[iel*3+kk]=0;
+            }
+        }
+        read_the_tree(mesh, Tree, x_new, n_nodes_new, elements_new, n_elements_new, boundaries_new, boundary_tags_new, n_boundaries_new);
+        
+        problem->map = (Tree_cell *) malloc(n_elements_new*sizeof(Tree_cell));
+        problem->map2Tree = (int *) malloc(n_elements_new*sizeof(Tree_cell));
+        get_the_mapping(Tree, problem->map, problem->map2Tree, problem->size_initial_mesh);
+        
+        n_boundaries_new += delta_boundaries_new;
+        
+        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);
+        free(elements_new);
+        free(x_new);
+        free(boundaries_new);
+        free(boundary_tags_new);
+        
+        problem->Tree = Tree;
+
+    }else{
+        int n_nodes_new = 0;
+        int n_elements_new = 0;
+        int n_boundaries_new = 0;
+        deep_tree_walk_tree(mesh, problem->Tree, problem->map, problem->map2Tree, problem->size_initial_mesh, lcmin, &n_nodes_new, &n_elements_new, &n_boundaries_new);
+        
+        int *elements_new = (int*) malloc(3*n_elements_new*sizeof(int));
+        double *x_new = (double*) malloc(3*n_nodes_new*sizeof(double));
+        int *boundaries_new = (int*) calloc(2*n_boundaries_new, sizeof(int));
+        int *boundary_tags_new = (int*) malloc(n_boundaries_new*sizeof(int));
+        
+        Tree_cell *map_new = (Tree_cell*) malloc(n_elements_new*sizeof(Tree_cell));
+        int *map2Tree_new = (int*) malloc(n_elements_new*sizeof(Tree_cell));
+        
+        // initialize to zero
+        for (int n=0; n<mesh->n_elements; n++) {
+            for (int ii=0; ii<3; ii++) {
+                x_new[n*3+ii]=mesh->x[n*3+ii];
+            }
+        }
+        for (int n=mesh->n_elements; n<n_elements_new; n++) {
+            for (int ii=0; ii<3; ii++) {
+                x_new[n*3+ii]=0.;
+            }
+        }
+        
+        deep_tree_create_the_data_structure(mesh, problem->Tree, x_new, n_nodes_new, elements_new, n_elements_new, boundaries_new, n_boundaries_new, boundary_tags_new, map_new, map2Tree_new, problem->size_initial_mesh);
+        
+        //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];
+        sprintf(filename, "test_%d.msh",global_ind); global_ind++;
+        FILE *ftmp = fopen(filename,"w");
+        mesh_write_msh(problem->mesh,ftmp);
+        fclose(ftmp);*/
+        
+        /*printf("HERE \n");
+        problem->map = map_new;
+        problem->map2Tree = map2Tree_new;*/
+        
+        /*free(elements_new);
+        free(x_new);
+        free(boundaries_new);
+        free(boundary_tags_new);*/
+        exit(1);
+        
+    }
+    /*Mesh *mesh = problem->mesh;
     int n_elements_new = 0;
     int n_nodes_new = 0;
     int tag_node = mesh->n_nodes;
@@ -1790,7 +1888,7 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
     Tree_cell *Tree = (Tree_cell*) malloc(mesh->n_elements*sizeof(Tree_cell));
     
     create_the_tree(mesh, Tree, lcmin, &n_nodes_new, &n_elements_new, &delta_boundaries_new);
-    
+
     for (int kk=0; kk<mesh->n_interfaces; kk++) {
         if (mesh->interfaces[kk*4+2]==-1) {
             n_boundaries_new++;
@@ -1828,6 +1926,7 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
         }
     }
     // Get the elements_new, x_new and boundaries_new arrays with the desired structure
+    printf("n_boundaries = %d\n", n_boundaries_new);
     read_the_tree(mesh, Tree, x_new, n_nodes_new, elements_new, n_elements_new, boundaries_new, boundary_tags_new, n_boundaries_new);
    
     // Creation of the new mesh and transfert the solution on the new mesh
@@ -1841,11 +1940,11 @@ void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double
     mesh_write_msh(problem->mesh,ftmp);
     fclose(ftmp);
     
-    free(Tree);
+    
     free(elements_new);
     free(x_new);
     free(boundaries_new);
-    free(boundary_tags_new);
+    free(boundary_tags_new);*/
 }
 
 void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {

fluid/fluid_problem.h

@@ -63,6 +63,10 @@ struct FluidProblem {
   double kinetic_energy;
   Mesh *mesh;
   MeshTree *mesh_tree;
+  Tree_cell *Tree;
+  Tree_cell *map;
+  int *map2Tree;
+  int size_initial_mesh;
   MeshBoundary *boundaries;
   double *porosity;
   double *oldporosity;
@@ -106,7 +110,7 @@ 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, double cmax, double cmin);
 void fluid_problem_estimator_evaluation(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int boolean_init);
-void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin);
+void fluid_problem_local_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, int bool_init_tree);
 
 int *fluid_problem_particle_element_id(FluidProblem *problem);
 int fluid_problem_n_particles(FluidProblem *problem);

fluid/mesh.h

@@ -126,17 +126,22 @@ void initialized_edges_array(Mesh *mesh, Edge *edges_vec, int n_edges);
 typedef struct Tree_cell
 {
     struct Tree_cell *children;
+    struct Tree_cell *parent;
+    struct Tree_cell *nghb0;
+    struct Tree_cell *nghb1;
+    struct Tree_cell *nghb2;
+    
+    int tag;
+    int bool_ref;
     int num_children;
-    int *sub_nodes;
-    int *sub_edges;
-    int up;
+    int sub_nodes[3];
     
+    double XY[6];           // to be optimized but I do not have the time
     int nodes[3];       // nodes tag of the current triangle
     int neighbours[3];  // tag of the three neighbors of the current triangle
     int level;          // level of refinement
-    double x[6];        // coordinates of the nodes defining the current triangle
-    int edges[3];       // informtation about the edges : tag of the three edges and boundary information
     int boundary;
+    int phys[3];        // -1 if the edge is not on the boundaries and phys_id if the edge is part of the boundaries.
     
     
 } Tree_cell ;
@@ -144,22 +149,46 @@ typedef struct Tree_cell
 void create_the_tree(Mesh *mesh, Tree_cell *Tree, double h_min, int *n_nodes_new, int *n_elements_new, int *delta_boundaries_new);
 void read_the_tree(Mesh *mesh, Tree_cell *Tree, double *x_new, int n_nodes_new, int *elements_new, int n_elements_new, int *boundaries_new, int *boundary_tags_new, int n_boundaries_new);
 
-int find_neighbours(Mesh *mesh, int current, int n_edges, int neighbours[3], int edges[3]);
+void deep_tree_walk_tree(Mesh *mesh, Tree_cell *Tree, Tree_cell *map, int *map2Tree, int size, double h_min, int *n_nodes_new, int *n_elements_new, int *n_boundaries_new);
+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);
+void deep_tree_get_map(Tree_cell *current, Tree_cell *map_new, int *map2Tree_new, int *elt, int root);
+
+void deep_tree_create_x_elt_bd(Mesh *mesh, Tree_cell *current, double *x_new, int *elements_new, int *boundaries_new, int *boundary_tags_new, int *ibd);
+void deep_tree_init_nghb_step1(Mesh *mesh, Tree_cell *map, Tree_cell *current);
+void deep_tree_init_nghb_step2(Mesh *mesh, Tree_cell *Tree, Tree_cell *map, int *map2Tree, Tree_cell *current);
+void deep_tree_transfert_nghb(Tree_cell *current, Tree_cell *tree_elt, int buffer, int bool_nghb0, int bool_nghb1);
+void deep_tree_transfer(Tree_cell *current, int iel);
+
+void deep_tree_set_bool_ref(Mesh *mesh, Tree_cell *current, double h_min);
+void deep_tree_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *current, double h_min, int *index, int *elt);
+void deep_tree_complete_the_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *current, int *index, int *elt);
+void deep_tree_remove_illegale_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *current, int *index, int *elt, int *ibd, int *num);
+
+void deep_tree_talk_2_neighbours(Tree_cell *map, Tree_cell *current, int iel, int *index);
+void deep_tree_cut_into_4(Mesh *mesh, Tree_cell *Tree, Tree_cell *current, int iel, int *index, int *elt, int bool_2_transfert);
+void deep_tree_cut_into_2(Tree_cell *current, int *index, int *elt);
+void deep_tree_transform_into_4(Mesh *mesh, Tree_cell *Tree, Tree_cell *current, int *index, int *elt);
+void deep_tree_transfer_for_transform(Tree_cell *buf, int ii, Tree_cell *buf_nghb, int *index, int cut_edge, int save_tag1, int save_tag2);
+void deep_tree_green_refinement(Mesh *mesh, Tree_cell *Tree, Tree_cell *map, Tree_cell *current, int *map2Tree, int *index, int *elt);
+
+
+int find_nghb(Mesh *mesh, int current, int n_edges, int neighbours[3], int nodes[6], int phys[3]);
+void order_edges(Tree_cell *current, int nodes[6]);
 
 int need_refinement(double x[6], double h_target, double h_min);
 
+void init_children(Tree_cell *current);
 void talk_2_neighbours(Tree_cell *Tree, int iel, int *index);
 
 int bool_2_nodes_no_children(Tree_cell *Tree, int iel);
 
-void add_boundaries(Mesh *mesh, Tree_cell *Tree, int iel, int *boundaries_new, int *boundary_tags_new, int *ibd);
-int find_tag_in_boundaries(int *boundaries_new, int tag1, int tag2, int n_boundaries);
+void cut_in_2(Tree_cell *Tree, int iel, int *elements_new, int *elt);
+void cut_in_4(Mesh *mesh, Tree_cell *Tree, int iel, double *x_new, int *elements_new, int *elt, int *boundaries_new, int *boundary_tags_new, int *ibd);
 
-void cut_in_2_interior(Mesh *mesh, Tree_cell *Tree, int iel);
-void cut_in_4_interior(Mesh *mesh, Tree_cell *Tree, int iel);
-void cut_in_4_boundary(Mesh *mesh, Tree_cell *Tree, int iel, int *index);
+void get_the_mapping(Tree_cell *Tree, Tree_cell *map, int *map2Tree, int size);
 
 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 print_current(Tree_cell *current, int bool_print);
 #endif
 
 

python/fluid.py

@@ -196,14 +196,14 @@ class FluidProblem :
         self.strong_cb_ref.append(bndcb)
         self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
     
-    def adapt_gen_local_mesh(self, lcmax, lcmin):
+    def adapt_gen_local_mesh(self, lcmax, lcmin, bool_init_tree):
         """ Generate the mesh by local adpative strategy of Bank (REGULAR refinement and GREEN refinement).
             
             Keyword arguments:
             lcmax -- maximum mesh radius
             lcmin -- minimum mesh radius
         """
-        self._lib.fluid_problem_local_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin));
+        self._lib.fluid_problem_local_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_int(bool_init_tree));
     
     def adapt_mesh(self, lcmax, lcmin, n_el, boolean_init, cmax=1, cmin=1) :
         """Create the h_target array linked to the gradient recovery estimator.