/*
 * MigFlow - Copyright (C) <2010-2018>
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
 * 	
 * List of the contributors to the development of MigFlow: see AUTHORS file.
 * Description and complete License: see LICENSE file.
 * 	
 * This program (MigFlow) is free software: 
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program (see COPYING and COPYING.LESSER files).  If not, 
 * see <http://www.gnu.org/licenses/>.
 */

#include "mesh.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "tools.h"
#define N_N (DIMENSION+1)

typedef struct {
  int n[DIMENSION];
  int element;
  int pos; // or tag if element == -1
}HalfEdge;

int hedgecmp(const void *pe0, const void *pe1) {
  int *e0 = ((HalfEdge*)pe0)->n;
  int *e1 = ((HalfEdge*)pe1)->n;
  if (e0[0] != e1[0]) return e0[0]-e1[0];
#if DIMENSION==2
  return e0[1]-e1[1];
#else
  if (e0[1] != e1[1]) return e0[1]-e1[1];
  return e0[2]-e1[2];
#endif
}

static void swapint(int *a, int *b) {
  int c = *a;
  *a = *b;
  *b = c;
}

static void sort_edge_nodes(int *n) {
  if(n[0]>n[1])
    swapint(n,n+1);
#if DIMENSION == 3
  if(n[0]>n[2])
    swapint(n,n+2);
  if(n[1]>n[2])
    swapint(n+1,n+2);
#endif
}

void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
  /*
  for (int i= 0; i < n_bnd_elements; ++i) {
    printf("%i %i : %i\n",bnd_elements[i*2+0],bnd_elements[i*2+1], bnd_tags[i]);
  }
  for (int i = 0; i < m->n_elements; ++i) {
    printf(" el %i : %i %i %i\n",i, m->elements[i*3+0], m->elements[i*3+1], m->elements[i*3+2]);
  }
  */
  int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements;
  HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges);
  HalfEdge *he = halfedges; // current
  int id = 0;
  for (int i=0; i<m->n_elements; ++i) {
    for (int j = 0; j < DIMENSION+1; ++j) {
      for (int k = 0; k < DIMENSION; ++k) {
        he->n[k] = m->elements[i*(DIMENSION+1)+elbnd[j][k]];
      }
      sort_edge_nodes(he->n);
      he->element = i;
      he->pos = j;
      he++;
    }
  }
  for (int i=0; i<n_bnd_elements; ++i) {
    for (int k = 0; k < DIMENSION; ++k){
      he->n[k] = bnd_elements[i*DIMENSION+k];
    }
    sort_edge_nodes(he->n);
    he->element = -1;
    he->pos = bnd_tags[i];
    he ++;
  }
  qsort(halfedges,n_half_edges,sizeof(HalfEdge),hedgecmp);
  int n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)) n_edges++;
  }
  int *edges = malloc(n_edges*4*sizeof(int));
  n_edges = 0;
  for (int i = 0; i < n_half_edges; ++i) {
    if (i == 0 || hedgecmp(halfedges+i-1,halfedges+i)){
      n_edges += 1;
      edges[(n_edges-1)*4+0] = halfedges[i].element;
      edges[(n_edges-1)*4+1] = halfedges[i].pos;
      edges[(n_edges-1)*4+2] = -1;
      edges[(n_edges-1)*4+3] = -1;
    }
    else {
      edges[(n_edges-1)*4+2] = halfedges[i].element;
      edges[(n_edges-1)*4+3] = halfedges[i].pos;
    }
  }
  m->n_interfaces = n_edges;
  m->interfaces = edges;
  free(halfedges);
  for (int i = 0; i < n_edges; ++i) {
    int *e = edges + i*4;
    // ensure boundary tag is on the second part of the closure
    if (e[0] < 0) {
      swapint(&e[0],&e[2]);
      swapint(&e[1],&e[3]);
    }
    // compute the closure id of the second element
    if (e[2] >= 0) {
#if DIMENSION==2
      e[3] += 3;
#else
      int k = 0;
      const int *e0 = &m->elements[e[0]*4];
      const int *e1 = &m->elements[e[2]*4];
      e[3] += 12;
      while (e0[elbnd[e[1]][0]] != e1[elbnd[e[3]][k]]) e[3] += 4;
#endif
    }
    else {
      if (e[3] == -1){
        const int *e0 = &m->elements[e[0]*(DIMENSION+1)];
        printf("boundary edge %i %i with tag -1 !!!\n",e0[elbnd[e[1]][0]],e0[elbnd[e[1]][1]]);
        exit(1);
      }
    }
  }
}


static int check_word_in_file(FILE *f, const char *w) {
  char word[256];
  if(fscanf(f,"%255s", word) != 1){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    return -1;
  }
  if(strcmp(word, w) != 0){
    printf("Invalid mesh file (\"%s\" expected).\n", w);
    return -1;
  }
  return 0;
}

static int read_quoted_string(FILE *f, char *w, size_t maxl) {
  char c;
  size_t pos = 0;
  c = fgetc(f);
  while (c != '\"' && !feof(f)) {
    c = fgetc(f);
  }
  do {
    c = fgetc(f);
    if (c != '\"')
      w[pos++] = c;
  }while (c != '\"' && !feof(f) && pos < maxl);
  if (pos == maxl)
    return printf("string too long in mesh file\n");
  if(feof(f))
    return printf("string not terminated in mesh file\n");
  w[pos] = '\0';
  return 0;
}

static int int_cmp(const void *p0, const void *p1)
{
  int i0 = *(int*)p0;
  int i1 = *(int*)p1;
  if (i0 == i1) return 0;
  if (i0 < i1) return -1;
  return 1;
}

static void sortDi(int i[DIMENSION]) {
#if DIMENSION == 2
   if (i[0]>i[1]) {
     int n = i[0];
     i[0] = i[1];
     i[1] = n;
   }
#else
   int n[3] = {i[0],i[1],i[2]};
   if (n[0]>n[1]) {
      i[1]=n[0];    
      i[0]=n[1];   
   } else {
      i[1]=n[1];  
      i[0]=n[0];  
   } 
   if (i[1]>n[2]) { 
      i[2]=i[1];    
      if(i[0]>n[2]){         
         i[1]=i[0];                
         i[0]=n[2];
      } else i[1]=n[2];      
   } else i[2]=n[2]; 
#endif
}

Mesh *mesh_load_msh(const char *filename)
{
  Mesh *m = malloc(sizeof(Mesh));
  m->interfaces = NULL;
  FILE *f = fopen(filename, "r");
  if (!f){
    printf("Cannot open file \"%s\".\n", filename);
    return NULL;
  }
  check_word_in_file(f,"$MeshFormat");
  float v0;
  int v1, v2;
  if(fscanf(f, "%g %d %d", &v0, &v1, &v2) != 3){
    printf("Cannot parse mesh file version\n");
    return NULL;
  }
  if ((int)10*v0 != 22){
    printf("Mesh file format is %g and not msh 2.2\n", v0);
    return NULL;
  }
  check_word_in_file(f,"$EndMeshFormat");
  check_word_in_file(f,"$PhysicalNames");
  int nphysfile;
  if (fscanf(f, "%d", &nphysfile) != 1){
    printf("Cannot read physical names\n");
    return NULL;
  }
  m->boundary_names = malloc(sizeof(char*)*nphysfile);
  int *physid = malloc(sizeof(int)*nphysfile);
  int nphys = 0;
  for (int i = 0; i < nphysfile; ++i) {
    int physdim, pid;
    if (fscanf(f, "%i %i", &physdim, &pid) != 2){
      printf("Cannot read physical names\n");
      return NULL;
    }
    char pname[256];
    if(read_quoted_string(f, pname, 256) < 0) return NULL;
    if (physdim == DIMENSION-1) {
      physid[nphys] = pid;
      m->boundary_names[nphys] = strdup(pname);
      nphys++;
    }
  }
  m->n_boundaries = nphys;
  check_word_in_file(f,"$EndPhysicalNames");
  check_word_in_file(f,"$Nodes");
  if (fscanf(f, "%d", &m->n_nodes) != 1){
      printf("Cannot read physical nodes.\n");
      return NULL;
  }
  m->x = malloc(sizeof(double)*m->n_nodes*3);

  for (int i = 0; i < m->n_nodes; ++i) {
    int nid;
    if (fscanf(f, "%d %le %le %le", &nid, &m->x[i*3], &m->x[i*3+1], &m->x[i*3+2]) != 4){
      printf("Cannot read nodes.\n");
      return NULL;
    }
    if (nid != i+1){
      printf("Nodes are not sequentialy numbered.\n");
      return NULL;
    }
  }
  check_word_in_file(f,"$EndNodes");
  check_word_in_file(f,"$Elements");
  int n_el;
  if (fscanf(f, "%d", &n_el) != 1){
    printf("Cannot read elements.\n");
    return NULL;
  }
  int *bnd_elements = malloc(sizeof(int)*n_el*DIMENSION);
  int *bnd_tags = malloc(sizeof(int)*n_el);
  size_t n_bnd_elements = 0;
  m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_el);
  m->n_elements = 0;
  for (int i = 0; i < n_el; ++i) {
    int eid, etype, nflags, flag;
    int physfile=0;
    if (fscanf(f, "%d %d %d", &eid, &etype, &nflags) != 3){
      printf("Cannot read elements.\n");
      return NULL;
    }
    for (int j = 0; j < nflags; ++j){
      if (fscanf(f, "%d", &flag) != 1){
        printf("Cannot read elements.\n");
        return NULL;
      }
      if (j == 0)
        physfile = flag;
    }
    int n[DIMENSION+1];
    int phys = -1;
    switch (etype) {
      case(1):
        if (fscanf(f, "%d %d", &n[0], &n[1]) != 2){
          printf("Cannot read elements.\n");
          return NULL;
        }
#if DIMENSION==2
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
          bnd_tags[n_bnd_elements] = phys;
          for (int iD = 0; iD < 2; ++iD) bnd_elements[n_bnd_elements*2+iD] = n[iD]-1;
          sortDi(&bnd_elements[n_bnd_elements*2]);
          n_bnd_elements++;
        }
#endif
        break;
      case(15):
        if (fscanf(f, "%d", &n[0]) != 1){
          printf("Cannot read elements.\n");
          return NULL;
        }
        break;
      case(2):
        if (fscanf(f, "%d %d %d", &n[0], &n[1], &n[2]) != 3){
          printf("Cannot read elements.\n");
          return NULL;
        }
#if DIMENSION == 2
        m->elements[m->n_elements*3+0] = n[0]-1;
        m->elements[m->n_elements*3+1] = n[1]-1;
        m->elements[m->n_elements*3+2] = n[2]-1;
        m->n_elements++;
        break;
#else
        phys = -1;
        for (int in = 0; in < m->n_boundaries; ++in) {
          if (physid[in] == physfile) phys = in;
        }
        if (phys != -1) {
          bnd_tags[n_bnd_elements] = phys;
          for (int iD = 0; iD < 3; ++iD) bnd_elements[n_bnd_elements*3+iD] = n[iD]-1;
          sortDi(&bnd_elements[n_bnd_elements*3]);
          n_bnd_elements++;
        }
        break;
      case(4):
        if (fscanf(f, "%d %d %d %d", &n[0], &n[1], &n[2], &n[3]) != 4){
          printf("Cannot read elements.\n");
          return NULL;
        }
        m->elements[m->n_elements*4+0] = n[0]-1;
        m->elements[m->n_elements*4+1] = n[1]-1;
        m->elements[m->n_elements*4+2] = n[2]-1;
        m->elements[m->n_elements*4+3] = n[3]-1;
        m->n_elements++;
        break;
#endif
      default:
        break;
    }
  }
  m->elements = realloc(m->elements, (DIMENSION+1)*sizeof(int)*m->n_elements);
  check_word_in_file(f,"$EndElements");
  fclose(f);
  free(physid);
  mesh_gen_edges(m,n_bnd_elements,bnd_elements,bnd_tags);
  free(bnd_elements);
  free(bnd_tags);
  return m;
}

