36    const bool only_if_improve_quality, 
const double move,
 
   37    const int line_search, 
Tag th, 
int verb) {
 
   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);
 
   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);
 
   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);
 
   74            "no common edge between nodes");
 
   75  } 
else if (common_edge.empty()) {
 
   77    if (tets_ptr != NULL) {
 
   78      seed_tets.merge(*tets_ptr);
 
   87  CHKERR m_field.
get_moab().get_adjacencies(common_edge, 3, 
true, edge_tets);
 
   89  if (tets_ptr != NULL) {
 
   90    edge_tets = intersect(edge_tets, *tets_ptr);
 
   93  mother_tets = subtract(mother_tets, edge_tets);
 
   94  father_tets = subtract(father_tets, edge_tets);
 
   97                               const int num_nodes) {
 
  100      CHKERR m_field.
get_moab().get_coords(conn, num_nodes, &coords[0]);
 
  107  auto get_tensor = [](
VectorDouble &coords, 
const int shift) {
 
  109        &coords[shift], &coords[shift + 1], &coords[shift + 2]);
 
  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));
 
  122  if (line_search > 0) {
 
  123    Range check_tests = unite(father_tets, mother_tets);
 
  127  if (only_if_improve_quality) {
 
  128    Range check_tests = mother_tets;
 
  130    if (move > 0 || line_search) {
 
  131      check_tests.merge(father_tets);
 
  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;
 
  153  if (move > 0 || line_search) {
 
  162                             int *ret_num_nodes = 
nullptr) {
 
  165    CHKERR m_field.
get_moab().get_connectivity(ent, conn, num_nodes, 
true);
 
  167      *ret_num_nodes = num_nodes;
 
  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) {
 
  180      CHKERR m_field.
get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
 
  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])
 
  191        if (
n == 4 && !tet_found) {
 
  195          THROW_MESSAGE(
"More that one tet with the same connectivity");
 
  201      CHKERR m_field.
get_moab().create_element(MBTET, new_conn, 4, tet);
 
  213  for (
auto m_tet : mother_tets) {
 
  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");
 
  223      if (conn[nn] == mother) {
 
  224        new_conn[nn] = father;
 
  227        new_conn[nn] = conn[nn];
 
  230    if (nb_mother_verts != 1) {
 
  232               "Tet should have only one vertex but have %d", nb_mother_verts);
 
  238    created_tets.insert(create_tet(new_conn, m_tet));
 
  242  auto get_adj_ents = [&](
const Range &ents, 
const bool create) {
 
  244    for (
int dd = 1; dd <= 2; dd++)
 
  246                                                moab::Interface::UNION);
 
  249  auto adj_crated_ents = get_adj_ents(created_tets, 
true);
 
  250  adj_crated_ents.erase(common_edge[0]);
 
  253  for (
auto ent : adj_crated_ents) {
 
  258    bool father_node = 
false;
 
  259    for (
int nn = 0; nn != num_nodes; nn++) {
 
  260      if (conn[nn] == father)
 
  263        small_conn[ii++] = conn[nn];
 
  267        std::sort(&small_conn[0], &small_conn[ii]);
 
  269        face_map.insert(
FaceMap(ent, small_conn[0], small_conn[1]));
 
  271        face_map.insert(
FaceMap(ent, small_conn[0], 0));
 
  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) {
 
  284    for (; nn != num_nodes; ++nn) {
 
  285      if (conn[nn] == mother) {
 
  288        small_conn[ii++] = conn[nn];
 
  291    if (nb_new_node > 0) {
 
  293        std::sort(&small_conn[0], &small_conn[ii]);
 
  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!");
 
  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);
 
  324  seed_tets.merge(created_tets);
 
  325  out_tets.swap(seed_tets);
 
  330    std::cerr << 
"Nodes merged" << endl;