Father is stays, mother is merged.
Move node on the edge, 0 not move, 1 move to mother side, 0.5 will be in the middle. 
   37                                             {
   41 
   42  
   44  CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1, 
false, father_edges);
 
   46  CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1, 
false, mother_edges);
 
   47 
   48  
   50  CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3, 
false, father_tets);
 
   52  CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 3, 
false, mother_tets);
 
   53  if (tets_ptr != NULL) {
   54    mother_tets = intersect(mother_tets, *tets_ptr);
   55    father_tets = intersect(father_tets, *tets_ptr);
   56  }
   57 
   58  
   60  common_edge = intersect(father_edges, mother_edges);
   61  if (tets_ptr != NULL) {
   62    Range tets = unite(father_tets, mother_tets);
 
   64    CHKERR m_field.get_moab().get_adjacencies(tets, 1, 
false, tets_edges,
 
   65                                              moab::Interface::UNION);
   66    common_edge = intersect(common_edge, tets_edges);
   67    father_edges = intersect(father_edges, tets_edges);
   68    mother_edges = intersect(mother_edges, tets_edges);
   69  }
   70 
   71  
   74            "no common edge between nodes");
   75  } else if (common_edge.empty()) {
   77    if (tets_ptr != NULL) {
   78      seed_tets.merge(*tets_ptr);
   79    }
   80    out_tets = seed_tets;
   83  }
   84 
   85  
   87  CHKERR m_field.get_moab().get_adjacencies(common_edge, 3, 
true, edge_tets);
 
   88  
   89  if (tets_ptr != NULL) {
   90    edge_tets = intersect(edge_tets, *tets_ptr);
   91  }
   92  
   93  mother_tets = subtract(mother_tets, edge_tets);
   94  father_tets = subtract(father_tets, edge_tets);
   95 
   97                               const int num_nodes) {
  100      CHKERR m_field.get_moab().get_coords(conn, num_nodes, &coords[0]);
 
  101    } else {
  102      CHKERR m_field.get_moab().tag_get_data(
th, conn, num_nodes, &coords[0]);
 
  103    }
  104    return coords;
  105  };
  106 
  107  auto get_tensor = [](
VectorDouble &coords, 
const int shift) {
 
  109        &coords[shift], &coords[shift + 1], &coords[shift + 2]);
  110  };
  111 
  112  
  114  if (move > 0) {
  117    auto t_n0 = get_tensor(coords, 0);
  118    auto t_n1 = get_tensor(coords, 3);
  119    t_move(
i) = t_n0(
i) + move * (t_n1(
i) - t_n0(
i));
 
  120  }
  121 
  122  if (line_search > 0) {
  123    Range check_tests = unite(father_tets, mother_tets);
 
  125  }
  126 
  127  if (only_if_improve_quality) {
  128    Range check_tests = mother_tets;
 
  129    
  130    if (move > 0 || line_search) {
  131      check_tests.merge(father_tets);
  132    }
  133 
  134    double min_quality0 = 1;
  139    double min_quality = 1;
  141                      ((move > 0) || line_search) ? &t_move(0) : NULL,
  143    if (min_quality < min_quality0) {
  144      if (tets_ptr != NULL) {
  145        out_tets = *tets_ptr;
  146      }
  149    }
  150  }
  151 
  152  
  153  if (move > 0 || line_search) {
  155      CHKERR m_field.get_moab().set_coords(&father, 1, &t_move(0));
 
  156    } else {
  157      CHKERR m_field.get_moab().tag_set_data(
th, &father, 1, &t_move(0));
 
  158    }
  159  }
  160 
  162                             int *ret_num_nodes = nullptr) {
  163    int num_nodes;
  165    CHKERR m_field.get_moab().get_connectivity(ent, conn, num_nodes, 
true);
 
  166    if (ret_num_nodes)
  167      *ret_num_nodes = num_nodes;
  168    return conn;
  169  };
  170 
  171  auto create_tet = [
this, &m_field](
const EntityHandle *new_conn,
 
  175    CHKERR m_field.get_moab().get_adjacencies(new_conn, 4, 3, 
false, tets);
 
  176    bool tet_found = false;
  177    for (auto it_tet : tets) {
  179      int num_nodes;
  180      CHKERR m_field.get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
 
  181                                                 true);
  182      const EntityHandle *p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
 
  183      if (p != tet_conn + 4) {
  184        int s = std::distance(tet_conn, p);
  186        for (; 
n != 4; ++
n) {
 
  187          const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
  188          if (tet_conn[idx[s + 
n]] != new_conn[
n])
 
  189            break;
  190        }
  191        if (
n == 4 && !tet_found) {
 
  192          tet = it_tet;
  193          tet_found = true;
  195          THROW_MESSAGE(
"More that one tet with the same connectivity");
 
  196        }
  197      }
  198    }
  199    if (!tet_found) {
  200      
  201      CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
 
  203                                             &tet, 1, &parent);
  205    }
  206    return tet;
  207  };
  208 
  209  
  211 
  213  for (auto m_tet : mother_tets) {
  216    
  217    int nb_mother_verts = 0;
  218    for (int nn = 0; nn != 4; ++nn) {
  219      if (conn[nn] == father) {
  221                "Tet has father vertex, impossible but here it is");
  222      }
  223      if (conn[nn] == mother) {
  224        new_conn[nn] = father;
  225        nb_mother_verts++;
  226      } else {
  227        new_conn[nn] = conn[nn];
  228      }
  229    }
  230    if (nb_mother_verts != 1) {
  232               "Tet should have only one vertex but have %d", nb_mother_verts);
  233    }
  234 
  236 
  237    
  238    created_tets.insert(create_tet(new_conn, m_tet));
  239  }
  240 
  241  
  244    for (
int dd = 1; 
dd <= 2; 
dd++)
 
  245      CHKERR m_field.get_moab().get_adjacencies(ents, dd, create, adj,
 
  246                                                moab::Interface::UNION);
  247    return adj;
  248  };
  249  auto adj_crated_ents = 