int mesh_write_msh(Mesh *mesh, FILE *f)
{
  int nbnd = 0;
  for (int i = 0; i< mesh->n_interfaces; ++i) {
    if (mesh->interfaces[i*4+2] == -1) nbnd++;
  }

  fprintf(f,"$MeshFormat\n2.2 0 0\n$EndMeshFormat\n");
  fprintf(f,"$PhysicalNames\n");
  fprintf(f,"%i\n", mesh->n_boundaries);
  for (int i = 0; i < mesh->n_boundaries; ++i) {
    fprintf(f,"%i %i \"%s\"\n", DIMENSION-1,i+1,mesh->boundary_names[i]);
  }
  fprintf(f,"$EndPhysicalNames\n");
  fprintf(f,"$Nodes\n");
  fprintf(f, "%d\n", mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%i %.16g %.16g %.16g\n", i+1, mesh->x[i*3], mesh->x[i*3+1], mesh->x[i*3+2]);
  }
  fprintf(f, "$EndNodes\n");
  fprintf(f, "$Elements\n");
  fprintf(f, "%i\n", mesh->n_elements+nbnd);
  for (int i= 0; i < mesh->n_elements; ++i) {
#if DIMENSION==2
    int *tri = mesh->elements+i*3;
    fprintf(f, "%i 2 2 1 1 %i %i %i\n", i+1, tri[0]+1, tri[1]+1, tri[2]+1);
#else
    int *tet = mesh->elements+i*4;
    fprintf(f, "%i 4 2 1 1 %i %i %i %i\n", i+1, tet[0]+1, tet[1]+1, tet[2]+1, tet[3]+1);
#endif
  }
  nbnd = 0;
  for (int i = 0; i< mesh->n_interfaces; ++i) {
    if (mesh->interfaces[i*4+2] != -1) continue;
    int *cl = elbnd[mesh->interfaces[i*4+1]];
    int iel = mesh->interfaces[i*4+0];
    int phys = mesh->interfaces[i*4+3];
    int imsh = nbnd+mesh->n_elements;
#if DIMENSION==2
    int *tri = mesh->elements+iel*3;
    fprintf(f, "%i 1 2 %i %i %i %i\n", imsh+1, phys+1,phys+1,tri[cl[0]]+1, tri[cl[1]]+1);
#else
    int *tet = mesh->elements+iel*4;
    fprintf(f, "%i 2 2 %i %i %i %i %i\n", imsh+1, phys+1,phys+1,tet[cl[0]]+1, tet[cl[1]]+1, tet[cl[2]]+1);
#endif
    nbnd++;
  }
  fprintf(f, "$EndElements\n");
  return 0;
}

void mesh_free(Mesh *m)
{
  free(m->elements);
  free(m->x);
  for (int i = 0; i < m->n_boundaries; ++i) {
    free(m->boundary_names[i]);
  }
  free(m->boundary_names);
  free(m->interfaces);
  free(m);
}

int mesh_write_msh_scalar(const Mesh *mesh, FILE *f, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id)
{
  fprintf(f,"$NodeData\n");
  fprintf(f, "1\n\"%s\"\n",field_name);
  fprintf(f, "1\n%.16g\n",t);
  fprintf(f, "3\n%i\n%i\n%i\n", iter, 1, mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
    fprintf(f, "%i %.16g\n", i+1, solution[i*n_fields+field_id]);
  }
  fprintf(f, "$EndNodeData\n");
  return 0;
}

int mesh_write_msh_vector(const Mesh *mesh, FILE *f, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id[DIMENSION])
{
  fprintf(f,"$NodeData\n");
  fprintf(f, "1\n\"%s\"\n",field_name);
  fprintf(f, "1\n%.16g\n",t);
  fprintf(f, "3\n%i\n%i\n%i\n", iter, 3, mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i) {
#if DIMENSION==2
    fprintf(f, "%i %.16g %.16g 0.\n", i+1, solution[i*n_fields+field_id[0]], solution[i*n_fields+field_id[1]]);
#else
    fprintf(f, "%i %.16g %.16g %.16g\n", i+1, solution[i*n_fields+field_id[0]], solution[i*n_fields+field_id[1]], solution[i*n_fields+field_id[2]]);
#endif
  }
  fprintf(f, "$EndNodeData\n");
  return 0;
}

int mesh_write_pos_scalar(const Mesh *mesh, const char *dir_name, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id)
{
  char file_name[1024];
  sprintf(file_name,"%s/%s_%05i.pos",dir_name, field_name, iter);
  FILE *f = fopen(file_name, "w");
  if (!f)
    return printf("Cannot open file \"%s\" for writing.\n", file_name);
  fprintf(f,"View\"%s\"{\nTIME {%g};\n", field_name, t);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION==2
    int *tri = mesh->elements+iel*3;
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    double V[3] = {solution[tri[0]*n_fields+field_id], solution[tri[1]*n_fields+field_id], solution[tri[2]*n_fields+field_id]};
    fprintf(f, "ST(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
        V[0],V[1],V[2]);
#else
    int *tet = mesh->elements+iel*4;
    const double x[4][3] = {
      {mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
      {mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
      {mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
      {mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]},
    };
    double V[4] = {solution[tet[0]*n_fields+field_id], solution[tet[1]*n_fields+field_id], solution[tet[2]*n_fields+field_id], solution[tet[3]*n_fields+field_id]};
    fprintf(f, "ST(%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g){%.8g, %.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],x[2][1],x[2][2],x[3][0],x[3][1],x[3][2],
        V[0],V[1],V[2],V[3]);
#endif
  }
  fprintf(f, "};\n");
  fclose(f);
  return 0;
}

int mesh_write_pos_vector(const Mesh *mesh, const char *dir_name, const char *field_name, double t, int iter, const double *solution, int n_fields, int field_id[DIMENSION])
{
  char file_name[1024];
  sprintf(file_name,"%s/%s_%05i.pos",dir_name, field_name, iter);
  FILE *f = fopen(file_name, "w");
  if (!f)
    return printf("Cannot open file \"%s\" for writing.\n", file_name);
  fprintf(f,"View\"%s\"{\nTIME {%g};\n", field_name, t);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
#if DIMENSION ==2
    int *tri = mesh->elements+iel*3;
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    double V0[3] = {solution[tri[0]*n_fields+field_id[0]], solution[tri[1]*n_fields+field_id[0]], solution[tri[2]*n_fields+field_id[0]]};
    double V1[3] = {solution[tri[0]*n_fields+field_id[1]], solution[tri[1]*n_fields+field_id[1]], solution[tri[2]*n_fields+field_id[1]]};
    fprintf(f, "VT(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0};\n",
        x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
        V0[0],V1[0], V0[1], V1[1], V0[2], V1[2]);
#else
    int *tet = mesh->elements+iel*4;
    const double x[4][3] = {
      {mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
      {mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
      {mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
      {mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]},
    };
    double V0[4] = {solution[tet[0]*n_fields+field_id[0]], solution[tet[1]*n_fields+field_id[0]], solution[tet[2]*n_fields+field_id[0]], solution[tet[3]*n_fields+field_id[0]]};
    double V1[4] = {solution[tet[0]*n_fields+field_id[1]], solution[tet[1]*n_fields+field_id[1]], solution[tet[2]*n_fields+field_id[1]], solution[tet[3]*n_fields+field_id[1]]};
    double V2[4] = {solution[tet[0]*n_fields+field_id[2]], solution[tet[1]*n_fields+field_id[2]], solution[tet[2]*n_fields+field_id[2]], solution[tet[3]*n_fields+field_id[2]]};
    fprintf(f, "VT(%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g){%.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g, %.8g};\n",
        x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],x[2][1],x[2][2],x[3][0],x[3][1],x[3][2],
        V0[0],V1[0],V2[0], V0[1], V1[1],V2[1], V0[2], V1[2],V2[2], V0[3], V1[3],V2[3]);
#endif
  }
  fprintf(f, "};\n");
  fclose(f);
  return 0;
}

static int intcmp(const void *p0, const void *p1) {
  return *(int*)p0 - *(int*)p1;
}

MeshBoundary *mesh_boundaries_create(const Mesh *m) {
  MeshBoundary *b = malloc(sizeof(MeshBoundary)*m->n_boundaries);
  for (int i = 0; i < m->n_boundaries; ++i) {
    b[i].n_nodes = 0;
    b[i].n_interfaces = 0;
  }
  for (int i = 0; i < m->n_interfaces; ++i) {
    int *interface = m->interfaces + i*4;
    if (interface[2] == -1) {
      b[interface[3]].n_interfaces++;
    }
  }
  for (int i = 0; i < m->n_boundaries; ++i) {
    b[i].interfaces = malloc(b[i].n_interfaces*sizeof(int));
    b[i].nodes = malloc(b[i].n_interfaces*sizeof(int)*DIMENSION);
    b[i].n_interfaces = 0;
  }
  for (int i = 0; i < m->n_interfaces; ++i) {
    int *interface = m->interfaces + i*4;
    if (interface[2] == -1) {
      int tag = interface[3];
      int *el = m->elements + interface[0]*(DIMENSION+1);
      int *cl = elbnd[interface[1]];
      b[tag].interfaces[b[tag].n_interfaces++] = i;
      for (int i = 0; i < DIMENSION; ++i) {
        b[tag].nodes[b[tag].n_nodes++] = el[cl[i]];
      }
    }
  }
  for (int i = 0; i < m->n_boundaries; ++i) {
    qsort(b[i].nodes,b[i].n_nodes,sizeof(int),intcmp);
    int end = b[i].n_nodes;
    b[i].n_nodes = 0;
    int p = -1;
    for (int j = 0; j < end; ++j) {
      if (b[i].nodes[j] == p){
        continue;
      }
      b[i].nodes[b[i].n_nodes++] = b[i].nodes[j];
      p = b[i].nodes[j];
    }
    b[i].nodes = realloc(b[i].nodes, b[i].n_nodes*sizeof(int));
  }
  return b;
}

void mesh_boundaries_free(MeshBoundary *bs, int n)
{
  for (int i =0; i < n; ++i) {
    free(bs[i].nodes);
    free(bs[i].interfaces);
  }
  free(bs);
}

// 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++) {
        const unsigned int *el = &mesh->elements[iel*N_N];
        for (int jj=0; jj<3; jj++){
            int tag_node1 = el[jj];
            int tag_node2 = el[(jj+1)%3];
            
            int ii = 0; int bool_stop = 0;
            
            while (ii<tag_edge && bool_stop==0) {
                
                int tag_node1_test = edges_vec[ii].tag_Node[0];
                int tag_node2_test = edges_vec[ii].tag_Node[1];
                
                if ((tag_node1_test==tag_node1 && tag_node2_test==tag_node2) || (tag_node1_test==tag_node2 && tag_node2_test==tag_node1)){
                    edges_vec[ii].tag_Tri[1]=iel+1; // the edge exits and we add its reference to the edge triangle
                    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];
                edges_vec[ii].tag_Tri[0]=iel+1; // such that 0 is a triangle that does not exist
                edges_vec[ii].tag=tag_edge+1;
                tag_edge++;
            }
        }
    }
}

