src/MEDLoader/CrackAlgo.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email :
  18 // webmaster.salome@opencascade.com
  19 //
  20 // Author : Aymeric SONOLET (CEA/DES)
  21
  22 #include "CrackAlgo.hxx"
  23
  24 #include <algorithm>
  25 #include <unordered_map>
  26 #include <unordered_set>
  27 #include <string>
  28 #include <memory>
  29 #include <utility>
  30
  31 #include "./MEDFileMesh.hxx"
  32
  33 #include <MEDCouplingUMesh.hxx>
  34 #include <MEDCouplingMemArray.hxx>
  35 #include <InterpKernelException.hxx>
  36 #include <MCType.hxx>
  37 #include <MCAuto.hxx>
  38
  39 using namespace MEDCoupling;
  40 using namespace std;
  41 using MCU = MCAuto<MEDCouplingUMesh>;
  42 using DAI = MCAuto<DataArrayIdType>;
  43 using DAD = MCAuto<DataArrayDouble>;
  44
  45
  46 CrackAlgo::Map2Map
  47 CrackAlgo::Compute(
  48     MEDFileUMesh* mm,
  49     const string& grp_name,
  50     bool grpMustBeFullyDup
  51 ) {
  52
  53     vector<int> levs = mm->getNonEmptyLevels();
  54     if (find(levs.begin(), levs.end(), 0) == levs.end() ||
  55             find(levs.begin(), levs.end(), -1) == levs.end())
  56         throw INTERP_KERNEL::Exception(
  57                 "MEDFileUMesh::duplicateFaces : This method "
  58                 "works only for mesh defined on level 0 and -1 !");
  59
  60     MCU m0 = mm->getMeshAtLevel(0);
  61
  62     MCU crackingMesh = mm->getGroup(-1, grp_name);
  63     MCU cleanCrackingMesh = CleanM1Mesh(*m0, *crackingMesh);
  64
  65     DAI c2f(DataArrayIdType::New()),
  66         c2fIdx(DataArrayIdType::New()),
  67         f2c(DataArrayIdType::New()),
  68         f2cIdx(DataArrayIdType::New());
  69     const MCU mf = m0->buildDescendingConnectivity(c2f, c2fIdx, f2c, f2cIdx);
  70
  71     DataArrayIdType * crackToMfP;
  72     const bool crackCellsInMf = mf->areCellsIncludedIn(cleanCrackingMesh, 2, crackToMfP);
  73     DAI f2dupIdInMf(crackToMfP);
  74     if (!crackCellsInMf)
  75         throw INTERP_KERNEL::Exception(
  76             "crackAlong: All cells in crack are not part of Mf.");
  77
  78     DAI n2c(DataArrayIdType::New()),
  79           n2cIdx(DataArrayIdType::New());
  80     m0->getReverseNodalConnectivity(n2c, n2cIdx);
  81
  82     const Map2Set n2c_dup = GetNode2CellMap(mf, n2cIdx, n2c, f2dupIdInMf);
  83
  84     const Set cTouchingN_dup = GetCellsTouchingNodesToDup(
  85         mf, n2cIdx, n2c, f2dupIdInMf);
  86
  87     const Graph c2c = BuildCutC2CGraph(
  88         c2fIdx, c2f, f2cIdx, f2c, cTouchingN_dup, f2dupIdInMf);
  89
  90     Map2Map cellOld2NewNode = CreateNewNodesInTopLevelMesh(n2c_dup, c2c, m0);
  91     // End of node creation, separation of M0
  92
  93     // Faces in M1 must be added/updated
  94
  95     MCU m1 = mm->getMeshAtLevel(-1);
  96
  97     // Before changing M1
  98     const DAI f2dupIdInM1 = GetFacesToDupInM1(*cleanCrackingMesh, *m1);
  99
 100     const auto f2changeM1Mf = GetFacesInM1TouchingDuplicatedNodes(
 101         n2c_dup, f2dupIdInM1, *mf, *m1);
 102     const DAI f2changeIdInM1 = f2changeM1Mf.first;
 103     const DAI f2changeIdInMf = f2changeM1Mf.second;
 104
 105     AddMissingElementsOnLevelM1AndChangeConnectivity(
 106         f2cIdx, f2c, f2dupIdInM1, f2dupIdInMf, cellOld2NewNode, m1, grpMustBeFullyDup);
 107
 108     ChangeConnectivityOfM1Elements(
 109             f2changeIdInM1, f2changeIdInMf, cellOld2NewNode, f2cIdx, f2c, m1);
 110
 111     // TODO(aymeric): If one wants to implement
 112     // AddMissingElementsOnLevelM2AndChangeConnectivity and
 113     // ChangeConnectivityOfM2Elements it should be there.
 114
 115     DAI famLev0 = CopyFamilyArrAtLev0(mm);
 116     DAI newFam = CopyAndCompleteFamilyArrAtLevelM1(mm, m1, f2dupIdInM1);
 117
 118     mm->setMeshAtLevel(0, m0);
 119     mm->setMeshAtLevel(-1, m1);
 120     // mm->setMeshAtLevel(-2, m2);
 121
 122     mm->setFamilyFieldArr(0, famLev0);
 123     mm->setFamilyFieldArr(-1, newFam);
 124     // mm->setFamilyFieldArr(-2, newFam2);
 125
 126     const Map2Set addedNodes = BuildMap2Set(cellOld2NewNode);
 127     CompleteFamilyArrAtNodeLevel(addedNodes, mm);
 128
 129     return cellOld2NewNode;
 130 }
 131
 132 CrackAlgo::Map2Set
 133 CrackAlgo::BuildMap2Set(
 134     const Map2Map& cellOld2NewNode
 135 ) {
 136     Map2Set res;
 137     for (const auto & keyElem : cellOld2NewNode) {
 138         const auto & old2NewMap = keyElem.second;
 139         for (const auto & old2New : old2NewMap) {
 140             res[old2New.first].insert(old2New.second);
 141         }
 142     }
 143     return res;
 144 }
 145
 146 DataArrayIdType *
 147 CrackAlgo::CopyAndCompleteFamilyArrAtLevelM1(
 148     const MEDFileUMesh * mm,
 149     const MEDCouplingUMesh * mf,
 150     const DataArrayIdType * f2dup
 151 ) {
 152     DataArrayIdType* newFam = nullptr;
 153     const DataArrayIdType *fam = mm->getFamilyFieldAtLevel(-1);
 154     if (fam != nullptr) {
 155         newFam = DataArrayIdType::New();
 156         const auto fam_size = static_cast<size_t>(fam->getNumberOfTuples());
 157         newFam->alloc(static_cast<size_t>(mf->getNumberOfCells()));
 158         copy(fam->begin(), fam->end(), newFam->getPointer());
 159         mcIdType * fam_ptr = newFam->getPointer();
 160         size_t i_elem_added = 0;
 161         for (const auto & face : *f2dup) {
 162             // NOTE(aymeric): assumes all faces are duplicated
 163             fam_ptr[fam_size + i_elem_added] = fam_ptr[face];
 164             i_elem_added++;
 165         }
 166     }
 167     return newFam;
 168 }
 169
 170 DataArrayIdType *
 171 CrackAlgo::CopyFamilyArrAtLev0(
 172     const MEDFileUMesh * mm
 173 ) {
 174     DataArrayIdType* famLev0 = nullptr;
 175     const DataArrayIdType *fam = mm->getFamilyFieldAtLevel(0);
 176     if (fam != nullptr) {
 177         famLev0 = DataArrayIdType::New();
 178         const auto fam_size = static_cast<size_t>(fam->getNumberOfTuples());
 179         famLev0->alloc(fam_size);
 180         copy(fam->begin(), fam->end(), famLev0->getPointer());
 181     }
 182     return famLev0;
 183 }
 184
 185 void
 186 CrackAlgo::CompleteFamilyArrAtNodeLevel(
 187     const Map2Set& addedNodes,
 188     MEDFileUMesh * mm
 189 ) {
 190     DataArrayIdType *node_fam = mm->getFamilyFieldAtLevel(1);
 191     if (node_fam != nullptr) {
 192         node_fam->reAlloc(static_cast<size_t>(mm->getCoords()->getNumberOfTuples()));
 193         mcIdType * node_fam_ptr = node_fam->getPointer();
 194         for (const auto & old2NewPair : addedNodes) {
 195             const auto & oldNode = old2NewPair.first;
 196             const auto & newNodes = old2NewPair.second;
 197             for (const auto & newNode : newNodes) {
 198                 node_fam_ptr[newNode] = node_fam_ptr[oldNode];
 199             }
 200         }
 201     }
 202 }
 203
 204 MEDCouplingUMesh *
 205 CrackAlgo::CleanM1Mesh(
 206     const MEDCouplingUMesh & m,
 207     const MEDCouplingUMesh & crackingMesh
 208 ) {
 209     MCU m0skin = m.computeSkin();
 210     DataArrayIdType * idsToKeepP;
 211     m0skin->areCellsIncludedIn(&crackingMesh, 2, idsToKeepP);
 212     DAI idsToKeep(idsToKeepP);
 213     // discard cells on the skin of M0
 214     DAI ids2 = idsToKeep->findIdsNotInRange(0, m0skin->getNumberOfCells());
 215     MCU otherDimM1OnSameCoords =
 216         crackingMesh.buildPartOfMySelf(ids2->begin(), ids2->end(), true);
 217     return otherDimM1OnSameCoords.retn();
 218 }
 219
 220 void CrackAlgo::Dfs(
 221     const Graph& graph,
 222     const mcIdType &node,
 223     const size_t &componentId,
 224     map<mcIdType, bool>& visited,
 225     map<mcIdType, size_t>& componentMap
 226 ) {
 227     visited[node] = true;
 228     // Assign a unique component ID to the node
 229     componentMap[node] = componentId;
 230
 231     for (const auto& neighbor : graph.at(node)) {
 232         if (!visited[neighbor]) {
 233             Dfs(graph, neighbor, componentId, visited, componentMap);
 234         }
 235     }
 236 }
 237
 238 // Function to find connected components
 239 vector<shared_ptr<vector<mcIdType>>>
 240 CrackAlgo::FindConnectedComponents(
 241     const Graph& graph
 242 ) {
 243     map<mcIdType, bool> visited;
 244     map<mcIdType, size_t> componentMap;
 245     size_t componentId = 0;
 246
 247     for (const auto &nodePair : graph)
 248         visited[nodePair.first] = false;
 249
 250     for (const auto& nodePair : graph) {
 251         if (!visited[nodePair.first]) {
 252             Dfs(graph, nodePair.first, componentId, visited, componentMap);
 253             componentId++;
 254         }
 255     }
 256
 257     using spv = shared_ptr<vector<mcIdType>>;
 258     vector<spv> components;
 259     for (size_t i = 0; i < componentId; i++)
 260         components.push_back(make_shared<vector<mcIdType>>());
 261     for (const auto& nodePair : componentMap) {
 262         components[nodePair.second]->push_back(nodePair.first);
 263     }
 264
 265     return components;
 266 }
 267
 268 CrackAlgo::Set CrackAlgo::DataArrayToSet(
 269     const DataArrayIdType& da
 270 ) {
 271     Set res;
 272     for (const auto &elem : da) {
 273         res.insert(elem);
 274     }
 275     return res;
 276 }
 277
 278 DataArrayIdType * CrackAlgo::SetToDataArray(
 279     const Set& s
 280 ) {
 281     DataArrayIdType * da(DataArrayIdType::New());
 282     da->alloc(s.size());
 283     mcIdType * da_p = da->rwBegin();
 284     size_t i = 0;
 285     for (const auto & elem : s) {
 286         da_p[i] = elem;
 287         i++;
 288     }
 289     return da;
 290 }
 291
 292 CrackAlgo::Graph
 293 CrackAlgo::BuildCutC2CGraph(
 294     const DataArrayIdType * c2fIdx,
 295     const DataArrayIdType * c2f,
 296     const DataArrayIdType * f2cIdx,
 297     const DataArrayIdType * f2c,
 298     const Set & cTouchingN_dup,
 299     const DataArrayIdType * f2dup
 300 ) {
 301
 302     Set f2dup_set = DataArrayToSet(*f2dup);
 303     Graph c2c;
 304
 305     const mcIdType * c2fIdx_p {c2fIdx->begin()};
 306     const mcIdType * c2f_p {c2f->begin()};
 307     const mcIdType * f2cIdx_p {f2cIdx->begin()};
 308     const mcIdType * f2c_p {f2c->begin()};
 309
 310     for (const auto &cell : cTouchingN_dup) {
 311         c2c[cell] = std::unordered_set<mcIdType>();
 312         for (auto faceIdx = c2fIdx_p[cell]; faceIdx < c2fIdx_p[cell + 1]; faceIdx++) {
 313             const mcIdType &face = c2f_p[faceIdx];
 314             // face is not in face to duplicate
 315             if (f2dup_set.find(face) == f2dup_set.end()) {
 316                 for (mcIdType cellIdx {f2cIdx_p[face]}; cellIdx < f2cIdx_p[face +1]; cellIdx++) {
 317                     // getting the id of the other cell (on face is
 318                     // connected to two cells)
 319                     const auto& cell2 { f2c_p[cellIdx] };
 320                     if (
 321                         (cell != cell2)
 322                         & (cTouchingN_dup.find(cell2) != cTouchingN_dup.end())
 323                     ) {
 324                         c2c.at(cell).insert(cell2);
 325                     }
 326                 }
 327             }
 328         }
 329     }
 330     return c2c;
 331 }
 332
 333 CrackAlgo::Map2Map
 334 CrackAlgo::CreateNewNodesInTopLevelMesh(
 335     const Map2Set & n2c_dup,
 336     const Graph & c2c,
 337     MEDCouplingUMesh* m0
 338 ) {
 339     DataArrayIdType* c2n = m0->getNodalConnectivity();
 340     DataArrayIdType* c2nIdx = m0->getNodalConnectivityIndex();
 341     mcIdType * c2n_ptr = c2n->rwBegin();
 342     const mcIdType * c2nIdx_ptr = c2nIdx->begin();
 343
 344     DataArrayDouble* coords = m0->getCoords();
 345     mcIdType i = coords->getNumberOfTuples();
 346     const size_t coordsDim = coords->getNumberOfComponents();
 347
 348     Map2Map cellOld2NewNode;
 349
 350     for (const auto &nc_pair : n2c_dup) {
 351         const auto &node = nc_pair.first;
 352         const auto &cells = nc_pair.second;
 353
 354         // Building local connection graph
 355         Graph c2c_local;
 356         for (const auto &cell : cells) {
 357             auto& val = c2c_local[cell];
 358             if (c2c.find(cell) == c2c.end())
 359                 throw INTERP_KERNEL::Exception(
 360                     "A cell touching a node is not part of c2c.");
 361             for (const auto &neighbor : c2c.at(cell)) {
 362                 if (cells.find(neighbor) != cells.end())
 363                     val.insert(neighbor);
 364             }
 365         }
 366
 367         const vector<shared_ptr<vector<mcIdType>>> compo_connex = FindConnectedComponents(c2c_local);
 368
 369         // Duplicate node for compo[1:]
 370         const size_t i_nodes_to_add = compo_connex.size() - 1;
 371         coords->reAlloc(static_cast<size_t>(i) + i_nodes_to_add);
 372
 373         int i_compo = 0;
 374         for (const auto &compo : compo_connex) {
 375             if (i_compo > 0) {
 376                 // Coords are copied at the end of the vector of coords, ie the
 377                 // node is duplicated
 378                 coords->getTuple(
 379                         node,
 380                         coords->getPointer() + static_cast<size_t>(i) * coordsDim);
 381                 for (const auto &cell : *compo) {
 382                     // The node number is replaced in all corresponding cells
 383                     std::replace(&c2n_ptr[c2nIdx_ptr[cell] + 1], &c2n_ptr[c2nIdx_ptr[cell + 1]], node, i);
 384                     // This map is build in order to assign later new nodes to
 385                     // the corresponding faces
 386                     cellOld2NewNode[cell][node] = i;
 387                 }
 388                 i++;
 389             }
 390             i_compo++;
 391         }
 392     }
 393     return cellOld2NewNode;
 394 }
 395
 396 void
 397 CrackAlgo::AddMissingElementsOnLevelM1AndChangeConnectivity(
 398     const DataArrayIdType * f2cIdx,
 399     const DataArrayIdType * f2c,
 400     const DataArrayIdType * f2dupIdInM1,
 401     const DataArrayIdType * f2dupIdInMf,
 402     const Map2Map & cellOld2NewNode,
 403     MEDCouplingUMesh * m1,
 404     bool grpMustBeFullyDup
 405 ) {
 406     DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex();
 407     DataArrayIdType* f2n = m1->getNodalConnectivity();
 408
 409     // Change connectivity index size
 410     const auto f2nIdx_size = static_cast<size_t>(f2nIdx->getNumberOfTuples());  // = nb_face + 1
 411     const auto f2dup_size = static_cast<size_t>(f2dupIdInM1->getNumberOfTuples());
 412     if (f2dup_size != static_cast<size_t>(f2dupIdInMf->getNumberOfTuples()))
 413         throw INTERP_KERNEL::Exception(
 414             "There is not as many faces to dup looking into Mf and into M1.");
 415
 416     f2nIdx->reAlloc(f2nIdx_size + f2dup_size);
 417
 418     // Change connectivity size
 419     const auto f2n_size = static_cast<size_t>(f2n->getNumberOfTuples());
 420     const mcIdType * const f2nIdx_p = f2nIdx->begin();
 421     {
 422         size_t new_connectivity_size = f2n_size;
 423         for (const auto& face : *f2dupIdInM1)
 424             new_connectivity_size += static_cast<size_t>(
 425                 f2nIdx_p[face + 1] - f2nIdx_p[face]);
 426         f2n->reAlloc(new_connectivity_size);
 427     }
 428
 429     mcIdType * const f2nIdx_pw = f2nIdx->rwBegin();
 430     const mcIdType * const f2n_p = f2n->begin();
 431     mcIdType * const f2n_pw = f2n->rwBegin();
 432     size_t lastFace = f2nIdx_size - 1;
 433
 434     const mcIdType * const f2cIdx_p = f2cIdx->begin();
 435     const mcIdType * const f2c_p = f2c->begin();
 436
 437     for (size_t iFace = 0; iFace < f2dup_size; iFace++) {
 438         const auto & faceM1 = f2dupIdInM1->begin()[iFace];
 439         const auto & faceMf = f2dupIdInMf->begin()[iFace];
 440
 441         // Copy the face appending it
 442         //
 443         //   1. Append last connectivity index number
 444         const mcIdType face_start = f2nIdx_p[faceM1];
 445         const mcIdType face_end = f2nIdx_p[faceM1 + 1];
 446         const mcIdType lenFace = face_end - face_start;
 447         f2nIdx_pw[lastFace + 1] = f2nIdx_p[lastFace] + lenFace;
 448
 449         //   2. Append last face connectivity
 450         mcIdType *const new_f2n_p_start = f2n_pw + f2nIdx_p[lastFace];
 451         copy(f2n_p + face_start, f2n_p + face_end, new_f2n_p_start);
 452
 453         // Change inplace the connectivity
 454         // Check that this face is an inner face, ie it is connected to two
 455         // cells
 456         const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf];
 457         if (d != 2)
 458             throw INTERP_KERNEL::Exception(
 459                 "MEDFileMesh::duplicateFaces, the face to cell (or node to"
 460                 "segment) DataArray does not always adress two cells.");
 461
 462         bool noNodesAreChanged = true;
 463
 464         //   1. In original face, attaching it to cell0
 465         const auto& cell0 = f2c_p[f2cIdx_p[faceMf]];
 466         // If nodes where changed in cell0
 467         if (cellOld2NewNode.find(cell0) != cellOld2NewNode.end()) {
 468             const auto & mapO2N0 = cellOld2NewNode.at(cell0);
 469             for (mcIdType i_node = f2nIdx_p[faceM1] + 1; i_node < f2nIdx_p[faceM1 + 1]; i_node++) {
 470                 const mcIdType node = f2n_p[i_node];
 471                 // If this node was modified
 472                 if (mapO2N0.find(node) != mapO2N0.end()) {
 473                     const auto& new_node = mapO2N0.at(node);
 474                     f2n_pw[i_node] = new_node;
 475                     noNodesAreChanged = false;
 476                 }
 477             }
 478         }
 479
 480         //   2. In newly created face, attaching it to cell1
 481         const auto& cell1 = f2c_p[f2cIdx_p[faceMf] + 1];
 482         // If nodes where changed in cell1
 483         if (cellOld2NewNode.find(cell1) != cellOld2NewNode.end()) {
 484             const auto & mapO2N1 = cellOld2NewNode.at(cell1);
 485             for (mcIdType i_node = f2nIdx_p[lastFace] + 1; i_node < f2nIdx_p[lastFace + 1]; i_node++) {
 486                 const mcIdType node = f2n_p[i_node];
 487                 // If this node was modified
 488                 if (mapO2N1.find(node) != mapO2N1.end()) {
 489                     const auto& new_node = mapO2N1.at(node);
 490                     f2n_pw[i_node] = new_node;
 491                     noNodesAreChanged = false;
 492                 }
 493             }
 494         }
 495
 496         if (noNodesAreChanged && grpMustBeFullyDup)
 497             throw INTERP_KERNEL::Exception(
 498                 "duplicateFaces: A face was supposed to be duplicated but could"
 499                 "not be. Please be aware that all M1 groups cannot be"
 500                 "duplicated, at least two adjecent inner faces are needed to"
 501                 "duplicate a node.");
 502
 503         lastFace++;
 504     }
 505 }
 506
 507 CrackAlgo::Set
 508 CrackAlgo::GetCellsTouchingNodesToDup(
 509     const MEDCouplingUMesh * mf,
 510     const DataArrayIdType * n2cIdx,
 511     const DataArrayIdType * n2c,
 512     const DataArrayIdType * f2dup
 513 ) {
 514     Set res{};
 515
 516     const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex();
 517     const DataArrayIdType* f2n = mf->getNodalConnectivity();
 518
 519     const mcIdType * f2nIdx_p {f2nIdx->begin()};
 520     const mcIdType * f2n_p {f2n->begin()};
 521     const mcIdType * n2cIdx_p {n2cIdx->begin()};
 522     const mcIdType * n2c_p {n2c->begin()};
 523
 524     for (const auto face : *f2dup)
 525         for (mcIdType nodeIdx {f2nIdx_p[face] + 1}; nodeIdx < f2nIdx_p[face + 1]; nodeIdx++) {
 526             const auto &node = f2n_p[nodeIdx];
 527             for (mcIdType cellIdx = n2cIdx_p[node]; cellIdx < n2cIdx_p[node + 1]; cellIdx++) {
 528                 const auto &cell = n2c_p[cellIdx];
 529                 res.insert(cell);
 530             }
 531         }
 532     return res;
 533 }
 534
 535 CrackAlgo::Map2Set CrackAlgo::GetNode2CellMap(const MEDCouplingUMesh * mf,
 536         const DataArrayIdType * n2cIdx,
 537         const DataArrayIdType * n2c,
 538         const DataArrayIdType * f2dup
 539 ) {
 540     const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex();
 541     const DataArrayIdType* f2n = mf->getNodalConnectivity();
 542
 543     Map2Set n2c_dup{};
 544     const mcIdType * f2nIdx_p {f2nIdx->begin()};
 545     const mcIdType * f2n_p {f2n->begin()};
 546     const mcIdType * n2cIdx_p {n2cIdx->begin()};
 547     const mcIdType * n2c_p {n2c->begin()};
 548
 549     for (const auto face : *f2dup)
 550         for (mcIdType nodeIdx {f2nIdx_p[face] + 1}; nodeIdx < f2nIdx_p[face + 1]; nodeIdx++) {
 551             const auto &node = f2n_p[nodeIdx];
 552             for (mcIdType cellIdx {n2cIdx_p[node]}; cellIdx < n2cIdx_p[node + 1]; cellIdx++) {
 553                 const auto &cell = n2c_p[cellIdx];
 554                 n2c_dup[node].insert(cell);
 555             }
 556         }
 557     return n2c_dup;
 558 }
 559
 560 void
 561 CrackAlgo::OpenCrack(
 562     MEDFileUMesh * mf,
 563     const Map2Map & cellOld2NewNode,
 564     const double & factor
 565 ) {
 566     if ((factor <= 0.0) || (factor >= 1.0))
 567         throw INTERP_KERNEL::Exception("factor should be between 0.0 and 1.0");
 568     DataArrayDouble * coords = mf->getCoords();
 569     const auto dim = static_cast<mcIdType>(coords->getNumberOfComponents());
 570     const double * coords_p = coords->begin();
 571     MCU m0 = mf->getMeshAtLevel(0);
 572     DAD barys = m0->computeCellCenterOfMass();
 573     const double * barys_p = barys->begin();
 574     for (const auto & cellPair : cellOld2NewNode) {
 575         const auto & cell = cellPair.first;
 576         const auto & mapO2N = cellPair.second;
 577         for (const auto & nodeO2N : mapO2N) {
 578             const auto & oldN = nodeO2N.first;
 579             const auto & newN = nodeO2N.second;
 580             double * pos_new_p = coords->rwBegin() + dim * newN;
 581             for (int i = 0; i < dim; i++) {
 582                 pos_new_p[i] += (1.0 - factor) * (barys_p[cell*dim+i] - coords_p[dim*oldN+i]);
 583             }
 584         }
 585     }
 586 }
 587
 588 DataArrayIdType *
 589 CrackAlgo::GetFacesToDupInM1(
 590     const MEDCouplingUMesh & crackMesh,
 591     const MEDCouplingUMesh & m1
 592 ) {
 593     DataArrayIdType * f2dupIdInM1P{};
 594     const bool allIn = m1.areCellsIncludedIn(&crackMesh, 2, f2dupIdInM1P);
 595     if (!allIn)
 596         throw INTERP_KERNEL::Exception(
 597             "crackAlong: cleanCrackMesh is not part of m1 mesh.");
 598     return f2dupIdInM1P;
 599 }
 600
 601 pair<DataArrayIdType*, DataArrayIdType*>
 602 CrackAlgo::GetFacesInM1TouchingDuplicatedNodes(
 603     const Map2Set & n2c_dup,
 604     const DataArrayIdType * f2dupIdInM1,
 605     const MEDCouplingUMesh & mf,
 606     const MEDCouplingUMesh & m1
 607 ) {
 608     Set f2dup_set = DataArrayToSet(*f2dupIdInM1);
 609
 610     DAI n2f(DataArrayIdType::New());
 611     DAI n2fIdx(DataArrayIdType::New());
 612     m1.getReverseNodalConnectivity(n2f, n2fIdx);
 613
 614     // 1. I retrieve all the faces of m1 which have a potentially duplicated
 615     // node, I just use the reverse nodal connectivity of m1 for this.
 616     Set f2changeIdInM1_set{};
 617
 618     for (const auto & nodeP : n2c_dup) {
 619         const auto & node = nodeP.first;
 620         for (mcIdType fId = n2fIdx->begin()[node]; fId < n2fIdx->begin()[node+1]; fId++) {
 621             const auto & f = n2f->begin()[fId];
 622             // f is not a face to duplicate
 623             if (f2dup_set.find(f) == f2dup_set.end())
 624                 f2changeIdInM1_set.insert(f);
 625         }
 626     }
 627     // 2. I retrieve M1 sub mesh
 628     DataArrayIdType * f2changeIdInM1 = SetToDataArray(f2changeIdInM1_set);
 629     MCU part = m1.buildPartOfMySelf(f2changeIdInM1->begin(), f2changeIdInM1->end());
 630     // 3. I search for the ids in mf
 631     DataArrayIdType * f2changeIdInMf{};
 632     const bool ok = mf.areCellsIncludedIn(part, 2, f2changeIdInMf);
 633     if (!ok)
 634         throw INTERP_KERNEL::Exception(
 635             "Some faces in -1 mesh next to the crack are not included in "
 636             "the face mesh issued from m0.");
 637
 638     // 4. I return the 2 DAI
 639     return pair<DataArrayIdType*, DataArrayIdType*>{f2changeIdInM1, f2changeIdInMf};
 640 }
 641
 642 void
 643 CrackAlgo::ChangeConnectivityOfM1Elements(
 644     const DataArrayIdType * f2changeIdInM1,
 645     const DataArrayIdType * f2changeIdInMf,
 646     const Map2Map & cellOld2NewNode,
 647     const DataArrayIdType * f2cIdx,
 648     const DataArrayIdType * f2c,
 649     MEDCouplingUMesh * m1
 650 ) {
 651     // 1. I go through the faces in question. I go through all the nodes on the
 652     // face. I will look in neighboring cells to see if there have been changes
 653     // in connectivity on the nodes of the face.
 654     // 2. if the changes are consistent, I change in place.
 655     // 3. otherwise I return a consistency error.
 656
 657     DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex();
 658     DataArrayIdType* f2n = m1->getNodalConnectivity();
 659
 660     // Change connectivity size
 661     const mcIdType * const f2nIdx_p = f2nIdx->begin();
 662     const mcIdType * const f2n_p = f2n->begin();
 663     mcIdType * const f2n_pw = f2n->rwBegin();
 664
 665     const mcIdType * const f2cIdx_p {f2cIdx->begin()};
 666     const mcIdType * const f2c_p {f2c->begin()};
 667
 668     const auto f2change_size = static_cast<size_t>(f2changeIdInM1->getNumberOfTuples());
 669
 670     for (size_t iFace = 0; iFace < f2change_size; iFace++) {
 671         const auto & faceM1 = f2changeIdInM1->begin()[iFace];
 672         const auto & faceMf = f2changeIdInMf->begin()[iFace];
 673
 674         // Change inplace the connectivity
 675         // Check that this face is an inner face, ie it is connected to two
 676         // cells
 677         const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf];
 678
 679         //   1. In original face, attaching it to cell0
 680         const auto& cell0 = f2c_p[f2cIdx_p[faceMf]];
 681         // const auto& cell1 = f2c_p[f2cIdx_p[faceMf] + 1];
 682
 683         // If nodes where changed in cell0
 684         if (cellOld2NewNode.find(cell0) != cellOld2NewNode.end()) {
 685             const auto & mapO2N0 = cellOld2NewNode.at(cell0);
 686             for (mcIdType i_node = f2nIdx_p[faceM1] + 1; i_node < f2nIdx_p[faceM1 + 1]; i_node++) {
 687                 const mcIdType node = f2n_p[i_node];
 688                 // If this node was modified
 689                 if (mapO2N0.find(node) != mapO2N0.end()) {
 690                     const auto& new_node = mapO2N0.at(node);
 691                     if ((d == 2)
 692                         && ((cellOld2NewNode.find(f2c_p[f2cIdx_p[faceMf] + 1])
 693                                 == cellOld2NewNode.end())
 694                             || (cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).find(node)
 695                                 == cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).end())
 696                             || (cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).at(node) != new_node)))
 697                             throw INTERP_KERNEL::Exception("There is an"
 698                                 "incoherent change of connectivity between both"
 699                                 "cell surrounding a face touching the crack but"
 700                                 "not duplicated.");
 701                     f2n_pw[i_node] = new_node;
 702                 }
 703             }
 704         }
 705     }
 706 }