Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "EFAFragment3D.h"
      11             : 
      12             : #include "EFAVolumeNode.h"
      13             : #include "EFANode.h"
      14             : #include "EFAEdge.h"
      15             : #include "EFAFace.h"
      16             : #include "EFAFuncs.h"
      17             : #include "EFAElement3D.h"
      18             : #include "EFAError.h"
      19             : 
      20       37441 : EFAFragment3D::EFAFragment3D(EFAElement3D * host,
      21             :                              bool create_faces,
      22             :                              const EFAElement3D * from_host,
      23       37441 :                              unsigned int frag_id)
      24       37441 :   : EFAFragment(), _host_elem(host)
      25             : {
      26       37441 :   if (create_faces)
      27             :   {
      28       34572 :     if (!from_host)
      29           0 :       EFAError("EFAfragment3D constructor must have a from_host to copy from");
      30       34572 :     if (frag_id == std::numeric_limits<unsigned int>::max()) // copy the from_host itself
      31             :     {
      32      139007 :       for (unsigned int i = 0; i < from_host->numFaces(); ++i)
      33      115686 :         _faces.push_back(new EFAFace(*from_host->getFace(i)));
      34             :     }
      35             :     else
      36             :     {
      37       11251 :       if (frag_id > from_host->numFragments() - 1)
      38           0 :         EFAError("In EFAfragment3D constructor fragment_copy_index out of bounds");
      39       71954 :       for (unsigned int i = 0; i < from_host->getFragment(frag_id)->numFaces(); ++i)
      40       60703 :         _faces.push_back(new EFAFace(*from_host->getFragmentFace(frag_id, i)));
      41             :     }
      42       34572 :     findFacesAdjacentToFaces(); // IMPORTANT
      43             :   }
      44       37441 : }
      45             : 
      46       74882 : EFAFragment3D::~EFAFragment3D()
      47             : {
      48      229271 :   for (unsigned int i = 0; i < _faces.size(); ++i)
      49             :   {
      50      191830 :     if (_faces[i])
      51             :     {
      52      191830 :       delete _faces[i];
      53      191830 :       _faces[i] = nullptr;
      54             :     }
      55             :   }
      56       74882 : }
      57             : 
      58             : void
      59      543415 : EFAFragment3D::switchNode(EFANode * new_node, EFANode * old_node)
      60             : {
      61     3199586 :   for (unsigned int i = 0; i < _faces.size(); ++i)
      62     2656171 :     _faces[i]->switchNode(new_node, old_node);
      63      543415 : }
      64             : 
      65             : bool
      66     2086584 : EFAFragment3D::containsNode(EFANode * node) const
      67             : {
      68             :   bool contains = false;
      69     6876053 :   for (unsigned int i = 0; i < _faces.size(); ++i)
      70             :   {
      71     6047391 :     if (_faces[i]->containsNode(node))
      72             :     {
      73             :       contains = true;
      74             :       break;
      75             :     }
      76             :   }
      77     2086584 :   return contains;
      78             : }
      79             : 
      80             : unsigned int
      81       26673 : EFAFragment3D::getNumCuts() const
      82             : {
      83             :   unsigned int num_cut_faces = 0;
      84      160211 :   for (unsigned int i = 0; i < _faces.size(); ++i)
      85      133538 :     if (_faces[i]->hasIntersection())
      86        4843 :       num_cut_faces += 1;
      87       26673 :   return num_cut_faces;
      88             : }
      89             : 
      90             : unsigned int
      91           0 : EFAFragment3D::getNumCutNodes() const
      92             : {
      93           0 :   mooseError("Not implemented yet for 3D.");
      94             : }
      95             : 
      96             : std::set<EFANode *>
      97       55932 : EFAFragment3D::getAllNodes() const
      98             : {
      99             :   std::set<EFANode *> all_nodes;
     100      339065 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     101     1310731 :     for (unsigned int j = 0; j < _faces[i]->numNodes(); ++j)
     102     1027598 :       all_nodes.insert(_faces[i]->getNode(j));
     103       55932 :   return all_nodes;
     104             : }
     105             : 
     106             : bool
     107       26277 : EFAFragment3D::isConnected(EFAFragment * other_fragment) const
     108             : {
     109       26277 :   EFAFragment3D * other_frag3d = dynamic_cast<EFAFragment3D *>(other_fragment);
     110       26277 :   if (!other_frag3d)
     111           0 :     EFAError("in isConnected other_fragment is not of type EFAfragement3D");
     112             : 
     113       96641 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     114      506800 :     for (unsigned int j = 0; j < other_frag3d->numFaces(); ++j)
     115      436436 :       if (_faces[i]->equivalent(other_frag3d->_faces[j]))
     116             :         return true;
     117             : 
     118             :   return false;
     119             : }
     120             : 
     121             : bool
     122       30676 : EFAFragment3D::isEdgeConnected(EFAFragment * other_fragment) const
     123             : {
     124       30676 :   EFAFragment3D * other_frag3d = dynamic_cast<EFAFragment3D *>(other_fragment);
     125       30676 :   if (!other_frag3d)
     126           0 :     EFAError("in isEdgeConnected other_fragment is not of type EFAfragement3D");
     127             : 
     128       69267 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     129      238872 :     for (unsigned int j = 0; j < _faces[i]->numEdges(); ++j)
     130     1080882 :       for (unsigned int k = 0; k < other_frag3d->numFaces(); ++k)
     131     4096716 :         for (unsigned int l = 0; l < other_frag3d->_faces[k]->numEdges(); ++l)
     132     3216115 :           if (_faces[i]->getEdge(j)->equivalent(*(other_frag3d->_faces[k]->getEdge(l))))
     133             :             return true;
     134             : 
     135             :   return false;
     136             : }
     137             : 
     138             : void
     139        3352 : EFAFragment3D::removeInvalidEmbeddedNodes(std::map<unsigned int, EFANode *> & EmbeddedNodes)
     140             : {
     141             :   // N.B. this method is only called before we update fragments
     142             :   // N.B. an embedded node is valid IF at least one of its host faces is exterior and has more than
     143             :   // 1 cuts
     144             :   // TODO: the invalid cases are generalized from 2D. The method may need improvements in 3D
     145        3352 :   if (hasFaceWithOneCut())
     146             :   {
     147             :     // build a local inverse map for all emb cut nodes in this fragment
     148             :     std::map<EFANode *, std::vector<EFAFace *>> emb_inverse_map;
     149        1394 :     for (unsigned int i = 0; i < _faces.size(); ++i)
     150        5540 :       for (unsigned int j = 0; j < _faces[i]->numEdges(); ++j)
     151        4368 :         if (_faces[i]->getEdge(j)->hasIntersection())
     152             :         {
     153         908 :           EFANode * emb_node = _faces[i]->getEdge(j)->getEmbeddedNode(0);
     154         908 :           emb_inverse_map[emb_node].push_back(_faces[i]);
     155             :         }
     156             : 
     157             :     // find all invalid embedded nodes
     158             :     std::vector<EFANode *> invalid_emb;
     159             :     std::map<EFANode *, std::vector<EFAFace *>>::iterator it;
     160         676 :     for (it = emb_inverse_map.begin(); it != emb_inverse_map.end(); ++it)
     161             :     {
     162         454 :       EFANode * emb_node = it->first;
     163             :       std::vector<EFAFace *> & emb_faces = it->second;
     164         454 :       if (emb_faces.size() != 2)
     165           0 :         EFAError("one embedded node must be owned by 2 faces");
     166             :       unsigned int counter = 0;
     167        1362 :       for (unsigned int i = 0; i < emb_faces.size(); ++i)
     168             :       {
     169         908 :         unsigned int face_id = getFaceID(emb_faces[i]);
     170         908 :         if (!isFaceInterior(face_id) && emb_faces[i]->hasIntersection())
     171         464 :           counter += 1; // count the appropriate emb's faces
     172             :       }
     173         454 :       if (counter == 0)
     174           0 :         invalid_emb.push_back(emb_node);
     175             :     }
     176             : 
     177             :     // delete all invalid emb nodes
     178         222 :     for (unsigned int i = 0; i < invalid_emb.size(); ++i)
     179             :     {
     180           0 :       Efa::deleteFromMap(EmbeddedNodes, invalid_emb[i]);
     181           0 :       _host_elem->removeEmbeddedNode(invalid_emb[i], true); // also remove from neighbors
     182             :     }                                                       // i
     183         222 :   }
     184        3352 : }
     185             : 
     186             : void
     187         254 : EFAFragment3D::combine_tip_faces()
     188             : {
     189         254 :   if (!_host_elem)
     190           0 :     EFAError("In combine_tip_faces() the frag must have host_elem");
     191             : 
     192        1618 :   for (unsigned int i = 0; i < _host_elem->numFaces(); ++i)
     193             :   {
     194             :     std::vector<unsigned int> frag_tip_face_id;
     195        9793 :     for (unsigned int j = 0; j < _faces.size(); ++j)
     196        8429 :       if (_host_elem->getFace(i)->containsFace(_faces[j]))
     197        1628 :         frag_tip_face_id.push_back(j);
     198             : 
     199        1364 :     if (frag_tip_face_id.size() == 2) // combine the two frag faces on this elem face
     200         264 :       combine_two_faces(frag_tip_face_id[0], frag_tip_face_id[1], _host_elem->getFace(i));
     201        1364 :   }
     202             :   // TODO: may need to combine other frag faces that have tip edges
     203         254 : }
     204             : 
     205             : bool
     206     4791281 : EFAFragment3D::isFaceInterior(unsigned int face_id) const
     207             : {
     208     4791281 :   if (!_host_elem)
     209           0 :     EFAError("in isFaceInterior() fragment must have host elem");
     210             : 
     211             :   bool face_in_elem_face = false;
     212    17651995 :   for (unsigned int i = 0; i < _host_elem->numFaces(); ++i)
     213    16753613 :     if (_host_elem->getFace(i)->containsFace(_faces[face_id]))
     214             :     {
     215             :       face_in_elem_face = true;
     216             :       break;
     217             :     }
     218             : 
     219             :   // the face is interior if it does not coincide with an element face
     220     4791281 :   return !face_in_elem_face;
     221             : }
     222             : 
     223             : std::vector<unsigned int>
     224           0 : EFAFragment3D::get_interior_face_id() const
     225             : {
     226             :   std::vector<unsigned int> interior_face_id;
     227           0 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     228           0 :     if (isFaceInterior(i))
     229           0 :       interior_face_id.push_back(i);
     230             : 
     231           0 :   return interior_face_id;
     232           0 : }
     233             : 
     234             : bool
     235           0 : EFAFragment3D::isThirdInteriorFace(unsigned int face_id) const
     236             : {
     237           0 :   if (!_host_elem)
     238           0 :     EFAError("in isThirdInteriorFace fragment must have host elem");
     239             : 
     240           0 :   for (unsigned int i = 0; i < _host_elem->numInteriorNodes(); ++i)
     241           0 :     if (_faces[face_id]->containsNode(_host_elem->getInteriorNode(i)->getNode()))
     242             :       return true;
     243             : 
     244             :   return false;
     245             : }
     246             : 
     247             : unsigned int
     248     7770563 : EFAFragment3D::numFaces() const
     249             : {
     250     7770563 :   return _faces.size();
     251             : }
     252             : 
     253             : EFAFace *
     254     1309199 : EFAFragment3D::getFace(unsigned int face_id) const
     255             : {
     256     1309199 :   if (face_id > _faces.size() - 1)
     257           0 :     EFAError("in EFAfragment3D::get_face, index out of bounds");
     258     1309199 :   return _faces[face_id];
     259             : }
     260             : 
     261             : unsigned int
     262         908 : EFAFragment3D::getFaceID(EFAFace * face) const
     263             : {
     264        2942 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     265        2942 :     if (_faces[i] == face)
     266         908 :       return i;
     267           0 :   EFAError("face not found in get_face_id()");
     268             : }
     269             : 
     270             : void
     271       15705 : EFAFragment3D::addFace(EFAFace * new_face)
     272             : {
     273       15705 :   _faces.push_back(new_face);
     274       15705 : }
     275             : 
     276             : std::set<EFANode *>
     277           0 : EFAFragment3D::getFaceNodes(unsigned int face_id) const
     278             : {
     279             :   std::set<EFANode *> face_nodes;
     280           0 :   for (unsigned int i = 0; i < _faces[face_id]->numNodes(); ++i)
     281           0 :     face_nodes.insert(_faces[face_id]->getNode(i));
     282           0 :   return face_nodes;
     283             : }
     284             : 
     285             : EFAElement3D *
     286           0 : EFAFragment3D::getHostElement() const
     287             : {
     288           0 :   return _host_elem;
     289             : }
     290             : 
     291             : std::vector<EFAFragment3D *>
     292        1657 : EFAFragment3D::split()
     293             : {
     294             :   // This method will split one existing fragment into one or two new fragments
     295             :   std::vector<EFAFragment3D *> new_fragments;
     296             :   std::vector<std::vector<EFAFace *>> all_subfaces;
     297       10095 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     298             :   {
     299        8438 :     std::vector<EFAFace *> subfaces = _faces[i]->split();
     300        8438 :     all_subfaces.push_back(subfaces);
     301        8438 :   }
     302             : 
     303             :   // build new frags
     304        1657 :   if (hasFaceWithOneCut()) // "fakely" cut fragment
     305             :   {
     306         445 :     EFAFragment3D * new_frag = new EFAFragment3D(_host_elem, false, nullptr);
     307        2827 :     for (unsigned int i = 0; i < all_subfaces.size(); ++i)
     308        5231 :       for (unsigned int j = 0; j < all_subfaces[i].size(); ++j)
     309        2849 :         new_frag->addFace(all_subfaces[i][j]);
     310         445 :     new_frag->findFacesAdjacentToFaces();
     311         445 :     new_fragments.push_back(new_frag);
     312             :   }
     313             :   else // thoroughly cut fragment
     314             :   {
     315             :     // find the first face with 2 sub-faces
     316             :     EFAFace * start_face1 = nullptr;
     317             :     EFAFace * start_face2 = nullptr;
     318             :     unsigned int startOldFaceID = 0;
     319        1385 :     for (unsigned int i = 0; i < _faces.size(); ++i)
     320             :     {
     321        1385 :       if (all_subfaces[i].size() == 2)
     322             :       {
     323        1212 :         start_face1 = all_subfaces[i][0];
     324        1212 :         start_face2 = all_subfaces[i][1];
     325             :         startOldFaceID = i;
     326        1212 :         break;
     327             :       }
     328             :     } // i
     329        1212 :     EFAFragment3D * new_frag1 = connectSubfaces(start_face1, startOldFaceID, all_subfaces);
     330        1212 :     EFAFragment3D * new_frag2 = connectSubfaces(start_face2, startOldFaceID, all_subfaces);
     331        1212 :     new_fragments.push_back(new_frag1);
     332        1212 :     new_fragments.push_back(new_frag2);
     333             :   }
     334        1657 :   return new_fragments;
     335        1657 : }
     336             : 
     337             : void
     338       40129 : EFAFragment3D::findFacesAdjacentToFaces()
     339             : {
     340       40129 :   _faces_adjacent_to_faces.clear();
     341      244089 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     342             :   {
     343      203960 :     std::vector<EFAFace *> face_adjacents(_faces[i]->numEdges(), nullptr);
     344     1281010 :     for (unsigned int j = 0; j < _faces.size(); ++j)
     345     1077050 :       if (_faces[j] != _faces[i] && _faces[i]->isAdjacent(_faces[j]))
     346             :       {
     347      739092 :         unsigned int adj_edge = _faces[i]->adjacentCommonEdge(_faces[j]);
     348      739092 :         face_adjacents[adj_edge] = _faces[j];
     349             :       }
     350             : 
     351      203960 :     _faces_adjacent_to_faces.push_back(face_adjacents);
     352      203960 :   } // i
     353       40129 : }
     354             : 
     355             : EFAFace *
     356           0 : EFAFragment3D::getAdjacentFace(unsigned int face_id, unsigned int edge_id) const
     357             : {
     358           0 :   return _faces_adjacent_to_faces[face_id][edge_id]; // possibly NULL
     359             : }
     360             : 
     361             : void
     362           0 : EFAFragment3D::removeEmbeddedNode(EFANode * emb_node)
     363             : {
     364           0 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     365           0 :     _faces[i]->removeEmbeddedNode(emb_node);
     366           0 : }
     367             : 
     368             : bool
     369        5009 : EFAFragment3D::hasFaceWithOneCut() const
     370             : {
     371             :   // N.B. this method can only be implemented when the fragment has just been marked
     372       28158 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     373       23816 :     if (_faces[i]->getNumCuts() == 1)
     374             :       return true;
     375             :   return false;
     376             : }
     377             : 
     378             : void
     379       54780 : EFAFragment3D::getNodeInfo(std::vector<std::vector<unsigned int>> & face_node_indices,
     380             :                            std::vector<EFANode *> & nodes) const
     381             : {
     382             :   // get all nodes' pointers - a vector
     383       54780 :   std::set<EFANode *> all_node_set = getAllNodes();
     384       54780 :   nodes.resize(all_node_set.size());
     385       54780 :   std::copy(all_node_set.begin(), all_node_set.end(), nodes.begin());
     386             : 
     387             :   // get face connectivity
     388       54780 :   face_node_indices.clear();
     389      331001 :   for (unsigned int i = 0; i < _faces.size(); ++i)
     390             :   {
     391             :     std::vector<unsigned int> line_face_indices;
     392     1276171 :     for (unsigned int j = 0; j < _faces[i]->numNodes(); ++j)
     393             :     {
     394      999950 :       EFANode * node = _faces[i]->getNode(j);
     395      999950 :       unsigned int vec_index = std::find(nodes.begin(), nodes.end(), node) - nodes.begin();
     396      999950 :       line_face_indices.push_back(vec_index);
     397             :     }
     398      276221 :     face_node_indices.push_back(line_face_indices);
     399      276221 :   }
     400       54780 : }
     401             : 
     402             : EFAFragment3D *
     403        2424 : EFAFragment3D::connectSubfaces(EFAFace * start_face,
     404             :                                unsigned int startOldFaceID,
     405             :                                std::vector<std::vector<EFAFace *>> & subfaces)
     406             : {
     407             :   // this method is only called in EFAfragment3D::split()
     408        2424 :   std::vector<bool> contributed(subfaces.size(), false);
     409             :   contributed[startOldFaceID] = true;
     410             :   unsigned int num_contrib_faces = 1;
     411             :   unsigned int old_num_contrib = 1;
     412        2424 :   std::vector<EFAFace *> frag_faces(1, start_face);
     413             : 
     414             :   // collect all subfaces connected to start_face
     415             :   do
     416             :   {
     417             :     old_num_contrib = num_contrib_faces;
     418       29156 :     for (unsigned int i = 0; i < subfaces.size(); ++i)
     419             :     {
     420       24296 :       if (!contributed[i]) // old face not contributed to new fragment yet
     421             :       {
     422             :         bool adjacent_found = false;
     423       17928 :         for (unsigned int j = 0; j < subfaces[i].size(); ++j)
     424             :         {
     425       33477 :           for (unsigned int k = 0; k < frag_faces.size(); ++k)
     426             :           {
     427       26929 :             if (subfaces[i][j]->isAdjacent(frag_faces[k]))
     428             :             {
     429             :               adjacent_found = true;
     430             :               contributed[i] = true;
     431        8008 :               frag_faces.push_back(subfaces[i][j]);
     432        8008 :               num_contrib_faces += 1;
     433             :               break;
     434             :             }
     435             :           } // k
     436             :           if (adjacent_found)
     437             :             break;
     438             :         } // j
     439             :       }
     440             :     } // i, loop over all old faces
     441        4860 :   } while (num_contrib_faces != old_num_contrib);
     442             : 
     443             :   // get the cut plane face
     444             :   std::vector<EFAEdge *> cut_plane_edges;
     445        2424 :   EFAFragment3D * new_frag = new EFAFragment3D(_host_elem, false, nullptr);
     446       12856 :   for (unsigned int i = 0; i < frag_faces.size(); ++i)
     447       10432 :     new_frag->addFace(frag_faces[i]);
     448        2424 :   new_frag->findFacesAdjacentToFaces();
     449             : 
     450       12856 :   for (unsigned int i = 0; i < new_frag->numFaces(); ++i)
     451             :   {
     452       10432 :     EFAEdge * lone_edge = new_frag->loneEdgeOnFace(i);
     453       10432 :     if (lone_edge != nullptr) // valid edge
     454        8752 :       cut_plane_edges.push_back(new EFAEdge(*lone_edge));
     455             :   }
     456             : 
     457        2424 :   EFAFace * cut_face = new EFAFace(cut_plane_edges.size());
     458       11176 :   for (unsigned int i = 0; i < cut_plane_edges.size(); ++i)
     459        8752 :     cut_face->setEdge(i, cut_plane_edges[i]);
     460        2424 :   cut_face->sortEdges();
     461        2424 :   cut_face->reverseEdges();
     462        2424 :   cut_face->createNodes();
     463             : 
     464             :   // finalize the new fragment
     465        2424 :   new_frag->addFace(cut_face);
     466        2424 :   new_frag->findFacesAdjacentToFaces();
     467        2424 :   return new_frag;
     468        2424 : }
     469             : 
     470             : EFAEdge *
     471       10432 : EFAFragment3D::loneEdgeOnFace(unsigned int face_id) const
     472             : {
     473             :   // if any face edge is not shared by any other face, we call it a lone edge
     474       34183 :   for (unsigned int i = 0; i < _faces[face_id]->numEdges(); ++i)
     475       32503 :     if (_faces_adjacent_to_faces[face_id][i] == nullptr)
     476        8752 :       return _faces[face_id]->getEdge(i);
     477             :   return nullptr;
     478             : }
     479             : 
     480             : void
     481         264 : EFAFragment3D::combine_two_faces(unsigned int face_id1,
     482             :                                  unsigned int face_id2,
     483             :                                  const EFAFace * elem_face)
     484             : {
     485             :   // get the new full face
     486         264 :   EFAFace * full_face = _faces[face_id1]->combineWithFace(_faces[face_id2]);
     487         264 :   full_face->resetEdgeIntersection(elem_face); // IMPORTANT
     488             : 
     489             :   // take care of the common adjacent faces (combine their tip edges)
     490             :   std::set<EFAFace *> face1_neigh;
     491         264 :   face1_neigh.insert(_faces_adjacent_to_faces[face_id1].begin(),
     492             :                      _faces_adjacent_to_faces[face_id1].end());
     493             :   std::set<EFAFace *> face2_neigh;
     494         264 :   face2_neigh.insert(_faces_adjacent_to_faces[face_id2].begin(),
     495             :                      _faces_adjacent_to_faces[face_id2].end());
     496         264 :   std::vector<EFAFace *> common_adjacent_faces = Efa::getCommonElems(face1_neigh, face2_neigh);
     497             : 
     498         782 :   for (unsigned int i = 0; i < common_adjacent_faces.size(); ++i)
     499             :   {
     500         518 :     EFAFace * comm_face = common_adjacent_faces[i];
     501         518 :     if (comm_face != nullptr)
     502             :     {
     503         508 :       unsigned int edge_id1 = comm_face->adjacentCommonEdge(_faces[face_id1]);
     504         508 :       unsigned int edge_id2 = comm_face->adjacentCommonEdge(_faces[face_id2]);
     505         508 :       comm_face->combineTwoEdges(edge_id1, edge_id2);
     506         508 :       comm_face->resetEdgeIntersection(elem_face); // IMPORTANT
     507             :     }
     508             :   }
     509             : 
     510             :   // delete old faces and update private members of EFAfragment3D
     511         264 :   delete _faces[face_id1];
     512         264 :   delete _faces[face_id2];
     513         264 :   _faces[face_id1] = full_face;
     514         264 :   _faces.erase(_faces.begin() + face_id2);
     515         264 :   findFacesAdjacentToFaces(); // rebuild _adjacent_face_ix: IMPORTANT
     516         528 : }