/*
 This function aims to initialize the attribut min_h_target of the structure 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;

    mesh->min_h_target = (double*) calloc(mesh->n_elements, sizeof(double));
    for (int iel=0; iel<mesh->n_elements; iel++) {
        const unsigned int *el = &mesh->elements[iel*N_N];
        X1 = mesh->x[el[0]*3+0]; Y1 = mesh->x[el[0]*3+1];
        X2 = mesh->x[el[1]*3+0]; Y2 = mesh->x[el[1]*3+1];
        X3 = mesh->x[el[2]*3+0]; Y3 = mesh->x[el[2]*3+1];
        J = (X2-X1)*(Y3-Y1) - (X3-X1)*(Y2-Y1);
        vol = J/2.0;
        a = sqrt((X1-X2)*(X1-X2)+(Y1-Y2)*(Y1-Y2));
        b = sqrt((X2-X3)*(X2-X3)+(Y2-Y3)*(Y2-Y3));
        c = sqrt((X3-X1)*(X3-X1)+(Y3-Y1)*(Y3-Y1));
        h_k = 0.5*a*b*c/vol;
        //double rrr = (rand() % 200)/100.0;
        //mesh->min_h_target[iel] = h_k*rrr;
        mesh->min_h_target[iel] = h_k;
        //mesh->min_h_target[iel] = 0.5*h_k*(iel==0) + 0.5*h_k*(iel==1) + 0.5*h_k*(iel==2) + 0.5*h_k*(iel==3);
    }
}

// ****************************************************************************************************************************************************************
// LOCAL MESH ADAPTATION USING RECURSION ==> KEEP AND MAINTAIN THE TREE STRUCTURE ITERATION AFTER ITERATION
// ****************************************************************************************************************************************************************

// INITIALIZATION OF THE TREE
void deep_tree_initialized(Mesh *mesh, Tree_cell *Tree, int size){
    for (int iel=0; iel<size; iel++) {
        Tree_cell *current = &Tree[iel];
        current->children = NULL;
        current->parent = NULL;
        current->num_children = 0;
        current->tag = iel;
        for (int ii=0; ii<3; ii++) {
            current->sub_nodes[ii]=-1;
            current->phys[ii]=-1;
            current->neighbours[ii]=-1;
        }
        
        const int *el = &mesh->elements[iel*N_N];
        for (int ii=0; ii<3; ii++) {
            current->nodes[ii] = el[ii];
        }
    }
}

void deep_tree_get_mapping(Tree_cell *Tree, Tree_cell **map, int size){
    for (int iel=0; iel<size; iel++) {
        deep_tree_go_to_leaves(&Tree[iel], map);
    }
}

void deep_tree_go_to_leaves(Tree_cell *current, Tree_cell **map){
    if (current->children==NULL) {
        map[current->tag] = current;
    }else{
        for (int ii=0; ii<current->num_children; ii++) {
            deep_tree_go_to_leaves(&current->children[ii], map);
        }
    }
}

void initialize_children(Tree_cell *current, int excess){
    for (int ii=0; ii<current->num_children+excess; ii++) {
        current->children[ii].parent=current;
        current->children[ii].children=NULL;
        current->children[ii].num_children = 0;
        current->children[ii].bool_ref = 0;
        current->children[ii].bool_def = 0;
        for (int jj=0; jj<3; jj++) {
            current->children[ii].neighbours[jj]=-1;
            current->children[ii].phys[jj]=-1;
            current->children[ii].sub_nodes[jj]=-1;
            current->children[ii].bdies[jj]=-1;
        }
    }
}

void deep_tree_check_for_refinement(Mesh *mesh, Tree_cell *current, double h_target, double h_min){
    double x[6] = {mesh->x[current->nodes[0]*3+0], mesh->x[current->nodes[0]*3+1], mesh->x[current->nodes[1]*3+0], mesh->x[current->nodes[1]*3+1], mesh->x[current->nodes[2]*3+0], mesh->x[current->nodes[2]*3+1]};
    double alpha = 1.5;
    double J = (x[2]-x[0])*(x[5]-x[1]) - (x[4]-x[0])*(x[3]-x[1]);
    double vol = J/2.0;
    double a = sqrt((x[0]-x[2])*(x[0]-x[2])+(x[1]-x[3])*(x[1]-x[3]));
    double b = sqrt((x[2]-x[4])*(x[2]-x[4])+(x[3]-x[5])*(x[3]-x[5]));
    double c = sqrt((x[4]-x[0])*(x[4]-x[0])+(x[5]-x[1])*(x[5]-x[1]));
    double h_k = 0.5*a*b*c/vol;
    if (h_k/1.4>h_target) {
        current->bool_ref = 1;
        current->bool_def = 0;
    }else if(h_k*alpha<h_target){
        current->bool_ref = 0;
        current->bool_def = 0; //1;
    }else{
        current->bool_def = 0;
        current->bool_ref = 0;
    }
}

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){
    int index = mesh->n_nodes;
    int elt = mesh->n_elements;
    int index_remove = 0;
    int elt_remove = 0;
    
    int ibd = 0;
    // NEED TO CARRY ABOUT THE NEIGHBOURS
    for (int kk=0; kk<mesh->n_interfaces; kk++) {
        int *interface = &mesh->interfaces[kk*4];
        int el0 = interface[0];
        int cl0 = interface[1];
        int el1 = interface[2];
        int cl1 = interface[3];
        
        if (el1==-1) {
            map[el0]->phys[cl0%3] = cl1;
            map[el0]->neighbours[cl0%3] = -1;
            map[el0]->bdies[cl0%3] = ibd;
            ibd++;
        }
        else {
            map[el0]->neighbours[cl0%3] = el1;
            map[el0]->bdies[cl0%3] = -1;
            map[el1]->neighbours[cl1%3] = el0;
            map[el1]->bdies[cl1%3] = -1;
        }
    }
    // ASK FOR REFINEMENT TO THE ESTIMATOR
    for (int iel=0; iel<mesh->n_elements; iel++) {
        deep_tree_check_for_refinement(mesh, map[iel], mesh->min_h_target[iel], h_min);
        Tree_cell *current = map[iel];
        if (current->children!=NULL && current->num_children==0) {
            free(current->children);
            current->children=NULL;
        }
        //printf("TAG %4d : (%4d, %4d, %4d)\n", current->tag, current->nodes[0], current->nodes[1], current->nodes[2]);
    }
    
    // CUT INTO 4 ELEMENTS
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        int bool_boundary = 0;
        if (current->parent!=NULL) {
            if (current->parent->num_children==4 && current->bool_ref==1 && bool_boundary==0) {
                //printf("====>CUT INTO 4\n");
                current->children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
                current->num_children=4;
                initialize_children(current, 0);
                deep_tree_cut_into_4(current, map, &index, &elt, &ibd);
            }
        }else if(current->bool_ref==1 && bool_boundary==0){
            //printf("====>CUT INTO 4\n");
            current->children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
            current->num_children=4;
            initialize_children(current, 0);
            deep_tree_cut_into_4(current, map, &index, &elt, &ibd);
        }
    }
    
    // TRANSFORM ELEMENT CUT INTO 2 OR 3 INTO ELEMENT CUT INTO 4
    printf("=============================> STARTING TRANSFORM ELEMENT <=============================\n");
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        if (current->children==NULL) {
            // WE ONLY LOOK FOR ELEMENT THAT HAVE A PARENT AND 1 OR 2 BROTHERS
            if (current->parent!=NULL) {
                // DEFINE BOOLEAN
                int bool_one_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_two_nodes = (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);
                
                if (current->parent->num_children==2 ) {
                    // TRANSFORM A 2-CUT INTO A 4-CUT
                    if (current->bool_ref || bool_one_node || bool_two_nodes) {
                        //printf("====>Transform 2 into 4\n");
                        deep_tree_transform2_into_4(current, map, &index, &elt, &ibd);
                    }
                }else if (current->parent->num_children==3){
                    // TRANSFORM A 3-CUT INTO A 4-CUT
                    if (current->bool_ref || bool_one_node || bool_two_nodes) {
                        //printf("====>Transform 3 into 4\n");
                        deep_tree_transform3_into_4(current, map, &index, &elt, &ibd);
                    }
                }
            }
        }
    }
    
    printf("=============================> STARTING COARSENING THE MESH <=============================\n");
    FILE *file1 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/Mesh_Remove.dat", "w");
    if(file1 == NULL){
        printf("Error at opening file \n");
        exit(1);
    }
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        if (current->bool_def==1 && current->tag!=-1) {
            if (current->parent!=NULL) {
                for (int ii=0; ii<3; ii++) {
                    int tag_node_remove = current->parent->nodes[ii];
                    fprintf(file1, "%f, %f, ", mesh->x[tag_node_remove*3+0], mesh->x[tag_node_remove*3+1]);
                }
                fprintf(file1, "\n");
            }
            deep_tree_remove_elt(current, map, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt);
        }
    }
    fclose(file1);
    
    FILE *file2 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/tab_node_before.dat", "w");
    if(file2 == NULL){
        printf("Error at opening file \n");
        exit(1);
    }
    for (int ii=0; ii<len_node; ii++) {
        fprintf(file2, "%4d, ", tab_node[ii]);
    }
    fclose(file2);
    
    printf("=============================> STARTING SMOOTHING THE REFINEMENT (0)<=============================\n");
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        if (current->parent!=NULL && current->tag!=-1 && current->children==NULL) {
            Tree_cell *parent = current->parent;
            int count_sub_nodes = 0;
            for (int ii=0; ii<3; ii++) {
                if (parent->sub_nodes[ii]!=-1) {
                    count_sub_nodes++;
                }
            }
            if (count_sub_nodes != parent->num_children-1 && len_node!=0) {
                if (count_sub_nodes==2 && parent->num_children==4) {
                    //printing_current(parent);
                    //printf("Transform 4->3\n");
                    deep_tree_transform4_into_3(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else if (count_sub_nodes==1 && parent->num_children==4){
                    //printing_current(parent);
                    //printf("Transform 4->2\n");
                    deep_tree_transform4_into_2(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else if (count_sub_nodes==0 && parent->num_children==4){
                    //printing_current(parent);
                    //printf("Transform 4->1\n");
                    deep_tree_transform4_into_1(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else if (count_sub_nodes==1 && parent->num_children==3){
                    //printf("Transform 3->2\n");
                    deep_tree_transform3_into_2(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else if (count_sub_nodes==0 && parent->num_children==3){
                    //printf("Transform 3->1\n");
                    deep_tree_transform3_into_1(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else if (count_sub_nodes==0 && parent->num_children==2){
                    //printf("Transform 2->1\n");
                    deep_tree_transform2_into_1(parent, len_elt, tab_elt, size_init_elt);
                    //printing_current(parent);
                }else{
                    printf("ERROR IN SMOOTHING PROCEDURE (0)\n");
                    printing_current(parent);
                    exit(1);
                }
            }
        }
    }
    
    printf("=============================> STARTING SMOOTHING THE REFINEMENT (1)<=============================\n");
    // COMPLETE THE REFINEMENT BY CUTTING INTO 2 OR 3 DEENDING ON THE NUMBER OF SUB_NODES CONTAINED ON THE ELEMENT EDGES
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        if (current->children==NULL && current->parent!=NULL) {
            if (current->parent->num_children==2) {
                int tag = 0;
                for(int ii=0; ii<3; ii++){
                    if(current->parent->sub_nodes[ii]==-1 && current->parent->sub_nodes[(ii+1)%3]==-1){
                        tag = (ii+2)%3;
                    }
                }
                if (tag == 0) {
                    current->parent->sub_nodes[1] = current->parent->children[1].sub_nodes[1];
                    current->parent->sub_nodes[2] = current->parent->children[0].sub_nodes[2];
                }else if(tag == 1){
                    current->parent->sub_nodes[2] = current->parent->children[1].sub_nodes[1];
                    current->parent->sub_nodes[0] = current->parent->children[0].sub_nodes[2];
                }else{
                    current->parent->sub_nodes[0] = current->parent->children[1].sub_nodes[1];
                    current->parent->sub_nodes[1] = current->parent->children[0].sub_nodes[2];
                }
                int bool_two_nodes = (current->parent->sub_nodes[0]!=-1 && current->parent->sub_nodes[1]!=-1 && current->parent->sub_nodes[2]==-1)||(current->parent->sub_nodes[0]==-1 && current->parent->sub_nodes[1]!=-1 && current->parent->sub_nodes[2]!=-1)||(current->parent->sub_nodes[0]!=-1 && current->parent->sub_nodes[1]==-1 && current->parent->sub_nodes[2]!=-1);
                int bool_three_nodes = (current->parent->sub_nodes[0]!=-1 && current->parent->sub_nodes[1]!=-1 && current->parent->sub_nodes[2]!=-1);
                int save_tag[2] = {current->parent->children[0].tag, current->parent->children[1].tag};
                
                if (bool_two_nodes) {
                    current->parent->num_children = 3;
                    initialize_children(current->parent, 1);
                    deep_tree_cut_into_3(current->parent, &elt);
                    elt-=2;
                    for (int ii=0; ii<2; ii++) {
                        current->parent->children[ii].tag = save_tag[ii];
                    }
                    current->parent->children[2].tag = elt;
                    elt+=1;
                    
                }else if(bool_three_nodes){
                    current->parent->num_children = 4;
                    initialize_children(current->parent, 0);
                    deep_tree_cut_into_4(current->parent, map, &index, &elt, &ibd);
                    elt-=3;
                    for (int ii=0; ii<2; ii++) {
                        current->parent->children[ii].tag = save_tag[ii];
                    }
                    current->parent->children[2].tag = elt; elt+=1;
                    current->parent->children[3].tag = elt; elt+=1;
                    
                }
            }else if (current->parent->num_children==3){
                int tag = 0;
                for(int ii=0; ii<3; ii++){
                    if(current->parent->sub_nodes[ii]==-1){
                        tag = ii;
                    }
                }
                current->parent->sub_nodes[tag] = current->parent->children[0].sub_nodes[2];
                
                int bool_three_nodes = (current->parent->sub_nodes[0]!=-1 && current->parent->sub_nodes[1]!=-1 && current->parent->sub_nodes[2]!=-1);
                int save_tag[3] = {current->parent->children[0].tag, current->parent->children[1].tag, current->parent->children[2].tag};
                if (bool_three_nodes) {
                    current->parent->num_children = 4;
                    initialize_children(current->parent, 0);
                    deep_tree_cut_into_4(current->parent, map, &index, &elt, &ibd);
                    elt-=3;
                    for (int ii=0; ii<3; ii++) {
                        current->parent->children[ii].tag = save_tag[ii];
                    }
                    current->parent->children[3].tag = elt; elt+=1;
                }
            }
        }
    }
    
    printf("=============================> STARTING SMOOTHING THE REFINEMENT (2)<=============================\n");
    // COMPLETE THE REFINEMENT BY CUTTING INTO 2 OR 3 DEENDING ON THE NUMBER OF SUB_NODES CONTAINED ON THE ELEMENT EDGES
    for (int iel=0; iel<mesh->n_elements; iel++) {
        Tree_cell *current = map[iel];
        int num_child = 0;
        if (current->parent!=NULL) {
            num_child = current->parent->num_children;
        }
        if ((current->children==NULL && current->parent==NULL) || (current->children==NULL && num_child==4)) {
            int bool_one_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_two_nodes = (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_three_nodes = (current->sub_nodes[0]!=-1 && current->sub_nodes[1]!=-1 && current->sub_nodes[2]!=-1);
            
            if (bool_one_node) {
                //printf("====>CUT INTO 2\n");
                current->children = (Tree_cell*) malloc(4*sizeof(Tree_cell));
                current->num_children = 2;
                initialize_children(current, 2);
                deep_tree_cut_into_2(current, &elt);
            }else if(bool_two_nodes){
                //printf("====>CUT INTO 3\n");
                current->children = (Tree_cell*) malloc(4*sizeof(Tree_cell));
                current->num_children = 3;
                initialize_children(current, 1);
                deep_tree_cut_into_3(current, &elt);
            }else if(bool_three_nodes){
                //printf("====>CUT INTO 4\n");
                current->children = (Tree_cell*) malloc(4*sizeof(Tree_cell));
                current->num_children = 4;
                initialize_children(current, 0);
                deep_tree_cut_into_4(current, map, &index, &elt, &ibd);
            }
        }
    }
    
    
    /*FILE *file3 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/tab_node_after.dat", "w");
    if(file3 == NULL){
        printf("Error at opening file \n");
        exit(1);
    }
    for (int ii=0; ii<len_node; ii++) {
        fprintf(file3, "%4d, ", tab_node[ii]);
    }
    fclose(file3);*/
    
    deep_tree_correspondence(tab_node, len_node, size_init_node, tab_elt, len_elt, size_init_elt, &index_remove, &elt_remove);
    
    /*FILE *file4 = fopen("/Users/margauxboxho/migflow/testcases/my_local_drop/tab_node_after_after.dat", "w");
    if(file4 == NULL){
        printf("Error at opening file \n");
        exit(1);
    }
    for (int ii=0; ii<len_node; ii++) {
        fprintf(file4, "%4d, ", tab_node[ii]);
    }
    fclose(file4);*/
    
    *n_nodes_new = index-index_remove;
    *n_elements_new = elt-elt_remove;
    *n_boundaries_new = ibd; // we do not touch to the boundaries
    
    printf("Number of nodes to be removed %4d\n", index_remove);
    printf("Number of elements to be removed %4d\n", elt_remove);

    *n_nodes_removed = index_remove;
    *n_elements_removed = elt_remove;
}

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){
    printf("====> CREATE THE DATA STRUCTURE <====\n");
    int ibd = 0;
    for (int iel=0; iel<size; iel++) {
        Tree_cell *current = &Tree[iel];
        deep_tree_create_data_structure(mesh, current, map_new, x_new, elements_new, boundaries_new, boundary_tags_new, n_nodes_new, n_elements_new, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, n_elements_removed);
    }
     for (int iel=0; iel<mesh->n_elements; iel++) {
         map[iel]->bool_boundary=1;
    }
    
    printf("====> CREATE THE BOUNDARY STRUCTURE <====\n");
    int count = 0;
    for (int iel=0; iel<n_elements_new; iel++) {
        Tree_cell *current = map_new[iel];
        for (int ii=0; ii<3; ii++) {
            int ibd = current->bdies[ii];
            if (ibd>-1) {
                count++;
                boundaries_new[2*ibd+0] = current->nodes[ii];
                boundaries_new[2*ibd+1] = current->nodes[(ii+1)%3];
                boundary_tags_new[ibd] = current->phys[ii];
            }
        }
    }
    
    printf("n_boundaries_new = %d \n", n_boundaries_new);
    
    if (n_nodes_removed!=0 && n_elements_removed!=0) {
        for (int iel=0; iel<size; iel++) {
            Tree_cell *current = &Tree[iel];
            deep_tree_travel_to_get_new_correspondence(current, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, n_elements_removed);
        }
    }
}

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){
    
    if (current->children==NULL || current->num_children==0) { // the second statement will be used when elements are removed
        
        int save_old_nodes[3];
        if (current->parent!=NULL) {
            if (current->parent->num_children==4){
                for (int ii=0; ii<3; ii++) {
                    save_old_nodes[ii] = current->parent->nodes[ii];
                }
            }
        }
        
        // Get the new tag for the nodes, sub_nodes, neighbours and the current element
        if (len_node!=0 && len_elt!=0 && n_nodes_removed!=0 && n_elements_removed!=0) {
            deep_tree_map_for_remove(current, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, n_elements_removed, 1);
        }
        
        if (current->tag>n_elements_new) {
            printf("ERROR IN THE NODE TAG\n");
            printing_current(current);
            exit(1);
        }
        
        if (current->tag == -1){
            printf("ERROR THE CURRENT ELEMENT HAS A TAG -1 !!!!\n");
            printing_current(current);
            printing_current(current->parent);
            exit(1);
        }
        
        current->bool_boundary = 1;
        map_new[current->tag] = current;
        
        elements_new[current->tag*3+0] = current->nodes[0];
        elements_new[current->tag*3+1] = current->nodes[1];
        elements_new[current->tag*3+2] = current->nodes[2];
        
        if (current->parent!=NULL) {
            int count_sub_nodes = 0;
            for (int ii=0; ii<3; ii++) {
                if (current->parent->sub_nodes[ii]!=-1) {
                    count_sub_nodes ++;
                }
            }
            if (current->parent->num_children==4 && count_sub_nodes==3) {
                if (save_old_nodes[0]>mesh->n_nodes || save_old_nodes[1]>mesh->n_nodes || save_old_nodes[2]>mesh->n_nodes) {
                    printf("(%d, %d, %d) \n", save_old_nodes[0], save_old_nodes[1], save_old_nodes[2]);
                    printf("ERROR IN THE PARENT NODES TAG\n");
                    printing_current(current->parent);
                    exit(1);
                }
                // Extract node using their old tag
                double X1 = mesh->x[save_old_nodes[0]*3+0]; double Y1 = mesh->x[save_old_nodes[0]*3+1];
                double X2 = mesh->x[save_old_nodes[1]*3+0]; double Y2 = mesh->x[save_old_nodes[1]*3+1];
                double X3 = mesh->x[save_old_nodes[2]*3+0]; double Y3 = mesh->x[save_old_nodes[2]*3+1];
                
                // Extract sub_nodes tag values in the old mesh numerotation
                int ind0 = current->parent->sub_nodes[0];
                int ind1 = current->parent->sub_nodes[1];
                int ind2 = current->parent->sub_nodes[2];
                
                // Finding the tag of the sub_nodes using the correspondence table
                int ind[3] = {ind0, ind1, ind2};
                for (int ii=0; ii<3; ii++) {
                    if (ind[ii]!=-1 && ind[ii]>=size_init_node && ind[ii]<len_node+size_init_node) {
                        ind[ii] = tab_node[ind[ii]-size_init_node];
                    }else if(ind[ii]>=len_node+size_init_node){
                        ind[ii] = ind[ii]-n_nodes_removed;
                    }
                }
                
                ind0 = ind[0]; ind1 = ind[1]; ind2 = ind[2];
    
                if (ind0>n_nodes_new || ind1>n_nodes_new || ind2>n_nodes_new) {
                    printf("ERROR IN THE SUB_NODES TAG \n");
                    printing_current(current->parent);
                    exit(1);
                }
                
                x_new[ind0*3+0] = 0.5*(X1+X2); x_new[ind0*3+1] = 0.5*(Y1+Y2); x_new[ind0*3+2] = 0.;
                x_new[ind1*3+0] = 0.5*(X2+X3); x_new[ind1*3+1] = 0.5*(Y2+Y3); x_new[ind1*3+2] = 0.;
                x_new[ind2*3+0] = 0.5*(X3+X1); x_new[ind2*3+1] = 0.5*(Y3+Y1); x_new[ind2*3+2] = 0.;
            }
        }
        
    }else{
        for (int ii=0; ii<current->num_children; ii++) {
            deep_tree_create_data_structure(mesh, &current->children[ii], map_new, x_new, elements_new, boundaries_new, boundary_tags_new, n_nodes_new, n_elements_new, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, 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){
    if (current->children==NULL || current->num_children==0) {
        // DO NOTHING ALREADY PERFORMED AT THE PREVIOUS STEP
    }else{
        if (len_node!=0 && len_elt!=0) {
            deep_tree_map_for_remove(current, len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, n_elements_removed, 1);
        }
        for (int ii=0; ii<current->num_children; ii++) {
            deep_tree_travel_to_get_new_correspondence(&current->children[ii], len_node, tab_node, size_init_node, len_elt, tab_elt, size_init_elt, n_nodes_removed, n_elements_removed);
        }
    }
}

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){
    // Change the tag of the element
    if (bool_tag==1) {
        if (current->tag>=size_init_elt && current->tag<len_elt+size_init_elt) {
            current->tag = tab_elt[current->tag-size_init_elt];
        }else if(current->tag>=len_elt+size_init_elt){
            current->tag = current->tag-n_elements_removed;
        }
    }
    // Change the tag of the nodes
    for (int ii=0; ii<3; ii++) {
        if (current->nodes[ii]>=size_init_node && current->nodes[ii]<len_node+size_init_node) {
            current->nodes[ii] = tab_node[current->nodes[ii]-size_init_node];
        }else if(current->nodes[ii]>=len_node+size_init_node){
            current->nodes[ii] = current->nodes[ii]-n_nodes_removed;
        }
    }
    // Change the tag of the sub_nodes
    for (int ii=0; ii<3; ii++) {
        if (current->sub_nodes[ii]!=-1 && current->sub_nodes[ii]>=size_init_node && current->sub_nodes[ii]<len_node+size_init_node) {
            current->sub_nodes[ii] = tab_node[current->sub_nodes[ii]-size_init_node];
        }else if(current->sub_nodes[ii]>=len_node+size_init_node){
            current->sub_nodes[ii] = current->sub_nodes[ii]-n_nodes_removed;
        }
    }
    // Change the tag of the neighbours
    for (int ii=0; ii<3; ii++) {
        if (current->neighbours[ii]!=-1 && current->neighbours[ii]>=size_init_elt && current->neighbours[ii]<len_elt+size_init_elt) {
            current->neighbours[ii] = tab_elt[current->neighbours[ii]-size_init_elt];
        }else if(current->neighbours[ii]>=len_elt+size_init_elt){
            current->neighbours[ii] = current->neighbours[ii]-n_elements_removed;
        }
    }
}

/*
 This function aims to suppress ONLY the nodes for elements that are pointed by the estimator to be removed.
 */
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){
    if (current->parent!=NULL) {
        if (current->parent->num_children==4) {
            int bool_children_ref = 1;
            int bool_ibd = 1;
            // Check if the brothers of the current element was not refined during the adaptation
            for (int ii=0; ii<4; ii++) {
                bool_children_ref = bool_children_ref & (current->parent->children[ii].children==NULL);
                for (int jj=0; jj<3; jj++) {
                    // Check if the current element is not on a boundary
                    bool_ibd = bool_ibd & (current->parent->children[ii].neighbours[jj]!=-1);
                }
            }
            // Check if the neighbours were refined during the adaptation OR transform into 4
            int bool_nghb[3] = {1, 1, 1}; int bool_tag[3] = {1, 1, 1};
            deep_tree_check_for_coarsening(current, map, bool_nghb, bool_tag, len_elt+size_init_elt);

            // If the current element has not borther that was refined and has no neighbours that was refined, we can remove the children of the current element parent.
            if (bool_children_ref && bool_ibd & bool_nghb[0] & bool_nghb[1] & bool_nghb[2]) {
                // We need to supress the four children of the parent as well as the middle nodes created during previous adaptation
                int i_nghb[3] = {0, 1, 1};
                for (int ii=0; ii<3; ii++) {
                    if (bool_tag[ii]) {
                        int i_node = current->parent->sub_nodes[ii];
                        if (i_node!=-1 && i_node>=size_init_node && i_node<len_node+size_init_node) {
                            tab_node[i_node-size_init_node] = i_node;
                            if (i_node==1158) {
                                printing_current(current->parent);
                            }
                            current->parent->sub_nodes[ii] = -1;
                        }
                        // Transfer the information to the neighbours
                        int index_nghb = current->parent->children[ii].neighbours[i_nghb[ii]];
                        Tree_cell *nghb = map[index_nghb];
                        if (i_node==1158) {
                            printing_current(nghb);
                        }
                        if (nghb->parent!=NULL) {
                            for (int jj=0; jj<3; jj++) {
                                if (nghb->parent->sub_nodes[jj]==i_node) {
                                    nghb->parent->sub_nodes[jj]=-1;
                                }
                            }
                        }
                    } // else the neighbours was transformed into 4
                }
            }
        }
    }
}

/* 
 This function checks if the neighbours of the current elements have some children. It also look after possible transformation through the bool_tag table.
 */
void deep_tree_check_for_coarsening(Tree_cell *current, Tree_cell **map, int bool_nghb[3], int bool_tag[3], int n_elements){
    // Check for the neighbour 0 of the child 0, the neighbour 1 of the child 1 and the neighbour 1 of the child 2
    int i_nghb[3] = {0,1,1};
    for (int ii=0; ii<3; ii++) {
        int index = current->parent->children[ii].neighbours[i_nghb[ii]];
        int bool_def_nghb = 0;
        int bool_sub_node = 1;
        if (index!=-1) {
            Tree_cell *nghb = map[index];
            if (nghb->parent!=NULL) {
                for (int jj=0; jj<nghb->parent->num_children; jj++) {
                    bool_nghb[ii] = bool_nghb[ii] && (nghb->parent->children[jj].children==NULL);
                    bool_def_nghb = bool_def_nghb || (nghb->parent->children[jj].bool_def==1);
                    for (int ii=0; ii<3; ii++) {
                        bool_sub_node = bool_sub_node && (nghb->parent->children[jj].sub_nodes[ii]==-1);
                    }
                    bool_tag[ii] = bool_tag[ii] && (nghb->parent->children[jj].tag<n_elements);
                }
            }
        }
        bool_nghb[ii] = bool_nghb[ii] && bool_def_nghb && bool_sub_node;
    }
    
}

void deep_tree_transform4_into_3(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    int edge = 0;
    // Find the edge that won't be cut
    for (int ii=0; ii<3; ii++) {
        if (parent->sub_nodes[ii]==-1) {
            edge = ii;
        }
    }
    
    int save_bdies[4]; int save_phys[4];
    // Save the tag of the possible boundaries and save the physical value of these boundaries
    if (edge==0) {
        save_bdies[0] = parent->children[1].bdies[1]; save_bdies[1] = parent->children[2].bdies[0];
        save_bdies[2] = parent->children[2].bdies[1]; save_bdies[3] = parent->children[0].bdies[2];
        save_phys[0] = parent->children[1].phys[1]; save_phys[1] = parent->children[2].phys[0];
        save_phys[2] = parent->children[2].phys[1]; save_phys[3] = parent->children[0].phys[2];
    }else if (edge==1){
        save_bdies[0] = parent->children[2].bdies[1]; save_bdies[1] = parent->children[0].bdies[2];
        save_bdies[2] = parent->children[0].bdies[0]; save_bdies[3] = parent->children[1].bdies[0];
        save_phys[0] = parent->children[2].phys[1]; save_phys[1] = parent->children[0].phys[2];
        save_phys[2] = parent->children[0].phys[0]; save_phys[3] = parent->children[1].phys[0];
    }else{
        save_bdies[0] = parent->children[0].bdies[0]; save_bdies[1] = parent->children[1].bdies[0];
        save_bdies[2] = parent->children[1].bdies[1]; save_bdies[3] = parent->children[2].bdies[0];
        save_phys[0] = parent->children[0].phys[0]; save_phys[1] = parent->children[1].phys[0];
        save_phys[2] = parent->children[1].phys[1]; save_phys[3] = parent->children[2].phys[0];
    }
    
    // Save the tag of the 4 children
    int save_tag[4] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag, parent->children[3].tag}; // they are already in ascending order

    // Start to remove the fourth child
    parent->num_children = 3;
    initialize_children(parent, 1);
    int a=0;
    deep_tree_cut_into_3(parent, &a); // we do not create 2 elements !!!!
    for (int ii=0; ii<3; ii++) {
        parent->children[ii].tag = save_tag[ii];
    }
    for (int ii=0; ii<3; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].bdies[jj]=-1;
            parent->children[ii].phys[jj]=-1;
            parent->children[ii].neighbours[jj]=-1;
        }
    }
    parent->children[3].tag = -1;
    parent->children[3].parent = NULL;
    for (int ii=0; ii<3; ii++) {
        parent->children[3].nodes[ii] = -1;
    }
    // Transfer information about the removed elements in the correspondence table
    tab_elt[save_tag[3]-size_init_elt] = save_tag[3];
    
    // Transfer the informations about the boundaries
    parent->children[1].bdies[0] = save_bdies[0]; parent->children[1].phys[0] = save_phys[0];
    parent->children[2].bdies[0] = save_bdies[1]; parent->children[2].phys[0] = save_phys[1];
    parent->children[2].bdies[1] = save_bdies[2]; parent->children[2].phys[1] = save_phys[2];
    parent->children[0].bdies[1] = save_bdies[3]; parent->children[0].phys[1] = save_phys[3];
    
}

