1 // Copyright (C) 2007-2024 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email :
18 // webmaster.salome@opencascade.com
20 // Author : Aymeric SONOLET (CEA/DES)
22 #include "CrackAlgo.hxx"
31 #include "./MEDFileMesh.hxx"
33 #include <MEDCouplingUMesh.hxx>
34 #include <MEDCouplingMemArray.hxx>
35 #include <InterpKernelException.hxx>
39 using namespace MEDCoupling;
41 using MCU = MCAuto<MEDCouplingUMesh>;
42 using DAI = MCAuto<DataArrayIdType>;
43 using DAD = MCAuto<DataArrayDouble>;
49 const string& grp_name,
50 bool grpMustBeFullyDup
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 !");
60 MCU m0 = mm->getMeshAtLevel(0);
62 MCU crackingMesh = mm->getGroup(-1, grp_name);
63 MCU cleanCrackingMesh = CleanM1Mesh(*m0, *crackingMesh);
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);
71 DataArrayIdType * crackToMfP;
72 const bool crackCellsInMf = mf->areCellsIncludedIn(cleanCrackingMesh, 2, crackToMfP);
73 DAI f2dupIdInMf(crackToMfP);
75 throw INTERP_KERNEL::Exception(
76 "crackAlong: All cells in crack are not part of Mf.");
78 DAI n2c(DataArrayIdType::New()),
79 n2cIdx(DataArrayIdType::New());
80 m0->getReverseNodalConnectivity(n2c, n2cIdx);
82 const Map2Set n2c_dup = GetNode2CellMap(mf, n2cIdx, n2c, f2dupIdInMf);
84 const Set cTouchingN_dup = GetCellsTouchingNodesToDup(
85 mf, n2cIdx, n2c, f2dupIdInMf);
87 const Graph c2c = BuildCutC2CGraph(
88 c2fIdx, c2f, f2cIdx, f2c, cTouchingN_dup, f2dupIdInMf);
90 Map2Map cellOld2NewNode = CreateNewNodesInTopLevelMesh(n2c_dup, c2c, m0);
91 // End of node creation, separation of M0
93 // Faces in M1 must be added/updated
95 MCU m1 = mm->getMeshAtLevel(-1);
98 const DAI f2dupIdInM1 = GetFacesToDupInM1(*cleanCrackingMesh, *m1);
100 const auto f2changeM1Mf = GetFacesInM1TouchingDuplicatedNodes(
101 n2c_dup, f2dupIdInM1, *mf, *m1);
102 const DAI f2changeIdInM1 = f2changeM1Mf.first;
103 const DAI f2changeIdInMf = f2changeM1Mf.second;
105 AddMissingElementsOnLevelM1AndChangeConnectivity(
106 f2cIdx, f2c, f2dupIdInM1, f2dupIdInMf, cellOld2NewNode, m1, grpMustBeFullyDup);
108 ChangeConnectivityOfM1Elements(
109 f2changeIdInM1, f2changeIdInMf, cellOld2NewNode, f2cIdx, f2c, m1);
111 // TODO(aymeric): If one wants to implement
112 // AddMissingElementsOnLevelM2AndChangeConnectivity and
113 // ChangeConnectivityOfM2Elements it should be there.
115 DAI famLev0 = CopyFamilyArrAtLev0(mm);
116 DAI newFam = CopyAndCompleteFamilyArrAtLevelM1(mm, m1, f2dupIdInM1);
118 mm->setMeshAtLevel(0, m0);
119 mm->setMeshAtLevel(-1, m1);
120 // mm->setMeshAtLevel(-2, m2);
122 mm->setFamilyFieldArr(0, famLev0);
123 mm->setFamilyFieldArr(-1, newFam);
124 // mm->setFamilyFieldArr(-2, newFam2);
126 const Map2Set addedNodes = BuildMap2Set(cellOld2NewNode);
127 CompleteFamilyArrAtNodeLevel(addedNodes, mm);
129 return cellOld2NewNode;
133 CrackAlgo::BuildMap2Set(
134 const Map2Map& cellOld2NewNode
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);
147 CrackAlgo::CopyAndCompleteFamilyArrAtLevelM1(
148 const MEDFileUMesh * mm,
149 const MEDCouplingUMesh * mf,
150 const DataArrayIdType * f2dup
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];
171 CrackAlgo::CopyFamilyArrAtLev0(
172 const MEDFileUMesh * mm
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());
186 CrackAlgo::CompleteFamilyArrAtNodeLevel(
187 const Map2Set& addedNodes,
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];
205 CrackAlgo::CleanM1Mesh(
206 const MEDCouplingUMesh & m,
207 const MEDCouplingUMesh & crackingMesh
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();
222 const mcIdType &node,
223 const size_t &componentId,
224 map<mcIdType, bool>& visited,
225 map<mcIdType, size_t>& componentMap
227 visited[node] = true;
228 // Assign a unique component ID to the node
229 componentMap[node] = componentId;
231 for (const auto& neighbor : graph.at(node)) {
232 if (!visited[neighbor]) {
233 Dfs(graph, neighbor, componentId, visited, componentMap);
238 // Function to find connected components
239 vector<shared_ptr<vector<mcIdType>>>
240 CrackAlgo::FindConnectedComponents(
243 map<mcIdType, bool> visited;
244 map<mcIdType, size_t> componentMap;
245 size_t componentId = 0;
247 for (const auto &nodePair : graph)
248 visited[nodePair.first] = false;
250 for (const auto& nodePair : graph) {
251 if (!visited[nodePair.first]) {
252 Dfs(graph, nodePair.first, componentId, visited, componentMap);
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);
268 CrackAlgo::Set CrackAlgo::DataArrayToSet(
269 const DataArrayIdType& da
272 for (const auto &elem : da) {
278 DataArrayIdType * CrackAlgo::SetToDataArray(
281 DataArrayIdType * da(DataArrayIdType::New());
283 mcIdType * da_p = da->rwBegin();
285 for (const auto & elem : s) {
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
302 Set f2dup_set = DataArrayToSet(*f2dup);
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()};
310 for (const auto &cell : cTouchingN_dup) {
311 c2c[cell] = 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] };
322 & (cTouchingN_dup.find(cell2) != cTouchingN_dup.end())
324 c2c.at(cell).insert(cell2);
334 CrackAlgo::CreateNewNodesInTopLevelMesh(
335 const Map2Set & n2c_dup,
339 DataArrayIdType* c2n = m0->getNodalConnectivity();
340 DataArrayIdType* c2nIdx = m0->getNodalConnectivityIndex();
341 mcIdType * c2n_ptr = c2n->rwBegin();
342 const mcIdType * c2nIdx_ptr = c2nIdx->begin();
344 DataArrayDouble* coords = m0->getCoords();
345 mcIdType i = coords->getNumberOfTuples();
346 const size_t coordsDim = coords->getNumberOfComponents();
348 Map2Map cellOld2NewNode;
350 for (const auto &nc_pair : n2c_dup) {
351 const auto &node = nc_pair.first;
352 const auto &cells = nc_pair.second;
354 // Building local connection graph
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);
367 const vector<shared_ptr<vector<mcIdType>>> compo_connex = FindConnectedComponents(c2c_local);
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);
374 for (const auto &compo : compo_connex) {
376 // Coords are copied at the end of the vector of coords, ie the
377 // node is duplicated
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 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;
393 return cellOld2NewNode;
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
406 DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex();
407 DataArrayIdType* f2n = m1->getNodalConnectivity();
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.");
416 f2nIdx->reAlloc(f2nIdx_size + f2dup_size);
418 // Change connectivity size
419 const auto f2n_size = static_cast<size_t>(f2n->getNumberOfTuples());
420 const mcIdType * const f2nIdx_p = f2nIdx->begin();
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);
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;
434 const mcIdType * const f2cIdx_p = f2cIdx->begin();
435 const mcIdType * const f2c_p = f2c->begin();
437 for (size_t iFace = 0; iFace < f2dup_size; iFace++) {
438 const auto & faceM1 = f2dupIdInM1->begin()[iFace];
439 const auto & faceMf = f2dupIdInMf->begin()[iFace];
441 // Copy the face appending it
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;
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);
453 // Change inplace the connectivity
454 // Check that this face is an inner face, ie it is connected to two
456 const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf];
458 throw INTERP_KERNEL::Exception(
459 "MEDFileMesh::duplicateFaces, the face to cell (or node to"
460 "segment) DataArray does not always adress two cells.");
462 bool noNodesAreChanged = true;
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;
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;
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.");
508 CrackAlgo::GetCellsTouchingNodesToDup(
509 const MEDCouplingUMesh * mf,
510 const DataArrayIdType * n2cIdx,
511 const DataArrayIdType * n2c,
512 const DataArrayIdType * f2dup
516 const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex();
517 const DataArrayIdType* f2n = mf->getNodalConnectivity();
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()};
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];
535 CrackAlgo::Map2Set CrackAlgo::GetNode2CellMap(const MEDCouplingUMesh * mf,
536 const DataArrayIdType * n2cIdx,
537 const DataArrayIdType * n2c,
538 const DataArrayIdType * f2dup
540 const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex();
541 const DataArrayIdType* f2n = mf->getNodalConnectivity();
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()};
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);
561 CrackAlgo::OpenCrack(
563 const Map2Map & cellOld2NewNode,
564 const double & factor
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]);
589 CrackAlgo::GetFacesToDupInM1(
590 const MEDCouplingUMesh & crackMesh,
591 const MEDCouplingUMesh & m1
593 DataArrayIdType * f2dupIdInM1P{};
594 const bool allIn = m1.areCellsIncludedIn(&crackMesh, 2, f2dupIdInM1P);
596 throw INTERP_KERNEL::Exception(
597 "crackAlong: cleanCrackMesh is not part of m1 mesh.");
601 pair<DataArrayIdType*, DataArrayIdType*>
602 CrackAlgo::GetFacesInM1TouchingDuplicatedNodes(
603 const Map2Set & n2c_dup,
604 const DataArrayIdType * f2dupIdInM1,
605 const MEDCouplingUMesh & mf,
606 const MEDCouplingUMesh & m1
608 Set f2dup_set = DataArrayToSet(*f2dupIdInM1);
610 DAI n2f(DataArrayIdType::New());
611 DAI n2fIdx(DataArrayIdType::New());
612 m1.getReverseNodalConnectivity(n2f, n2fIdx);
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{};
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);
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);
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.");
638 // 4. I return the 2 DAI
639 return pair<DataArrayIdType*, DataArrayIdType*>{f2changeIdInM1, f2changeIdInMf};
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
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.
657 DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex();
658 DataArrayIdType* f2n = m1->getNodalConnectivity();
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();
665 const mcIdType * const f2cIdx_p {f2cIdx->begin()};
666 const mcIdType * const f2c_p {f2c->begin()};
668 const auto f2change_size = static_cast<size_t>(f2changeIdInM1->getNumberOfTuples());
670 for (size_t iFace = 0; iFace < f2change_size; iFace++) {
671 const auto & faceM1 = f2changeIdInM1->begin()[iFace];
672 const auto & faceMf = f2changeIdInMf->begin()[iFace];
674 // Change inplace the connectivity
675 // Check that this face is an inner face, ie it is connected to two
677 const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf];
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];
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);
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"
701 f2n_pw[i_node] = new_node;