get_adj_ents(created_tets, 
true);
 
  250  adj_crated_ents.erase(common_edge[0]);
  251 
  253  for (auto ent : adj_crated_ents) {
  254    int num_nodes;
  257    int ii = 0;
  258    bool father_node = false;
  259    for (int nn = 0; nn != num_nodes; nn++) {
  260      if (conn[nn] == father)
  261        father_node = true;
  262      else
  263        small_conn[ii++] = conn[nn];
  264    }
  265    if (father_node) {
  266      if (ii > 1) 
  267        std::sort(&small_conn[0], &small_conn[ii]);
  268      if (ii == 2) 
  269        face_map.insert(FaceMap(ent, small_conn[0], small_conn[1]));
  270       else 
  271        face_map.insert(FaceMap(ent, small_conn[0], 0));
  272    }
  273  }
  274 
  275  auto adj_mother_ents = 
get_adj_ents(mother_tets, 
false);
 
  276  adj_mother_ents.erase(common_edge[0]);
  277  for (auto ent : adj_mother_ents) {
  278    int num_nodes;
  281    int nb_new_node = 0;
  282    int nn = 0;
  283    int ii = 0;
  284    for (; nn != num_nodes; ++nn) {
  285      if (conn[nn] == mother) {
  286        nb_new_node++;
  287      } else {
  288        small_conn[ii++] = conn[nn];
  289      }
  290    }
  291    if (nb_new_node > 0) {
  292      if (ii > 1) 
  293        std::sort(&small_conn[0], &small_conn[ii]);
  294      
  296      if (ii == 2) 
  297        n1 = small_conn[1];
  298      
  299      FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
  300      if (fit != face_map.end()) {
  303        if (m_field.get_moab().dimension_from_handle(parent) !=
  304            m_field.get_moab().dimension_from_handle(child))
  306                  "Huston we have a problem!");
  307 
  308        
  310                                               &child, 1, &parent);
  311        
  313      }
  314    }
  315  }
  316 
  317  
  319  if (tets_ptr != NULL) {
  320    seed_tets.merge(*tets_ptr);
  321    seed_tets = subtract(seed_tets, edge_tets);
  322    seed_tets = subtract(seed_tets, mother_tets);
  323  }
  324  seed_tets.merge(created_tets);
  325  out_tets.swap(seed_tets);
  326 
  328 
  330    std::cerr << "Nodes merged" << endl;
  331 
  333}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
const double n
refractive index of diffusive medium
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
UBlasVector< double > VectorDouble
Tag get_th_RefParentHandle() const
multi_index_container< FaceMap, indexed_by< hashed_unique< composite_key< FaceMap, member< FaceMap, EntityHandle, &FaceMap::n0 >, member< FaceMap, EntityHandle, &FaceMap::n1 > > > > > FaceMapIdx
MoFEMErrorCode lineSearch(Range &check_tests, EntityHandle father, EntityHandle mother, int line_search, FTensor::Tensor1< double, 3 > &t_move, Tag th=NULL)
Use bisection method to find point of edge collapse.