void deep_tree_transform4_into_2(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    int edge = 0;
    
    for (int ii=0; ii<3; ii++) {
        if (parent->sub_nodes[ii]!=-1) {
            edge = ii;
        }
    }
    
    int save_bdies[2]; int save_phys[2];
    if (edge==0) {
        save_bdies[0] = parent->children[0].bdies[0]; save_bdies[1] = parent->children[1].bdies[0];
        save_phys[0] = parent->children[0].phys[0]; save_phys[1] = parent->children[1].phys[0];
    }else if (edge==1){
        save_bdies[0] = parent->children[1].bdies[1]; save_bdies[1] = parent->children[2].bdies[0];
        save_phys[0] = parent->children[1].phys[1]; save_phys[1] = parent->children[2].phys[0];
    }else{
        save_bdies[0] = parent->children[2].bdies[1]; save_bdies[1] = parent->children[0].bdies[2];
        save_phys[0] = parent->children[2].phys[1]; save_phys[1] = parent->children[0].phys[2];
    }
    // Save the tag of the 4 children
    int save_tag[4] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag, parent->children[3].tag}; // they are already in ascending order

    // Start to remove the third and fourth children
    parent->num_children = 2;
    initialize_children(parent, 2);
    int a = 0;
    deep_tree_cut_into_2(parent, &a); // we do not create 1 element
    for (int ii=0; ii<2; ii++) {
        parent->children[ii].tag = save_tag[ii];
    }
    
    for (int ii=0; ii<2; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].bdies[jj]=-1;
            parent->children[ii].phys[jj]=-1;
            parent->children[ii].neighbours[jj]=-1;
        }
    }
    for (int ii=2; ii<4; ii++) {
        tab_elt[save_tag[ii]-size_init_elt] = save_tag[ii];
        parent->children[ii].tag = -1;
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = -1;
        }
        parent->children[ii].parent = NULL;
    }
    
    // Transfer the informations about the boundaries
    parent->children[0].bdies[0] = save_bdies[0]; parent->children[1].bdies[0] = save_bdies[1];
    parent->children[0].phys[0] = save_phys[0]; parent->children[1].phys[0] = save_phys[1];
    
}
void deep_tree_transform4_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    // we do not need to treat boundaries in this case
    int save_tag[4] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag, parent->children[3].tag}; // they are already in ascending order

    parent->num_children = 0;
    parent->tag = save_tag[0];
    initialize_children(parent, 4);
    for (int ii=0; ii<4; ii++) {
        parent->children[ii].tag = -1;
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = -1;
        }
        parent->children[ii].parent = NULL;
    }
    
    for (int ii=1; ii<4; ii++) {
        tab_elt[save_tag[ii]-size_init_elt] = save_tag[ii];
    }
}

void deep_tree_transform3_into_2(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    int edge = 0; // the cutting edge for the 2-cut
    
    for (int ii=0; ii<3; ii++) {
        if (parent->sub_nodes[ii]!=-1) {
            edge = ii;
        }
    }
    int edge_cut = 0; // the none cutting edge for the 3-cut (previously)
    for (int ii=0; ii<3; ii++) {
        if (parent->nodes[ii] == parent->children[0].nodes[2] && parent->nodes[(ii+1)%3] == parent->children[0].nodes[0]) {
            edge_cut = ii;
        }
    }
    int save_bdies[6] = {parent->children[0].bdies[1], parent->children[2].bdies[1], parent->children[2].bdies[0],
                         parent->children[0].bdies[2], parent->children[1].bdies[0], parent->children[2].bdies[0]};
    int save_phys[6] = {parent->children[0].phys[1], parent->children[2].phys[1], parent->children[2].phys[0],
                        parent->children[0].phys[2], parent->children[1].phys[0], parent->children[2].phys[0]};
    int save_tag[3] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag};
    
    // Start removing the third element
    parent->num_children = 2;
    initialize_children(parent, 2);
    int a = 0;
    deep_tree_cut_into_2(parent, &a); // we do not create 1 element
    for (int ii=0; ii<2; ii++) {
        parent->children[ii].tag = save_tag[ii];
    }
    for (int ii=0; ii<2; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].bdies[jj]=-1;
            parent->children[ii].phys[jj]=-1;
            parent->children[ii].neighbours[jj]=-1;
        }
    }
    for (int ii=2; ii<3; ii++) {
        tab_elt[save_tag[ii]-size_init_elt] = save_tag[ii];
        parent->children[ii].tag = -1;
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = -1;
        }
        parent->children[ii].parent = NULL;
    }
    // Transfer the informations about the boundaries
    if ((edge+2)%3 == edge_cut) {
        parent->children[1].bdies[0] = save_bdies[0]; parent->children[1].phys[0] = save_phys[0];
        parent->children[0].bdies[0] = save_bdies[1]; parent->children[0].phys[0] = save_phys[1];
        parent->children[1].bdies[1] = save_bdies[2]; parent->children[1].phys[1] = save_phys[2];
    }else if((edge+1)%3 == edge_cut){
        parent->children[0].bdies[2] = save_bdies[3]; parent->children[0].phys[2] = save_phys[3];
        parent->children[0].bdies[0] = save_bdies[4]; parent->children[0].phys[0] = save_phys[4];
        parent->children[1].bdies[0] = save_bdies[5]; parent->children[1].phys[0] = save_phys[5];
    }
    
    
}
void deep_tree_transform3_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    int edge = 0;
    for (int ii=0; ii<3; ii++) {
        if (parent->nodes[ii] == parent->children[0].nodes[2] && parent->nodes[(ii+1)%3] == parent->children[0].nodes[0]) {
            edge = ii;
        }
    }
    int save_bdies = parent->children[0].bdies[2];
    int save_phys = parent->children[0].phys[2];
    int save_tag[3] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag};

    parent->num_children = 0;
    parent->tag = save_tag[0];
    initialize_children(parent, 4);
    for (int ii=0; ii<3; ii++) {
        parent->children[ii].tag = -1;
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = -1;
        }
        parent->children[ii].parent = NULL;
    }
    parent->bdies[edge] = save_bdies;
    parent->phys[edge] = save_phys;
    
    for (int ii=1; ii<3; ii++) {
        tab_elt[save_tag[ii]-size_init_elt] = save_tag[ii];
    }
}

void deep_tree_transform2_into_1(Tree_cell *parent, int len_elt, int *tab_elt, int size_init_elt){
    int edge = 0;
    for (int ii=0; ii<3; ii++) {
        int bool_edge0 = (parent->nodes[ii] == parent->children[0].nodes[0]) && (parent->nodes[(ii+2)%3] == parent->children[0].nodes[2]);
        int bool_edge1 = (parent->nodes[(ii+1)%3] == parent->children[1].nodes[1]) && (parent->nodes[(ii+2)%3] == parent->children[1].nodes[2]);
        if (bool_edge0 && bool_edge1) {
            edge = ii;
        }
    }
    int save_bdies[2] = {parent->children[0].bdies[2], parent->children[1].bdies[1]};
    int save_phys[2] = {parent->children[0].phys[2], parent->children[1].phys[1]};
    int save_tag[2] = {parent->children[0].tag, parent->children[1].tag};

    parent->num_children = 0;
    parent->tag = save_tag[0];
    initialize_children(parent, 4);
    for (int ii=0; ii<2; ii++) {
        parent->children[ii].tag = -1;
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = -1;
        }
        parent->children[ii].parent = NULL;
    }
    
    // Transfer the information about the boundaries
    parent->bdies[(edge+2)%3] = save_bdies[0]; parent->phys[(edge+2)%3] = save_phys[0];
    parent->bdies[(edge+1)%3] = save_bdies[1]; parent->phys[(edge+1)%3] = save_phys[1];
    
    tab_elt[save_tag[1]-size_init_elt] = save_tag[1];
}

void deep_tree_cut_into_4(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd){
    // SET THE BOOLEAN FOR REFINEMENT TO 0
    current->bool_ref = 0;
    // CREATE INDEX FOR THE THREE NEW NODES
    for (int ii=0; ii<3; ii++) {
        if (current->sub_nodes[ii]==-1) {
            current->sub_nodes[ii] = *index;
            *index+=1;
            if (current->neighbours[ii]==-1) {
                if (ii==0) {
                    current->children[0].bdies[0] = current->bdies[0];
                    current->children[0].phys[0] = current->phys[0];
                    current->children[1].bdies[0] = *ibd;
                    current->children[1].phys[0] = current->phys[0];
                }else if(ii==1){
                    current->children[1].bdies[1] = current->bdies[1];
                    current->children[1].phys[1] = current->phys[1];
                    current->children[2].bdies[0] = *ibd;
                    current->children[2].phys[0] = current->phys[1];
                }else{
                    current->children[2].bdies[1] = current->bdies[2];
                    current->children[2].phys[1] = current->phys[2];
                    current->children[0].bdies[2] = *ibd;
                    current->children[0].phys[2] = current->phys[2];
                }
                *ibd+=1;
            }
        }
    }
    // CREATE THE FOUR CHILDREN AND TRANSFERT INFORMATION ABOUT THE NEIGHBOURS
    int nodes[6] = {current->nodes[0], current->sub_nodes[0], current->nodes[1], current->sub_nodes[1], current->nodes[2], current->sub_nodes[2]};
    int tab[4][3] = {{0,1,5}, {1,2,3}, {3,4,5}, {1,3,5}};
    
    for (int ii=0; ii<4; ii++) {
        for (int jj=0; jj<3; jj++) {
            current->children[ii].nodes[jj] = nodes[tab[ii][jj]];
        }
    }
    
    // TRANSFERT INFORMATION ABOUT THE NODES TAG
    if (current->tag!=-1) {
        for (int ii=0; ii<3; ii++) {
            int buf = current->neighbours[ii];
            if (buf>-1 && current->sub_nodes[ii]!=-1) {
                if (map[buf]->children==NULL) {
                    for (int jj=0; jj<3; jj++) {
                        if (map[buf]->neighbours[jj]==current->tag) {
                            map[buf]->sub_nodes[jj]=current->sub_nodes[ii];
                        }
                    }
                }
            }
        }
    }
    
    // GIVES NEW TAG TO CHILDREN
    current->children[0].tag = current->tag;
    current->tag = -1; // BECAUSE IT BECOMES A PARENT
    current->children[1].tag = *elt; *elt+=1;
    current->children[2].tag = *elt; *elt+=1;
    current->children[3].tag = *elt; *elt+=1;
}

void deep_tree_cut_into_2(Tree_cell *current, int *elt){
    for (int ii=0; ii<3; ii++) {
        if (current->sub_nodes[ii]!=-1) {
            current->children[0].nodes[0] = current->nodes[ii];
            current->children[0].nodes[1] = current->sub_nodes[ii];
            current->children[0].nodes[2] = current->nodes[(ii+2)%3];
            current->children[0].bdies[2] = current->bdies[(ii+2)%3];
            current->children[0].phys[2] = current->phys[(ii+2)%3];
            
            current->children[1].nodes[0] = current->sub_nodes[ii];
            current->children[1].nodes[1] = current->nodes[(ii+1)%3];
            current->children[1].nodes[2] = current->nodes[(ii+2)%3];
            current->children[1].bdies[1] = current->bdies[(ii+1)%3];
            current->children[1].phys[1] = current->phys[(ii+1)%3];
        }
    }
    current->children[0].tag = current->tag;
    current->tag = -1; // BECAUSE IT BECOMES A PARENT
    current->children[1].tag = *elt; *elt+=1;
}

void deep_tree_cut_into_3(Tree_cell *current, int *elt){
    for (int ii=0; ii<3; ii++) {
        if (current->sub_nodes[ii]==-1) {
            
            current->children[0].nodes[0] = current->nodes[(ii+1)%3];
            current->children[0].nodes[1] = current->sub_nodes[(ii+2)%3];
            current->children[0].nodes[2] = current->nodes[ii];
            current->children[0].bdies[2] = current->bdies[ii];
            current->children[0].phys[2] = current->phys[ii];
            
            current->children[2].nodes[0] = current->sub_nodes[(ii+1)%3];
            current->children[2].nodes[1] = current->nodes[(ii+2)%3];
            current->children[2].nodes[2] = current->sub_nodes[(ii+2)%3];
            
            current->children[1].nodes[0] = current->nodes[(ii+1)%3];
            current->children[1].nodes[1] = current->sub_nodes[(ii+1)%3];
            current->children[1].nodes[2] = current->sub_nodes[(ii+2)%3];
        }
    }
    current->children[0].tag = current->tag;
    current->tag = -1; // BECAUSE IT BECOMES A PARENT
    current->children[1].tag = *elt; *elt+=1;
    current->children[2].tag = *elt; *elt+=1;
}

void deep_tree_transform2_into_4(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd){
    Tree_cell *parent = current->parent;
    
    // SAVE THE TAG OF THE CUTTING EDGE (0,1,2)
    int e_tag;
    for (int ii=0; ii<3; ii++) {
        if (parent->sub_nodes[ii]!=-1) {
            e_tag = ii;
        }
    }
    
    // SAVE SUB_NODES, BDIES AND PHYS FROM THE CHILDREN (WE ONLY TRUST ABOUT LOCAL INFORMATION)
    int sub_nodes_chld[2][3] = {{parent->children[0].sub_nodes[0], parent->children[0].sub_nodes[1], parent->children[0].sub_nodes[2]},
        {parent->children[1].sub_nodes[0], parent->children[1].sub_nodes[1], parent->children[1].sub_nodes[2]}};
    
    int bdies_chld[2][3] = {{parent->children[0].bdies[0], parent->children[0].bdies[1], parent->children[0].bdies[2]},
        {parent->children[1].bdies[0], parent->children[1].bdies[1], parent->children[1].bdies[2]}};
    
    int phys_chld[2][3] = {{parent->children[0].phys[0], parent->children[0].phys[1], parent->children[0].phys[2]},
            {parent->children[1].phys[0], parent->children[1].phys[1], parent->children[1].phys[2]}};
    
    for (int ii=0; ii<2; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].bdies[jj] = -1;
            parent->children[ii].phys[jj] = -1;
            parent->children[ii].sub_nodes[jj] = -1;
        }
    }
    
    // SAVE THE TAG OF THE TWO CHILDREN
    int save_tag[2] = {parent->children[0].tag, parent->children[1].tag};
    
    // WE CHECK IF THERE IS ALREADY OR NOT A SUB_NODE ON THE TWO NONE CUTTING EDGES
    if (sub_nodes_chld[0][2]==-1) {
        parent->sub_nodes[(e_tag+2)%3] = *index;
        if (parent->children[0].neighbours[2]==-1) {
            if ((e_tag+2)%3==0) {
                parent->children[0].bdies[0] = bdies_chld[0][2];
                parent->children[0].phys[0] = phys_chld[0][2];
                parent->children[1].bdies[0] = *ibd;
                parent->children[1].phys[0] = phys_chld[0][2];
            }else if((e_tag+2)%3==1){
                parent->children[1].bdies[1] = bdies_chld[0][2];
                parent->children[1].phys[1] = phys_chld[0][2];
                parent->children[2].bdies[0] = *ibd;
                parent->children[2].phys[0] = phys_chld[0][2];
            }else{
                parent->children[2].bdies[1] = bdies_chld[0][2];
                parent->children[2].phys[1] = phys_chld[0][2];
                parent->children[0].bdies[2] = *ibd;
                parent->children[0].phys[2] = phys_chld[0][2];
            }
            *ibd+=1;
        }
        int buf = parent->children[0].neighbours[2];
        if (buf>-1){
            if(map[buf]->children==NULL) {
                for (int jj=0; jj<3; jj++){
                    if (map[buf]->neighbours[jj]==save_tag[0]) {
                        map[buf]->sub_nodes[jj]=*index;
                    }
                }
            }
        }
        *index+=1;
    }else{
        parent->sub_nodes[(e_tag+2)%3] = sub_nodes_chld[0][2];
    }
    if (sub_nodes_chld[1][1]==-1) {
        parent->sub_nodes[(e_tag+1)%3] = *index;
        if (parent->children[1].neighbours[1]==-1) {
            if ((e_tag+1)%3==0) {
                parent->children[0].bdies[0] = bdies_chld[1][1];
                parent->children[0].phys[0] = phys_chld[1][1];
                parent->children[1].bdies[0] = *ibd;
                parent->children[1].phys[0] = phys_chld[1][1];
            }else if((e_tag+1)%3==1){
                parent->children[1].bdies[1] = bdies_chld[1][1];
                parent->children[1].phys[1] = phys_chld[1][1];
                parent->children[2].bdies[0] = *ibd;
                parent->children[2].phys[0] = phys_chld[1][1];
            }else{
                parent->children[2].bdies[1] = bdies_chld[1][1];
                parent->children[2].phys[1] = phys_chld[1][1];
                parent->children[0].bdies[2] = *ibd;
                parent->children[0].phys[2] = phys_chld[1][1];
            }
            *ibd+=1;
        }
        int buf = parent->children[1].neighbours[1];
        if (buf>-1){
            if(map[buf]->children==NULL) {
                for (int jj=0; jj<3; jj++){
                    if (map[buf]->neighbours[jj]==save_tag[1]) {
                        map[buf]->sub_nodes[jj]=*index;
                    }
                }
            }
        }
        *index+=1;
    }else{
        parent->sub_nodes[(e_tag+1)%3] = sub_nodes_chld[1][1];
    }
    
    // CREATE THE NEW PARENT AND ITS FOUR CHILDREN
    parent->num_children = 4;
    initialize_children_bis(parent);

    for (int ii=0; ii<2; ii++) {
        parent->children[ii].tag = save_tag[ii];
    }
    parent->children[2].tag = *elt; *elt+=1;
    parent->children[3].tag = *elt; *elt+=1;
    
    // CUT THE PARENT INTO FOUR CHILDREN
    int nodes[6] = {parent->nodes[0], parent->sub_nodes[0], parent->nodes[1], parent->sub_nodes[1], parent->nodes[2], parent->sub_nodes[2]};
    int tab[4][3] = {{0,1,5}, {1,2,3}, {3,4,5}, {1,3,5}};
    
    for (int ii=0; ii<4; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = nodes[tab[ii][jj]];
        }
    }
    
    // LOOK FOR POSSIBLE REFINEMENT OF THE CHILDREN 0, 1 and 2 (those one which are on the borders of the parent element)
    if (e_tag == 0) {
        parent->children[0].sub_nodes[0] = sub_nodes_chld[0][0];
        parent->children[1].sub_nodes[0] = sub_nodes_chld[1][0];
    }else if(e_tag == 1){
        parent->children[1].sub_nodes[1] = sub_nodes_chld[0][0];
        parent->children[2].sub_nodes[0] = sub_nodes_chld[1][0];
    }else{
        parent->children[2].sub_nodes[1] = sub_nodes_chld[0][0];
        parent->children[0].sub_nodes[2] = sub_nodes_chld[1][0];
    }
    
    // SMOOTH THE REFINEMENT OF THE CHILDREN BY CUTTING INTO 2 OR 3
    for (int ii=0; ii<2; ii++) {
        int bool_one_node = (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]!=-1);
        int bool_two_nodes = (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]!=-1) || (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]!=-1);
        if (bool_one_node) {
            // CUT_INTO_2
            parent->children[ii].children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
            parent->children[ii].num_children = 2;
            initialize_children(&parent->children[ii], 2);
            deep_tree_cut_into_2(&parent->children[ii], elt);
        }else if(bool_two_nodes){
            // CUT_INTO_3
            parent->children[ii].children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
            parent->children[ii].num_children = 3;
            initialize_children(&parent->children[ii], 1);
            deep_tree_cut_into_3(&parent->children[ii], elt);
        }
    }
}

void deep_tree_transform3_into_4(Tree_cell *current, Tree_cell **map, int *index, int *elt, int *ibd){
    Tree_cell *parent = current->parent;
    // SAVE THE LOCAL TAG OF THE NONE CUTTING EDGE (0,1,2)
    int e_tag;
    for (int ii=0; ii<3; ii++) {
        if (parent->sub_nodes[ii]==-1) {
            e_tag = ii;
        }
    }
    // SAVE SUB_NODES, BDIES AND PHYS FROM THE THREE CHILDREN
    int sub_nodes_chld[3][3] = {{parent->children[0].sub_nodes[0], parent->children[0].sub_nodes[1], parent->children[0].sub_nodes[2]},
                                {parent->children[1].sub_nodes[0], parent->children[1].sub_nodes[1], parent->children[1].sub_nodes[2]},
                                {parent->children[2].sub_nodes[0], parent->children[2].sub_nodes[1], parent->children[2].sub_nodes[2]}};
    
    int bdies_chld[3][3] = {{parent->children[0].bdies[0], parent->children[0].bdies[1], parent->children[0].bdies[2]},
                            {parent->children[1].bdies[0], parent->children[1].bdies[1], parent->children[1].bdies[2]},
                            {parent->children[2].bdies[0], parent->children[2].bdies[1], parent->children[2].bdies[2]}};
    
    int phys_chld[3][3] = {{parent->children[0].phys[0], parent->children[0].phys[1], parent->children[0].phys[2]},
                           {parent->children[1].phys[0], parent->children[1].phys[1], parent->children[1].phys[2]},
                           {parent->children[2].phys[0], parent->children[2].phys[1], parent->children[2].phys[2]}};
    
    for (int ii=0; ii<3; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].bdies[jj] = -1;
            parent->children[ii].phys[jj] = -1;
            parent->children[ii].sub_nodes[jj] = -1;
        }
    }
    // SAVE THE TAG OF THE THREE CHILDREN
    int save_tag[3] = {parent->children[0].tag, parent->children[1].tag, parent->children[2].tag};
    
    // WE CHECK IF THERE IS ALREADY OR NOT A SUB_NODE ON THE THIRD EDGE OF THE FIRST CHILD
    if (sub_nodes_chld[0][2]==-1) {
        parent->sub_nodes[e_tag] = *index;
        if (parent->children[0].neighbours[2]==-1) {
            if (e_tag==0) {
                parent->children[0].bdies[0] = bdies_chld[0][2];
                parent->children[0].phys[0] = phys_chld[0][2];
                parent->children[1].bdies[0] = *ibd;
                parent->children[1].phys[0] = phys_chld[0][2];
            }else if(e_tag==1){
                parent->children[1].bdies[1] = bdies_chld[0][2];
                parent->children[1].phys[1] = phys_chld[0][2];
                parent->children[2].bdies[0] = *ibd;
                parent->children[2].phys[0] = phys_chld[0][2];
            }else{
                parent->children[2].bdies[1] = bdies_chld[0][2];
                parent->children[2].phys[1] = phys_chld[0][2];
                parent->children[0].bdies[2] = *ibd;
                parent->children[0].phys[2] = phys_chld[0][2];
            }
            *ibd+=1;
        }
        int buf = parent->children[0].neighbours[2];
        if (buf>-1){
            if(map[buf]->children==NULL) {
                for (int jj=0; jj<3; jj++) {
                    if (map[buf]->neighbours[jj]==parent->children[0].tag) {
                        map[buf]->sub_nodes[jj]=*index;
                    }
                }
            }
        }
        *index+=1;
    }else{
        parent->sub_nodes[e_tag] = sub_nodes_chld[0][2];
    }
    
    // CREATE THE NEW PARENT AND ITS FOUR CHILDREN
    parent->num_children = 4;
    initialize_children_bis(parent);
    
    for (int ii=0; ii<3; ii++) {
        parent->children[ii].tag = save_tag[ii];
    }
    parent->children[3].tag = *elt; *elt+=1;
    
    
    // CUT THE PARENT INTO FOUR CHILDREN
    int nodes[6] = {parent->nodes[0], parent->sub_nodes[0], parent->nodes[1], parent->sub_nodes[1], parent->nodes[2], parent->sub_nodes[2]};
    int tab[4][3] = {{0,1,5}, {1,2,3}, {3,4,5}, {1,3,5}};
    
    for (int ii=0; ii<4; ii++) {
        for (int jj=0; jj<3; jj++) {
            parent->children[ii].nodes[jj] = nodes[tab[ii][jj]];
        }
    }
    
    // LOOK FOR POSSIBLE REFINEMENT OF THREE OF THE CHILREN, THE CHILD 0, 1 AND 2 (NOTHING ABOUT THE THIRD ONE)
    if (e_tag==0) {
        parent->children[0].sub_nodes[2] = sub_nodes_chld[0][1];
        parent->children[1].sub_nodes[1] = sub_nodes_chld[1][0];
        parent->children[2].sub_nodes[0] = sub_nodes_chld[2][0];
        parent->children[2].sub_nodes[1] = sub_nodes_chld[2][1];
    }else if(e_tag==1){
        parent->children[1].sub_nodes[1] = sub_nodes_chld[0][1];
        parent->children[2].sub_nodes[1] = sub_nodes_chld[1][0];
        parent->children[0].sub_nodes[2] = sub_nodes_chld[2][0];
        parent->children[0].sub_nodes[0] = sub_nodes_chld[2][1];
    }else{
        parent->children[2].sub_nodes[0] = sub_nodes_chld[0][1];
        parent->children[0].sub_nodes[0] = sub_nodes_chld[1][0];
        parent->children[1].sub_nodes[0] = sub_nodes_chld[2][0];
        parent->children[1].sub_nodes[1] = sub_nodes_chld[2][1];
    }
    
    // SMOOTH THE MESH BY CUTING INTO 2 OR 3 (ONLY FOR THE CHILDREN 0,1 and 2)
    for (int ii=0; ii<2; ii++) {
        int bool_one_node = (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]!=-1);
        int bool_two_nodes = (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]==-1) || (parent->children[ii].sub_nodes[0]==-1 && parent->children[ii].sub_nodes[1]!=-1 && parent->children[ii].sub_nodes[2]!=-1) || (parent->children[ii].sub_nodes[0]!=-1 && parent->children[ii].sub_nodes[1]==-1 && parent->children[ii].sub_nodes[2]!=-1);
        if (bool_one_node) {
            // CUT_INTO_2
            parent->children[ii].children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
            parent->children[ii].num_children = 2;
            initialize_children(&parent->children[ii], 2);
            deep_tree_cut_into_2(&parent->children[ii], elt);
        }else if(bool_two_nodes){
            // CUT_INTO_3
            parent->children[ii].children = (Tree_cell *) malloc(4*sizeof(Tree_cell));
            parent->children[ii].num_children = 3;
            initialize_children(&parent->children[ii], 1);
            deep_tree_cut_into_3(&parent->children[ii], elt);
        }
    }
}

void initialize_children_bis(Tree_cell *current){
    for (int ii=0; ii<current->num_children; ii++) {
        current->children[ii].parent=current;
        current->children[ii].children=NULL;
        current->children[ii].num_children = 0;
        current->children[ii].bool_ref = 0;
        current->children[ii].bool_def = 0;
        for (int jj=0; jj<3; jj++) {
            current->children[ii].neighbours[jj]=-1;
            current->children[ii].sub_nodes[jj]=-1;
        }
    }
}

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){
    int count = 0;
    for (int ii=0; ii<len_node; ii++) {
        if (tab_node[ii]==-1) {
            tab_node[ii] = ii+size_init_node-count;
        }else{
            if (tab_node[ii]==ii+size_init_node) {
                tab_node[ii] = -1;
                count++;
            }else{
                tab_node[ii] = ii+size_init_node-count;
            }
        }
    }
    *index_remove = count;
    count = 0;
    for (int ii=0; ii<len_elt; ii++) {
        if (tab_elt[ii]==-1) {
            tab_elt[ii] = ii+size_init_elt-count;
        }else{
            if (tab_elt[ii]==ii+size_init_elt) {
                tab_elt[ii] = -1;
                count++;
            }else{
                tab_elt[ii]= ii+size_init_elt-count;
            }
            
        }
    }
    *elt_remove = count;
}

void printing_current(Tree_cell *current){
    printf("----------TAG %d---------- \n", current->tag);
    printf("NODES : (%d, %d, %d)\n", current->nodes[0], current->nodes[1], current->nodes[2]);
    printf("SUB_NODES : (%d, %d, %d)\n", current->sub_nodes[0], current->sub_nodes[1], current->sub_nodes[2]);
    printf("NEIGHBOURS : (%d, %d, %d)\n", current->neighbours[0], current->neighbours[1], current->neighbours[2]);
    printf("BDIES : (%d, %d, %d)\n", current->bdies[0], current->bdies[1], current->bdies[2]);
    printf("PHYS : (%d, %d, %d)\n", current->phys[0], current->phys[1], current->phys[2]);
    printf("NUM_CHILDREN %d \n", current->num_children);
    printf("CHILDREN ? %d\n", current->children!=NULL);
    printf("BOOL_REF %d \n", current->bool_ref);
    printf("BOOL_DEF %d \n", current->bool_def);
    printf("BOOL_BOUNDARY %d \n", current->bool_boundary);
    printf("PARENT ? %d \n", current->parent!=NULL);
    if (current->parent!=NULL) {
        printf("PARENT-NUM_CHILDREN  %d \n", current->parent->num_children);
    }
    
}

void printing_post_processing_files(Mesh *mesh, int *elements_new, double *x_new, int *boundaries_new, int n_elements_new, int n_boundaries_new){
    char mesh_filename1[128] = "";
    sprintf(mesh_filename1, "/Users/margauxboxho/migflow/testcases/my_local_drop/New_Mesh.dat");
    
    FILE *file1 = fopen(mesh_filename1, "w");
    if(file1 == NULL){
        printf("Error at opening file %s!\n", mesh_filename1);
        exit(1);
    }
    
    char mesh_filename2[128] = "";
    sprintf(mesh_filename2, "/Users/margauxboxho/migflow/testcases/my_local_drop/Old_Mesh.dat");

    FILE *file2 = fopen(mesh_filename2, "w");
    if(file2 == NULL){
        printf("Error at opening file %s!\n", mesh_filename2);
        exit(1);
    }
    
    char mesh_filename3[128] = "";
    sprintf(mesh_filename3, "/Users/margauxboxho/migflow/testcases/my_local_drop/Edges.dat");
    
    FILE *file3 = fopen(mesh_filename3, "w");
    if(file3 == NULL){
        printf("Error at opening file %s!\n", mesh_filename3);
        exit(1);
    }
    
    // Copy old mesh
    for (int iel=0; iel<mesh->n_elements; iel++) {
        const int *el = &mesh->elements[iel*N_N];
        int tag1 = el[0];
        int tag2 = el[1];
        int tag3 = el[2];
        fprintf(file2, "%d, %d, %f, %f, %d, %f, %f, %d, %f, %f \n", iel, tag1, mesh->x[tag1*3+0], mesh->x[tag1*3+1], tag2, mesh->x[tag2*3+0], mesh->x[tag2*3+1], tag3, mesh->x[tag3*3+0], mesh->x[tag3*3+1]);
    }
    
    // Copy new mesh
    for (int iel=0; iel<n_elements_new; iel++) {
        int tag1 = elements_new[iel*3+0];
        int tag2 = elements_new[iel*3+1];
        int tag3 = elements_new[iel*3+2];
        fprintf(file1, "%d, %d, %f, %f, %d, %f, %f, %d, %f, %f \n", iel, tag1, x_new[tag1*3+0], x_new[tag1*3+1], tag2, x_new[tag2*3+0], x_new[tag2*3+1], tag3, x_new[tag3*3+0], x_new[tag3*3+1]);
    }
    
    // Copy boundaries
    for (int kk=0; kk<n_boundaries_new; kk++) {
        int tag1 = boundaries_new[kk*2+0];
        int tag2 = boundaries_new[kk*2+1];
        fprintf(file3, "%d, %f, %f, %f, %f \n", kk, x_new[tag1*3+0], x_new[tag1*3+1], x_new[tag2*3+0], x_new[tag2*3+1]);
    }
    
    fclose(file1);
    fclose(file2);
    fclose(file3);
}

/*
for (int iel=0; iel<mesh->n_elements; iel++) {
    Tree_cell *current = map[iel];
    if (current->tag==-1) {
        // THE CURRENT ELEMENT IS A PARENT OF 2, 3 OR 4 CHILDREN
        if (current->num_children==4) {
            for (int ii=0; ii<4; ii++) {
                for (int jj=0; jj<3; jj++) {
                    int ibd = current->children[ii].bdies[jj];
                    if (ibd>-1) {
                        count++;
                        boundaries_new[2*ibd+0] = current->children[ii].nodes[jj];
                        boundaries_new[2*ibd+1] = current->children[ii].nodes[(jj+1)%3];
                        boundary_tags_new[ibd] = current->children[ii].phys[jj];
                        //printf("%4d, %d, (%d, %d) \n", ibd, boundary_tags_new[ibd], boundaries_new[2*ibd+0], boundaries_new[2*ibd+1]);
                        //printf("%4d \n", ibd);//, boundary_tags_new[ibd], boundaries_new[2*ibd+0], boundaries_new[2*ibd+1]);
                    }
                }
            }
        }else{
            for (int ii=0; ii<3; ii++) {
                int ibd = current->bdies[ii];
                if (ibd>-1) {
                    count++;
                    boundaries_new[2*ibd+0] = current->nodes[ii];
                    boundaries_new[2*ibd+1] = current->nodes[(ii+1)%3];
                    boundary_tags_new[ibd] = current->phys[ii];
                    //printf("%4d, %d, (%d, %d) \n", ibd, boundary_tags_new[ibd], boundaries_new[2*ibd+0], boundaries_new[2*ibd+1]);
                    //printf("%4d \n", ibd);
                }
            }
        }
    }else{
        // THE CURRENT ELEMENT IS EITHER AND OLD ELEMENT OR A TRANSFORM ELEMENT
        if (current->parent!=NULL) {
            for (int ii=0; ii<current->parent->num_children; ii++) {
                for (int jj=0; jj<3; jj++) {
                    int ibd = current->parent->children[ii].bdies[jj];
                    if (ibd>-1) {
                        boundaries_new[2*ibd+0] = current->parent->children[ii].nodes[jj];
                        boundaries_new[2*ibd+1] = current->parent->children[ii].nodes[(jj+1)%3];
                        boundary_tags_new[ibd] = current->parent->children[ii].phys[jj];
                        //printf("%4d, %d, (%d, %d)\n", ibd, boundary_tags_new[ibd], boundaries_new[2*ibd+0], boundaries_new[2*ibd+1]);
                        //printf("%4d \n", ibd);
                    }
                }
            }
        }else{
            for (int ii=0; ii<3; ii++) {
                int ibd = current->bdies[ii];
                if (ibd>-1) {
                    count++;
                    boundaries_new[2*ibd+0] = current->nodes[ii];
                    boundaries_new[2*ibd+1] = current->nodes[(ii+1)%3];
                    boundary_tags_new[ibd] = current->phys[ii];
                    //printf("%4d, %d, (%d, %d)\n", ibd, boundary_tags_new[ibd], boundaries_new[2*ibd+0], boundaries_new[2*ibd+1]);
                    //printf("%4d \n", ibd);
                }
            }
        }
    }
}*/