From 3ccbd6672d4cab9ddd873774deb1276b43ccb621 Mon Sep 17 00:00:00 2001 From: SONOLET Aymeric Date: Mon, 27 May 2024 08:20:33 +0200 Subject: [PATCH] feat: new crackAlong method This new method is a more powerfull version of the buildInnerBoundariesAlongM1Group method because it naturally handle non connex and complex M1 groups wip: fix cycling over nodal connectivity fix shared_ptr> allocation wip: no more memleak ! TODO: for now this method only treat the mesh at level 0. The mesh at level -1 and the groups lying on it should be taken care of. wip: Duplicate faces inplace and replace old node number with new ones wip: removing useless shared_ptr, adding doc wip: moving duplicateFaces to CrackAlgo.cxx wip: refactor CrackAlgo and treat elements families missing node families feat: functional new duplicateFaces method fix: changing API and preparing methods for test feat: separate tests on crackAlong method wip: using M1 and not Mf fix: case when cell has no neighbor feat: test mesh at level -1 This test checks that the crack mesh has twice as many cells than the starting crack mesh, and that they are both equivalent after merging the nodes. wip: fixing connectivity of faces which are not duplicated fix: manage connectivity of m1 elements touching crack wip: add python tests Still figuring out why using M1 non connex group erease the families at lev -2 fix: crackAlong exceptions with boolean param A boolean is added to change the method behavior in the case when a face is duplicated but not modified (there is no real separation between volumes) feat: new tests and OpenCrack binding fix: Copy family at lev0 and set them back Because suprisingly they tend to disapear. refactor: OpenCrack -> openCrack fix: remove french comments, better use of stl refactor: better documentation, less auto, more const --- src/MEDLoader/CMakeLists.txt | 1 + src/MEDLoader/CrackAlgo.cxx | 706 +++++++++++++++++++++++++++ src/MEDLoader/CrackAlgo.hxx | 193 ++++++++ src/MEDLoader/MEDFileMesh.cxx | 51 ++ src/MEDLoader/MEDFileMesh.hxx | 16 + src/MEDLoader/Swig/CMakeLists.txt | 3 +- src/MEDLoader/Swig/CrackAlongTest.py | 648 ++++++++++++++++++++++++ src/MEDLoader/Swig/MEDLoaderCommon.i | 54 ++ src/MEDLoader/Swig/tests.set | 1 + src/MEDLoader/Test/CMakeLists.txt | 1 + src/MEDLoader/Test/CrackAlgoTest.cxx | 657 +++++++++++++++++++++++++ src/MEDLoader/Test/CrackAlgoTest.hxx | 108 ++++ src/MEDLoader/Test/MEDLoaderTest.cxx | 7 +- src/MEDLoader/Test/MEDLoaderTest.hxx | 2 + src/MEDLoader/Test/TestMEDLoader.cxx | 2 + 15 files changed, 2445 insertions(+), 5 deletions(-) create mode 100644 src/MEDLoader/CrackAlgo.cxx create mode 100644 src/MEDLoader/CrackAlgo.hxx create mode 100644 src/MEDLoader/Swig/CrackAlongTest.py create mode 100644 src/MEDLoader/Test/CrackAlgoTest.cxx create mode 100644 src/MEDLoader/Test/CrackAlgoTest.hxx diff --git a/src/MEDLoader/CMakeLists.txt b/src/MEDLoader/CMakeLists.txt index aec75954d..6062d3062 100644 --- a/src/MEDLoader/CMakeLists.txt +++ b/src/MEDLoader/CMakeLists.txt @@ -59,6 +59,7 @@ INCLUDE_DIRECTORIES( ) SET(medloader_SOURCES + CrackAlgo.cxx MEDLoader.cxx MEDLoaderBase.cxx MEDLoaderTraits.cxx diff --git a/src/MEDLoader/CrackAlgo.cxx b/src/MEDLoader/CrackAlgo.cxx new file mode 100644 index 000000000..836c14f18 --- /dev/null +++ b/src/MEDLoader/CrackAlgo.cxx @@ -0,0 +1,706 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : +// webmaster.salome@opencascade.com +// +// Author : Aymeric SONOLET (CEA/DES) + +#include "CrackAlgo.hxx" + +#include +#include +#include +#include +#include +#include + +#include "./MEDFileMesh.hxx" + +#include +#include +#include +#include +#include + +using namespace MEDCoupling; +using namespace std; +using MCU = MCAuto; +using DAI = MCAuto; +using DAD = MCAuto; + + +CrackAlgo::Map2Map +CrackAlgo::Compute( + MEDFileUMesh* mm, + const string& grp_name, + bool grpMustBeFullyDup +) { + + vector levs = mm->getNonEmptyLevels(); + if (find(levs.begin(), levs.end(), 0) == levs.end() || + find(levs.begin(), levs.end(), -1) == levs.end()) + throw INTERP_KERNEL::Exception( + "MEDFileUMesh::duplicateFaces : This method " + "works only for mesh defined on level 0 and -1 !"); + + MCU m0 = mm->getMeshAtLevel(0); + + MCU crackingMesh = mm->getGroup(-1, grp_name); + MCU cleanCrackingMesh = CleanM1Mesh(*m0, *crackingMesh); + + DAI c2f(DataArrayIdType::New()), + c2fIdx(DataArrayIdType::New()), + f2c(DataArrayIdType::New()), + f2cIdx(DataArrayIdType::New()); + const MCU mf = m0->buildDescendingConnectivity(c2f, c2fIdx, f2c, f2cIdx); + + DataArrayIdType * crackToMfP; + const bool crackCellsInMf = mf->areCellsIncludedIn(cleanCrackingMesh, 2, crackToMfP); + DAI f2dupIdInMf(crackToMfP); + if (!crackCellsInMf) + throw INTERP_KERNEL::Exception( + "crackAlong: All cells in crack are not part of Mf."); + + DAI n2c(DataArrayIdType::New()), + n2cIdx(DataArrayIdType::New()); + m0->getReverseNodalConnectivity(n2c, n2cIdx); + + const Map2Set n2c_dup = GetNode2CellMap(mf, n2cIdx, n2c, f2dupIdInMf); + + const Set cTouchingN_dup = GetCellsTouchingNodesToDup( + mf, n2cIdx, n2c, f2dupIdInMf); + + const Graph c2c = BuildCutC2CGraph( + c2fIdx, c2f, f2cIdx, f2c, cTouchingN_dup, f2dupIdInMf); + + Map2Map cellOld2NewNode = CreateNewNodesInTopLevelMesh(n2c_dup, c2c, m0); + // End of node creation, separation of M0 + + // Faces in M1 must be added/updated + + MCU m1 = mm->getMeshAtLevel(-1); + + // Before changing M1 + const DAI f2dupIdInM1 = GetFacesToDupInM1(*cleanCrackingMesh, *m1); + + const auto f2changeM1Mf = GetFacesInM1TouchingDuplicatedNodes( + n2c_dup, f2dupIdInM1, *mf, *m1); + const DAI f2changeIdInM1 = f2changeM1Mf.first; + const DAI f2changeIdInMf = f2changeM1Mf.second; + + AddMissingElementsOnLevelM1AndChangeConnectivity( + f2cIdx, f2c, f2dupIdInM1, f2dupIdInMf, cellOld2NewNode, m1, grpMustBeFullyDup); + + ChangeConnectivityOfM1Elements( + f2changeIdInM1, f2changeIdInMf, cellOld2NewNode, f2cIdx, f2c, m1); + + // TODO(aymeric): If one wants to implement + // AddMissingElementsOnLevelM2AndChangeConnectivity and + // ChangeConnectivityOfM2Elements it should be there. + + DAI famLev0 = CopyFamilyArrAtLev0(mm); + DAI newFam = CopyAndCompleteFamilyArrAtLevelM1(mm, m1, f2dupIdInM1); + + mm->setMeshAtLevel(0, m0); + mm->setMeshAtLevel(-1, m1); + // mm->setMeshAtLevel(-2, m2); + + mm->setFamilyFieldArr(0, famLev0); + mm->setFamilyFieldArr(-1, newFam); + // mm->setFamilyFieldArr(-2, newFam2); + + const Map2Set addedNodes = BuildMap2Set(cellOld2NewNode); + CompleteFamilyArrAtNodeLevel(addedNodes, mm); + + return cellOld2NewNode; +} + +CrackAlgo::Map2Set +CrackAlgo::BuildMap2Set( + const Map2Map& cellOld2NewNode +) { + Map2Set res; + for (const auto & keyElem : cellOld2NewNode) { + const auto & old2NewMap = keyElem.second; + for (const auto & old2New : old2NewMap) { + res[old2New.first].insert(old2New.second); + } + } + return res; +} + +DataArrayIdType * +CrackAlgo::CopyAndCompleteFamilyArrAtLevelM1( + const MEDFileUMesh * mm, + const MEDCouplingUMesh * mf, + const DataArrayIdType * f2dup +) { + DataArrayIdType* newFam = nullptr; + const DataArrayIdType *fam = mm->getFamilyFieldAtLevel(-1); + if (fam != nullptr) { + newFam = DataArrayIdType::New(); + const auto fam_size = static_cast(fam->getNumberOfTuples()); + newFam->alloc(static_cast(mf->getNumberOfCells())); + copy(fam->begin(), fam->end(), newFam->getPointer()); + mcIdType * fam_ptr = newFam->getPointer(); + size_t i_elem_added = 0; + for (const auto & face : *f2dup) { + // NOTE(aymeric): assumes all faces are duplicated + fam_ptr[fam_size + i_elem_added] = fam_ptr[face]; + i_elem_added++; + } + } + return newFam; +} + +DataArrayIdType * +CrackAlgo::CopyFamilyArrAtLev0( + const MEDFileUMesh * mm +) { + DataArrayIdType* famLev0 = nullptr; + const DataArrayIdType *fam = mm->getFamilyFieldAtLevel(0); + if (fam != nullptr) { + famLev0 = DataArrayIdType::New(); + const auto fam_size = static_cast(fam->getNumberOfTuples()); + famLev0->alloc(fam_size); + copy(fam->begin(), fam->end(), famLev0->getPointer()); + } + return famLev0; +} + +void +CrackAlgo::CompleteFamilyArrAtNodeLevel( + const Map2Set& addedNodes, + MEDFileUMesh * mm +) { + DataArrayIdType *node_fam = mm->getFamilyFieldAtLevel(1); + if (node_fam != nullptr) { + node_fam->reAlloc(static_cast(mm->getCoords()->getNumberOfTuples())); + mcIdType * node_fam_ptr = node_fam->getPointer(); + for (const auto & old2NewPair : addedNodes) { + const auto & oldNode = old2NewPair.first; + const auto & newNodes = old2NewPair.second; + for (const auto & newNode : newNodes) { + node_fam_ptr[newNode] = node_fam_ptr[oldNode]; + } + } + } +} + +MEDCouplingUMesh * +CrackAlgo::CleanM1Mesh( + const MEDCouplingUMesh & m, + const MEDCouplingUMesh & crackingMesh +) { + MCU m0skin = m.computeSkin(); + DataArrayIdType * idsToKeepP; + m0skin->areCellsIncludedIn(&crackingMesh, 2, idsToKeepP); + DAI idsToKeep(idsToKeepP); + // discard cells on the skin of M0 + DAI ids2 = idsToKeep->findIdsNotInRange(0, m0skin->getNumberOfCells()); + MCU otherDimM1OnSameCoords = + crackingMesh.buildPartOfMySelf(ids2->begin(), ids2->end(), true); + return otherDimM1OnSameCoords.retn(); +} + +void CrackAlgo::Dfs( + const Graph& graph, + const mcIdType &node, + const size_t &componentId, + map& visited, + map& componentMap +) { + visited[node] = true; + // Assign a unique component ID to the node + componentMap[node] = componentId; + + for (const auto& neighbor : graph.at(node)) { + if (!visited[neighbor]) { + Dfs(graph, neighbor, componentId, visited, componentMap); + } + } +} + +// Function to find connected components +vector>> +CrackAlgo::FindConnectedComponents( + const Graph& graph +) { + map visited; + map componentMap; + size_t componentId = 0; + + for (const auto &nodePair : graph) + visited[nodePair.first] = false; + + for (const auto& nodePair : graph) { + if (!visited[nodePair.first]) { + Dfs(graph, nodePair.first, componentId, visited, componentMap); + componentId++; + } + } + + using spv = shared_ptr>; + vector components; + for (size_t i = 0; i < componentId; i++) + components.push_back(make_shared>()); + for (const auto& nodePair : componentMap) { + components[nodePair.second]->push_back(nodePair.first); + } + + return components; +} + +CrackAlgo::Set CrackAlgo::DataArrayToSet( + const DataArrayIdType& da +) { + Set res; + for (const auto &elem : da) { + res.insert(elem); + } + return res; +} + +DataArrayIdType * CrackAlgo::SetToDataArray( + const Set& s +) { + DataArrayIdType * da(DataArrayIdType::New()); + da->alloc(s.size()); + mcIdType * da_p = da->rwBegin(); + size_t i = 0; + for (const auto & elem : s) { + da_p[i] = elem; + i++; + } + return da; +} + +CrackAlgo::Graph +CrackAlgo::BuildCutC2CGraph( + const DataArrayIdType * c2fIdx, + const DataArrayIdType * c2f, + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + const Set & cTouchingN_dup, + const DataArrayIdType * f2dup +) { + + Set f2dup_set = DataArrayToSet(*f2dup); + Graph c2c; + + const mcIdType * c2fIdx_p {c2fIdx->begin()}; + const mcIdType * c2f_p {c2f->begin()}; + const mcIdType * f2cIdx_p {f2cIdx->begin()}; + const mcIdType * f2c_p {f2c->begin()}; + + for (const auto &cell : cTouchingN_dup) { + c2c[cell] = std::unordered_set(); + for (auto faceIdx = c2fIdx_p[cell]; faceIdx < c2fIdx_p[cell + 1]; faceIdx++) { + const mcIdType &face = c2f_p[faceIdx]; + // face is not in face to duplicate + if (f2dup_set.find(face) == f2dup_set.end()) { + for (mcIdType cellIdx {f2cIdx_p[face]}; cellIdx < f2cIdx_p[face +1]; cellIdx++) { + // getting the id of the other cell (on face is + // connected to two cells) + const auto& cell2 { f2c_p[cellIdx] }; + if ( + (cell != cell2) + & (cTouchingN_dup.find(cell2) != cTouchingN_dup.end()) + ) { + c2c.at(cell).insert(cell2); + } + } + } + } + } + return c2c; +} + +CrackAlgo::Map2Map +CrackAlgo::CreateNewNodesInTopLevelMesh( + const Map2Set & n2c_dup, + const Graph & c2c, + MEDCouplingUMesh* m0 +) { + DataArrayIdType* c2n = m0->getNodalConnectivity(); + DataArrayIdType* c2nIdx = m0->getNodalConnectivityIndex(); + mcIdType * c2n_ptr = c2n->rwBegin(); + const mcIdType * c2nIdx_ptr = c2nIdx->begin(); + + DataArrayDouble* coords = m0->getCoords(); + mcIdType i = coords->getNumberOfTuples(); + const size_t coordsDim = coords->getNumberOfComponents(); + + Map2Map cellOld2NewNode; + + for (const auto &nc_pair : n2c_dup) { + const auto &node = nc_pair.first; + const auto &cells = nc_pair.second; + + // Building local connection graph + Graph c2c_local; + for (const auto &cell : cells) { + auto& val = c2c_local[cell]; + if (c2c.find(cell) == c2c.end()) + throw INTERP_KERNEL::Exception( + "A cell touching a node is not part of c2c."); + for (const auto &neighbor : c2c.at(cell)) { + if (cells.find(neighbor) != cells.end()) + val.insert(neighbor); + } + } + + const vector>> compo_connex = FindConnectedComponents(c2c_local); + + // Duplicate node for compo[1:] + const size_t i_nodes_to_add = compo_connex.size() - 1; + coords->reAlloc(static_cast(i) + i_nodes_to_add); + + int i_compo = 0; + for (const auto &compo : compo_connex) { + if (i_compo > 0) { + // Coords are copied at the end of the vector of coords, ie the + // node is duplicated + coords->getTuple( + node, + coords->getPointer() + static_cast(i) * coordsDim); + for (const auto &cell : *compo) { + // The node number is replaced in all corresponding cells + std::replace(&c2n_ptr[c2nIdx_ptr[cell] + 1], &c2n_ptr[c2nIdx_ptr[cell + 1]], node, i); + // This map is build in order to assign later new nodes to + // the corresponding faces + cellOld2NewNode[cell][node] = i; + } + i++; + } + i_compo++; + } + } + return cellOld2NewNode; +} + +void +CrackAlgo::AddMissingElementsOnLevelM1AndChangeConnectivity( + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + const DataArrayIdType * f2dupIdInM1, + const DataArrayIdType * f2dupIdInMf, + const Map2Map & cellOld2NewNode, + MEDCouplingUMesh * m1, + bool grpMustBeFullyDup +) { + DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex(); + DataArrayIdType* f2n = m1->getNodalConnectivity(); + + // Change connectivity index size + const auto f2nIdx_size = static_cast(f2nIdx->getNumberOfTuples()); // = nb_face + 1 + const auto f2dup_size = static_cast(f2dupIdInM1->getNumberOfTuples()); + if (f2dup_size != static_cast(f2dupIdInMf->getNumberOfTuples())) + throw INTERP_KERNEL::Exception( + "There is not as many faces to dup looking into Mf and into M1."); + + f2nIdx->reAlloc(f2nIdx_size + f2dup_size); + + // Change connectivity size + const auto f2n_size = static_cast(f2n->getNumberOfTuples()); + const mcIdType * const f2nIdx_p = f2nIdx->begin(); + { + size_t new_connectivity_size = f2n_size; + for (const auto& face : *f2dupIdInM1) + new_connectivity_size += static_cast( + f2nIdx_p[face + 1] - f2nIdx_p[face]); + f2n->reAlloc(new_connectivity_size); + } + + mcIdType * const f2nIdx_pw = f2nIdx->rwBegin(); + const mcIdType * const f2n_p = f2n->begin(); + mcIdType * const f2n_pw = f2n->rwBegin(); + size_t lastFace = f2nIdx_size - 1; + + const mcIdType * const f2cIdx_p = f2cIdx->begin(); + const mcIdType * const f2c_p = f2c->begin(); + + for (size_t iFace = 0; iFace < f2dup_size; iFace++) { + const auto & faceM1 = f2dupIdInM1->begin()[iFace]; + const auto & faceMf = f2dupIdInMf->begin()[iFace]; + + // Copy the face appending it + // + // 1. Append last connectivity index number + const mcIdType face_start = f2nIdx_p[faceM1]; + const mcIdType face_end = f2nIdx_p[faceM1 + 1]; + const mcIdType lenFace = face_end - face_start; + f2nIdx_pw[lastFace + 1] = f2nIdx_p[lastFace] + lenFace; + + // 2. Append last face connectivity + mcIdType *const new_f2n_p_start = f2n_pw + f2nIdx_p[lastFace]; + copy(f2n_p + face_start, f2n_p + face_end, new_f2n_p_start); + + // Change inplace the connectivity + // Check that this face is an inner face, ie it is connected to two + // cells + const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf]; + if (d != 2) + throw INTERP_KERNEL::Exception( + "MEDFileMesh::duplicateFaces, the face to cell (or node to" + "segment) DataArray does not always adress two cells."); + + bool noNodesAreChanged = true; + + // 1. In original face, attaching it to cell0 + const auto& cell0 = f2c_p[f2cIdx_p[faceMf]]; + // If nodes where changed in cell0 + if (cellOld2NewNode.find(cell0) != cellOld2NewNode.end()) { + const auto & mapO2N0 = cellOld2NewNode.at(cell0); + for (mcIdType i_node = f2nIdx_p[faceM1] + 1; i_node < f2nIdx_p[faceM1 + 1]; i_node++) { + const mcIdType node = f2n_p[i_node]; + // If this node was modified + if (mapO2N0.find(node) != mapO2N0.end()) { + const auto& new_node = mapO2N0.at(node); + f2n_pw[i_node] = new_node; + noNodesAreChanged = false; + } + } + } + + // 2. In newly created face, attaching it to cell1 + const auto& cell1 = f2c_p[f2cIdx_p[faceMf] + 1]; + // If nodes where changed in cell1 + if (cellOld2NewNode.find(cell1) != cellOld2NewNode.end()) { + const auto & mapO2N1 = cellOld2NewNode.at(cell1); + for (mcIdType i_node = f2nIdx_p[lastFace] + 1; i_node < f2nIdx_p[lastFace + 1]; i_node++) { + const mcIdType node = f2n_p[i_node]; + // If this node was modified + if (mapO2N1.find(node) != mapO2N1.end()) { + const auto& new_node = mapO2N1.at(node); + f2n_pw[i_node] = new_node; + noNodesAreChanged = false; + } + } + } + + if (noNodesAreChanged && grpMustBeFullyDup) + throw INTERP_KERNEL::Exception( + "duplicateFaces: A face was supposed to be duplicated but could" + "not be. Please be aware that all M1 groups cannot be" + "duplicated, at least two adjecent inner faces are needed to" + "duplicate a node."); + + lastFace++; + } +} + +CrackAlgo::Set +CrackAlgo::GetCellsTouchingNodesToDup( + const MEDCouplingUMesh * mf, + const DataArrayIdType * n2cIdx, + const DataArrayIdType * n2c, + const DataArrayIdType * f2dup +) { + Set res{}; + + const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex(); + const DataArrayIdType* f2n = mf->getNodalConnectivity(); + + const mcIdType * f2nIdx_p {f2nIdx->begin()}; + const mcIdType * f2n_p {f2n->begin()}; + const mcIdType * n2cIdx_p {n2cIdx->begin()}; + const mcIdType * n2c_p {n2c->begin()}; + + for (const auto face : *f2dup) + for (mcIdType nodeIdx {f2nIdx_p[face] + 1}; nodeIdx < f2nIdx_p[face + 1]; nodeIdx++) { + const auto &node = f2n_p[nodeIdx]; + for (mcIdType cellIdx = n2cIdx_p[node]; cellIdx < n2cIdx_p[node + 1]; cellIdx++) { + const auto &cell = n2c_p[cellIdx]; + res.insert(cell); + } + } + return res; +} + +CrackAlgo::Map2Set CrackAlgo::GetNode2CellMap(const MEDCouplingUMesh * mf, + const DataArrayIdType * n2cIdx, + const DataArrayIdType * n2c, + const DataArrayIdType * f2dup +) { + const DataArrayIdType* f2nIdx = mf->getNodalConnectivityIndex(); + const DataArrayIdType* f2n = mf->getNodalConnectivity(); + + Map2Set n2c_dup{}; + const mcIdType * f2nIdx_p {f2nIdx->begin()}; + const mcIdType * f2n_p {f2n->begin()}; + const mcIdType * n2cIdx_p {n2cIdx->begin()}; + const mcIdType * n2c_p {n2c->begin()}; + + for (const auto face : *f2dup) + for (mcIdType nodeIdx {f2nIdx_p[face] + 1}; nodeIdx < f2nIdx_p[face + 1]; nodeIdx++) { + const auto &node = f2n_p[nodeIdx]; + for (mcIdType cellIdx {n2cIdx_p[node]}; cellIdx < n2cIdx_p[node + 1]; cellIdx++) { + const auto &cell = n2c_p[cellIdx]; + n2c_dup[node].insert(cell); + } + } + return n2c_dup; +} + +void +CrackAlgo::OpenCrack( + MEDFileUMesh * mf, + const Map2Map & cellOld2NewNode, + const double & factor +) { + if ((factor <= 0.0) || (factor >= 1.0)) + throw INTERP_KERNEL::Exception("factor should be between 0.0 and 1.0"); + DataArrayDouble * coords = mf->getCoords(); + const auto dim = static_cast(coords->getNumberOfComponents()); + const double * coords_p = coords->begin(); + MCU m0 = mf->getMeshAtLevel(0); + DAD barys = m0->computeCellCenterOfMass(); + const double * barys_p = barys->begin(); + for (const auto & cellPair : cellOld2NewNode) { + const auto & cell = cellPair.first; + const auto & mapO2N = cellPair.second; + for (const auto & nodeO2N : mapO2N) { + const auto & oldN = nodeO2N.first; + const auto & newN = nodeO2N.second; + double * pos_new_p = coords->rwBegin() + dim * newN; + for (int i = 0; i < dim; i++) { + pos_new_p[i] += (1.0 - factor) * (barys_p[cell*dim+i] - coords_p[dim*oldN+i]); + } + } + } +} + +DataArrayIdType * +CrackAlgo::GetFacesToDupInM1( + const MEDCouplingUMesh & crackMesh, + const MEDCouplingUMesh & m1 +) { + DataArrayIdType * f2dupIdInM1P{}; + const bool allIn = m1.areCellsIncludedIn(&crackMesh, 2, f2dupIdInM1P); + if (!allIn) + throw INTERP_KERNEL::Exception( + "crackAlong: cleanCrackMesh is not part of m1 mesh."); + return f2dupIdInM1P; +} + +pair +CrackAlgo::GetFacesInM1TouchingDuplicatedNodes( + const Map2Set & n2c_dup, + const DataArrayIdType * f2dupIdInM1, + const MEDCouplingUMesh & mf, + const MEDCouplingUMesh & m1 +) { + Set f2dup_set = DataArrayToSet(*f2dupIdInM1); + + DAI n2f(DataArrayIdType::New()); + DAI n2fIdx(DataArrayIdType::New()); + m1.getReverseNodalConnectivity(n2f, n2fIdx); + + // 1. I retrieve all the faces of m1 which have a potentially duplicated + // node, I just use the reverse nodal connectivity of m1 for this. + Set f2changeIdInM1_set{}; + + for (const auto & nodeP : n2c_dup) { + const auto & node = nodeP.first; + for (mcIdType fId = n2fIdx->begin()[node]; fId < n2fIdx->begin()[node+1]; fId++) { + const auto & f = n2f->begin()[fId]; + // f is not a face to duplicate + if (f2dup_set.find(f) == f2dup_set.end()) + f2changeIdInM1_set.insert(f); + } + } + // 2. I retrieve M1 sub mesh + DataArrayIdType * f2changeIdInM1 = SetToDataArray(f2changeIdInM1_set); + MCU part = m1.buildPartOfMySelf(f2changeIdInM1->begin(), f2changeIdInM1->end()); + // 3. I search for the ids in mf + DataArrayIdType * f2changeIdInMf{}; + const bool ok = mf.areCellsIncludedIn(part, 2, f2changeIdInMf); + if (!ok) + throw INTERP_KERNEL::Exception( + "Some faces in -1 mesh next to the crack are not included in " + "the face mesh issued from m0."); + + // 4. I return the 2 DAI + return pair{f2changeIdInM1, f2changeIdInMf}; +} + +void +CrackAlgo::ChangeConnectivityOfM1Elements( + const DataArrayIdType * f2changeIdInM1, + const DataArrayIdType * f2changeIdInMf, + const Map2Map & cellOld2NewNode, + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + MEDCouplingUMesh * m1 +) { + // 1. I go through the faces in question. I go through all the nodes on the + // face. I will look in neighboring cells to see if there have been changes + // in connectivity on the nodes of the face. + // 2. if the changes are consistent, I change in place. + // 3. otherwise I return a consistency error. + + DataArrayIdType* f2nIdx = m1->getNodalConnectivityIndex(); + DataArrayIdType* f2n = m1->getNodalConnectivity(); + + // Change connectivity size + const mcIdType * const f2nIdx_p = f2nIdx->begin(); + const mcIdType * const f2n_p = f2n->begin(); + mcIdType * const f2n_pw = f2n->rwBegin(); + + const mcIdType * const f2cIdx_p {f2cIdx->begin()}; + const mcIdType * const f2c_p {f2c->begin()}; + + const auto f2change_size = static_cast(f2changeIdInM1->getNumberOfTuples()); + + for (size_t iFace = 0; iFace < f2change_size; iFace++) { + const auto & faceM1 = f2changeIdInM1->begin()[iFace]; + const auto & faceMf = f2changeIdInMf->begin()[iFace]; + + // Change inplace the connectivity + // Check that this face is an inner face, ie it is connected to two + // cells + const mcIdType d = f2cIdx_p[faceMf + 1] - f2cIdx_p[faceMf]; + + // 1. In original face, attaching it to cell0 + const auto& cell0 = f2c_p[f2cIdx_p[faceMf]]; + // const auto& cell1 = f2c_p[f2cIdx_p[faceMf] + 1]; + + // If nodes where changed in cell0 + if (cellOld2NewNode.find(cell0) != cellOld2NewNode.end()) { + const auto & mapO2N0 = cellOld2NewNode.at(cell0); + for (mcIdType i_node = f2nIdx_p[faceM1] + 1; i_node < f2nIdx_p[faceM1 + 1]; i_node++) { + const mcIdType node = f2n_p[i_node]; + // If this node was modified + if (mapO2N0.find(node) != mapO2N0.end()) { + const auto& new_node = mapO2N0.at(node); + if ((d == 2) + && ((cellOld2NewNode.find(f2c_p[f2cIdx_p[faceMf] + 1]) + == cellOld2NewNode.end()) + || (cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).find(node) + == cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).end()) + || (cellOld2NewNode.at(f2c_p[f2cIdx_p[faceMf] + 1]).at(node) != new_node))) + throw INTERP_KERNEL::Exception("There is an" + "incoherent change of connectivity between both" + "cell surrounding a face touching the crack but" + "not duplicated."); + f2n_pw[i_node] = new_node; + } + } + } + } +} diff --git a/src/MEDLoader/CrackAlgo.hxx b/src/MEDLoader/CrackAlgo.hxx new file mode 100644 index 000000000..462b771df --- /dev/null +++ b/src/MEDLoader/CrackAlgo.hxx @@ -0,0 +1,193 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : +// webmaster.salome@opencascade.com +// +// Author : Aymeric SONOLET (CEA/DES) + +#ifndef SRC_MEDLOADER_CRACKALGO_HXX_ +#define SRC_MEDLOADER_CRACKALGO_HXX_ + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace MEDCoupling { + +class MEDFileUMesh; +class DataArrayIdType; +class MEDCouplingUMesh; + +class CrackAlgo { + public: + using Set = std::unordered_set; + using Map = std::unordered_map; + using Graph = std::unordered_map; + using Map2Set = std::unordered_map; + using Map2Map = std::unordered_map; + + CrackAlgo() {} + ~CrackAlgo() {} + static Map2Map + Compute(MEDFileUMesh* mm, const std::string& grp_name, bool grpMustBeFullyDup = true); + + static void + OpenCrack( + MEDFileUMesh * mf, + const Map2Map & cellOld2NewNode, + const double & factor = 0.9); + + private: + /* Find connected components using a node to node graph. + */ + static std::vector>> + FindConnectedComponents( + const Graph& graph); + + /* Depth First Search function to find connected components. + */ + static void Dfs( + const Graph& graph, + const mcIdType &node, + const std::size_t &componentId, + std::map& visited, + std::map& componentMap); + + /* Converts DataArrayIdType to Set. + * + * Usefull for doing multiple search in the Set efficiently. + */ + static Set + DataArrayToSet( + const DataArrayIdType &da); + + static DataArrayIdType * + SetToDataArray( + const Set& s); + + static Set + GetCellsTouchingNodesToDup( + const MEDCouplingUMesh * mf, + const DataArrayIdType * n2cIdx, + const DataArrayIdType * n2c, + const DataArrayIdType * f2dup); + + static Map2Set + GetNode2CellMap( + const MEDCouplingUMesh * mf, + const DataArrayIdType * n2cIdx, + const DataArrayIdType * n2c, + const DataArrayIdType * f2dup); + + /* Building the cell to cell graph. + * + * This graph concerns only cells touching nodes touching faces to duplicate. + * The connection between cells sharing a face to dup is cut in this graph. + */ + static Graph + BuildCutC2CGraph( + const DataArrayIdType * c2fIdx, + const DataArrayIdType * c2f, + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + const Set & cTouchingN_dup, + const DataArrayIdType * f2dup); + + static Map2Map + CreateNewNodesInTopLevelMesh( + const Map2Set & n2c_dup, + const Graph & c2c, + MEDCouplingUMesh * m0); + + static void + AddMissingElementsOnLevelM1AndChangeConnectivity( + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + const DataArrayIdType * f2dupIdInM1, + const DataArrayIdType * f2dupIdInMf, + const Map2Map & cellOld2NewNode, + MEDCouplingUMesh * m1, + bool grpMustBeFullyDup = true); + + /* Create new family array for level -1, which is the extended copy of the + * original as new elements are appended. Caller owns newFamily DAI. + * Supposes that the elements to duplicate are all duplicated and appended + * in the order at the end of mf. + */ + static DataArrayIdType * + CopyAndCompleteFamilyArrAtLevelM1( + const MEDFileUMesh * mm, + const MEDCouplingUMesh * mf, + const DataArrayIdType * f2dup); + + static DataArrayIdType * + CopyFamilyArrAtLev0(const MEDFileUMesh * mm); + + /* Manage node families. + * + * Extend the family size inplace and set new nodes family to their + * predecessor. + */ + static void + CompleteFamilyArrAtNodeLevel( + const Map2Set& addedNodes, + MEDFileUMesh * mm); + + /* Aggregate CellOld2NewNodes into oldNode to newNodes. + */ + static Map2Set + BuildMap2Set( + const Map2Map& cellOld2NewNode); + + /* Remomves cells from crackingMesh which are part of m skin. + * + * Returns a new MEDCouplingUMesh. User is responsible to delete it. + */ + static MEDCouplingUMesh * + CleanM1Mesh( + const MEDCouplingUMesh & m, + const MEDCouplingUMesh & crackingMesh); + + static std::pair + GetFacesInM1TouchingDuplicatedNodes( + const Map2Set & n2c_dup, + const DataArrayIdType * f2dupIdInM1, + const MEDCouplingUMesh & mf, + const MEDCouplingUMesh & m1); + + static DataArrayIdType * + GetFacesToDupInM1( + const MEDCouplingUMesh & crackMesh, + const MEDCouplingUMesh & m1); + + static void + ChangeConnectivityOfM1Elements( + const DataArrayIdType * f2changeIdInM1, + const DataArrayIdType * f2changeIdInMf, + const Map2Map & cellOld2New, + const DataArrayIdType * f2cIdx, + const DataArrayIdType * f2c, + MEDCouplingUMesh * m1); +}; +} // namespace MEDCoupling +#endif // SRC_MEDLOADER_CRACKALGO_HXX_ diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index 7a0bd5f97..59791f7cf 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -25,6 +25,7 @@ #include "MEDLoaderNS.hxx" #include "MEDFileSafeCaller.txx" #include "MEDLoaderBase.hxx" +#include "CrackAlgo.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCouplingMappedExtrudedMesh.hxx" @@ -32,8 +33,13 @@ #include "InterpKernelAutoPtr.hxx" +#include #include #include +#include +#include +#include +#include // From MEDLOader.cxx TU extern med_geometry_type typmai[MED_N_CELL_FIXED_GEO]; @@ -4201,6 +4207,51 @@ void MEDFileUMesh::optimizeFamilies() _groups.erase(*it); } +/** + * \b this must be filled at level 0 and -1, the -1 level being part + * of the descending connectivity of the top level. This method build a + * "crack", or an inner boundary, in \b this along the group of level -1 named + * grpNameM1. The boundary is built according to the following method: + * - all nodes which are not connected anymore to its neighboring cells due to + * the crack are duplicated/tripled/... as needed. + * - new (-1)-level cells are built lying on those new nodes. So the + * edges/faces along the group are duplicated. The duplicated + * (-1)-level cells are added to the grpNameM1. + * - connectivity modification at lev -1 according to the new connectivity at + * lev 0. + * Note that some M1 group might lead to no some cells at level 0 being still + * connected despite sharing only a face in M1. The behavior of the method + * depends on the grpMustBeFullyDup boolean input parameter in this case. + * Groups at lev -1 are managed by adding the new duplicated faces to the same + * groups of their twin. + * Connectivity of mesh at lev -2 is not managed, so it could be not valid + * anymore. + * After this operation: + * - a top-level cell bordering the group will loose some + * neighbors (typically the cell which is on the other side of the group is no + * more a neighbor) + * - the connectivity of (part of) the top level-cells bordering the group is + * also modified so that some cells bordering the newly created boundary use + * the newly computed nodes. + * + * \param[in] grpNameM1 name of the (-1)-level group defining the boundary + * \param[in] grpMustBeFullyDup boolean indicating if the method should throw + * an exeption in the case the M1 group leads to some cells at lev 0 not being + * seperated. If set to false, duplicated faces are created but are identical. + * They share the same connectivity. + * \param[out] cell2old2newNode a dict of dict indicating for each modified + * cells at lev 0 the new node connectivity. + */ +std::unordered_map> +MEDFileUMesh::crackAlong(const std::string &grp_name, bool grpMustBeFullyDup) { + return CrackAlgo::Compute(this, grp_name, grpMustBeFullyDup); +} + +void MEDFileUMesh::openCrack(const std::unordered_map> & c2o2nN, const double & factor) +{ + CrackAlgo::OpenCrack(this, c2o2nN, factor); +} + /** * \b this must be filled at level 0 and -1, typically the -1 level being (part of) the descending connectivity * of the top level. This method build a "crack", or an inner boundary, in \b this along the group of level -1 named grpNameM1. diff --git a/src/MEDLoader/MEDFileMesh.hxx b/src/MEDLoader/MEDFileMesh.hxx index 91b377d82..4c2adbc4c 100644 --- a/src/MEDLoader/MEDFileMesh.hxx +++ b/src/MEDLoader/MEDFileMesh.hxx @@ -31,6 +31,9 @@ #include #include +#include +#include +#include namespace MEDCoupling { @@ -363,6 +366,8 @@ MCAuto& coords, MCAuto& partCoords, MCAuto& ms, bool renum=false); MEDLOADER_EXPORT void optimizeFamilies(); // tools + MEDLOADER_EXPORT std::unordered_map> crackAlong(const std::string &grpNameM1, bool grpMustBeFullyDup=true); + MEDLOADER_EXPORT void openCrack(const std::unordered_map> & c2o2nN, const double & factor); MEDLOADER_EXPORT void buildInnerBoundaryAlongM1Group(const std::string& grpNameM1, DataArrayIdType *&nodesDuplicated, DataArrayIdType *&cellsModified, DataArrayIdType *&cellsNotModified); MEDLOADER_EXPORT bool unPolyze(std::vector& oldCode, std::vector& newCode, DataArrayIdType *& o2nRenumCell); MEDLOADER_EXPORT DataArrayIdType *zipCoords(); @@ -398,6 +403,17 @@ MCAuto& coords, MCAuto& partCoords, MCAuto > getAllNonNullFamilyIds() const; MCAuto& checkAndGiveEntryInSplitL1(int meshDimRelToMax, MEDCouplingPointSet *m); + static std::vector>> + findConnectedComponents( + const std::map> &graph + ); + static void dfs( + const std::map> &graph, + const mcIdType &node, + const std::size_t& componentId, + std::map& visited, + std::map& componentMap + ); private: static const char SPE_FAM_STR_EXTRUDED_MESH[]; diff --git a/src/MEDLoader/Swig/CMakeLists.txt b/src/MEDLoader/Swig/CMakeLists.txt index 2cc09ced6..479afcf42 100644 --- a/src/MEDLoader/Swig/CMakeLists.txt +++ b/src/MEDLoader/Swig/CMakeLists.txt @@ -99,7 +99,7 @@ SALOME_INSTALL_SCRIPTS(${CMAKE_CURRENT_BINARY_DIR}/MEDLoader.py ${MEDCOUPLING_IN INSTALL(FILES MEDLoaderSplitter.py MEDLoaderFinalize.py DESTINATION ${MEDCOUPLING_INSTALL_PYTHON}) -INSTALL(FILES ${ALL_TESTS} MEDLoaderDataForTest.py MEDLoaderTest1.py MEDLoaderTest2.py MEDLoaderTest3.py DESTINATION ${MEDCOUPLING_INSTALL_SCRIPT_SCRIPTS}) +INSTALL(FILES ${ALL_TESTS} MEDLoaderDataForTest.py MEDLoaderTest1.py MEDLoaderTest2.py MEDLoaderTest3.py CrackAlongTest.py DESTINATION ${MEDCOUPLING_INSTALL_SCRIPT_SCRIPTS}) INSTALL(FILES CaseIO.py CaseReader.py CaseWriter.py VTKReader.py DESTINATION ${MEDCOUPLING_INSTALL_PYTHON}) @@ -138,6 +138,7 @@ SET(MEDLOADER_TEST_FILES MEDLoaderTest1.py MEDLoaderTest2.py MEDLoaderTest3.py + CrackAlongTest.py CaseIO.py CaseReader.py CaseWriter.py diff --git a/src/MEDLoader/Swig/CrackAlongTest.py b/src/MEDLoader/Swig/CrackAlongTest.py new file mode 100644 index 000000000..d93ffae73 --- /dev/null +++ b/src/MEDLoader/Swig/CrackAlongTest.py @@ -0,0 +1,648 @@ +# -*- coding: iso-8859-1 -*- +# Copyright (C) 2007-2024 CEA, EDF +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author : Aymeric SONOLET (CEA) + +import unittest + +from medcoupling import ( + DataArrayDouble, + DataArrayInt, + MEDCouplingCMesh, + MEDCouplingUMesh, +) +from MEDLoader import MEDFileUMesh + +from MEDLoaderDataForTest import WriteInTmpDir + + +class CrackAlongTest(unittest.TestCase): + @WriteInTmpDir + def testBuildInnerBoundaryAlongM1Group1(self): + fname = "Pyfile44.med" + m = MEDCouplingCMesh.New() + m.setCoordsAt(0, DataArrayDouble.New([0.0, 1.1, 2.3, 3.6, 5.0, 6.5])) + m.setCoordsAt(1, DataArrayDouble.New([0.0, 1.1, 2.3, 3.6, 5.0])) + m = m.buildUnstructured() + m.setName("AnthonyDuplicate") + m.getCoords().setInfoOnComponents(["X [km]", "Z [mm]"]) + m2 = m.buildDescendingConnectivity()[0][ + [8, 11, 14, 20, 21, 22, 23, 24, 25, 26, 31, 32, 33, 34, 35, 36, 37] + ] + m2.setName(m.getName()) + grp = DataArrayInt.New([4, 6, 8]) + grp.setName("Grp") + grp2 = DataArrayInt.New([9, 16]) + grp2.setName("Grp2") + mm = MEDFileUMesh.New() + mm.setMeshAtLevel(0, m) + mm.setMeshAtLevel(-1, m2) + mm.setGroupsAtLevel(-1, [grp, grp2]) + grpNode = DataArrayInt.New([4, 21, 23]) + grpNode.setName("GrpNode") + mm.setGroupsAtLevel(1, [grpNode]) + ref0 = [4, 15, 14, 20, 21, 4, 16, 15, 21, 22, 4, 17, 16, 22, 23] + ref1 = [4, 9, 8, 14, 15, 4, 10, 9, 15, 16, 4, 11, 10, 16, 17] + ref2 = [4, 32, 14, 20, 21, 4, 31, 32, 21, 22, 4, 30, 31, 22, 23] + # + self.assertEqual(30, mm.getNumberOfNodes()) + self.assertEqual( + ref0, mm.getMeshAtLevel(0)[[12, 13, 14]].getNodalConnectivity().getValues() + ) + self.assertEqual( + ref1, mm.getMeshAtLevel(0)[[7, 8, 9]].getNodalConnectivity().getValues() + ) + # + c2o2nN = mm.crackAlong("Grp") + self.assertEqual( + {12: {15: 32}, 13: {15: 32, 16: 31}, 14: {16: 31, 17: 30}}, c2o2nN + ) + self.assertEqual(33, mm.getNumberOfNodes()) + self.assertEqual([4, 6, 8, 17, 18, 19], mm.getGroupArr(-1, "Grp").getValues()) + self.assertEqual([9, 16], mm.getGroupArr(-1, "Grp2").getValues()) + self.assertEqual([4, 21, 23], mm.getGroupArr(1, "GrpNode").getValues()) + self.assertEqual( + ref2, mm.getMeshAtLevel(0)[[12, 13, 14]].getNodalConnectivity().getValues() + ) # cells 7,8,9 and 12,13,14 are lying on "Grp" but only 12,13,14 are renumbered + self.assertEqual( + ref1, mm.getMeshAtLevel(0)[[7, 8, 9]].getNodalConnectivity().getValues() + ) # + # fmt: off + refValues = DataArrayDouble.New( + [1.21, 1.32, 1.43, 1.54, 1.65, 1.32, 1.44, 1.56, 1.68, 1.80, + 1.43, 1.56, 1.69, 1.82, 1.95, 1.54, 1.68, 1.82, 1.96, 2.10] + ) + # fmt: on + valsToTest = mm.getMeshAtLevel(0).getMeasureField(True).getArray() + delta = valsToTest - refValues + delta.abs() + self.assertTrue(delta.getMaxValue()[0] < 1e-12) + + @WriteInTmpDir + def testBuildInnerBoundaryAlongM1Group2(self): + fname = "Pyfile45.med" + m = MEDCouplingCMesh.New() + m.setCoordsAt(0, DataArrayDouble.New([0.0, 1.1, 2.3, 3.6, 5.0, 6.5])) + m.setCoordsAt(1, DataArrayDouble.New([0.0, 1.1, 2.3, 3.6, 5.0])) + m = m.buildUnstructured() + m.setName("AnthonyDuplicate") + m.getCoords().setInfoOnComponents(["X [km]", "Z [mm]"]) + m2 = m.buildDescendingConnectivity()[0][ + [8, 11, 14, 20, 21, 22, 23, 24, 25, 26, 31, 32, 33, 34, 35, 36, 37] + ] + m2.setName(m.getName()) + grp = DataArrayInt.New([4, 6]) + grp.setName("Grp") + grp2 = DataArrayInt.New([9, 16]) + grp2.setName("Grp2") + mm = MEDFileUMesh.New() + mm.setMeshAtLevel(0, m) + mm.setMeshAtLevel(-1, m2) + mm.setGroupsAtLevel(-1, [grp, grp2]) + grpNode = DataArrayInt.New([4, 21, 23]) + grpNode.setName("GrpNode") + mm.setGroupsAtLevel(1, [grpNode]) + ref0 = [4, 15, 14, 20, 21, 4, 16, 15, 21, 22, 4, 17, 16, 22, 23] + ref1 = [4, 9, 8, 14, 15, 4, 10, 9, 15, 16] + ref2 = [4, 30, 14, 20, 21, 4, 16, 30, 21, 22, 4, 17, 16, 22, 23] + # + self.assertEqual(30, mm.getNumberOfNodes()) + self.assertEqual( + ref0, mm.getMeshAtLevel(0)[[12, 13, 14]].getNodalConnectivity().getValues() + ) + self.assertEqual( + ref1, mm.getMeshAtLevel(0)[[7, 8]].getNodalConnectivity().getValues() + ) + # + c2o2nN = mm.crackAlong("Grp") + self.assertEqual({13: {15: 30}, 12: {15: 30}}, c2o2nN) + self.assertEqual(31, mm.getNumberOfNodes()) + self.assertEqual([4, 6, 17, 18], mm.getGroupArr(-1, "Grp").getValues()) + self.assertEqual([9, 16], mm.getGroupArr(-1, "Grp2").getValues()) + self.assertEqual([4, 21, 23], mm.getGroupArr(1, "GrpNode").getValues()) + self.assertEqual( + ref2, mm.getMeshAtLevel(0)[[12, 13, 14]].getNodalConnectivity().getValues() + ) # cells 7,8,9 and 12,13,14 are lying on "Grp" but only 12 and 13 are renumbered + self.assertEqual( + ref1, mm.getMeshAtLevel(0)[[7, 8]].getNodalConnectivity().getValues() + ) + # fmt: off + refValues = DataArrayDouble.New( + [1.21, 1.32, 1.43, 1.54, 1.65, 1.32, 1.44, 1.56, 1.68, 1.80, + 1.43, 1.56, 1.69, 1.82, 1.95, 1.54, 1.68, 1.82, 1.96, 2.10] + ) + # fmt: on + valsToTest = mm.getMeshAtLevel(0).getMeasureField(True).getArray() + delta = valsToTest - refValues + delta.abs() + self.assertTrue(delta.getMaxValue()[0] < 1e-12) + + @WriteInTmpDir + def testBuildInnerBoundaryAlongM1Group3(self): + """Test buildInnerBoundaryAlongM1Group() with *non-connex* cracks""" + fname = "Pyfile73.med" + m = MEDCouplingCMesh.New() + m.setCoordsAt(0, DataArrayDouble([0.0, 1.1, 2.3, 3.6, 5.0])) + m.setCoordsAt(1, DataArrayDouble([0.0, 1.0, 2.0])) + m = m.buildUnstructured() + m.setName("simple") + m2 = m.buildDescendingConnectivity()[0] + m2.setName(m.getName()) + + # A crack in two non connected parts of the mesh: + grpSeg = DataArrayInt([3, 19]) + grpSeg.setName("Grp") + + mm = MEDFileUMesh.New() + mm.setMeshAtLevel(0, m) + mm.setMeshAtLevel(-1, m2) + mm.setGroupsAtLevel(-1, [grpSeg]) + c2o2nN = mm.crackAlong("Grp") + self.assertEqual({1: {1: 16}, 7: {13: 15}}, c2o2nN) + # self.assertEqual([1, 13], nodes.getValues()) + # self.assertEqual([0, 6], cellsMod.getValues()) + # self.assertEqual([1, 7], cellsNotMod.getValues()) + self.assertEqual(17, mm.getNumberOfNodes()) + self.assertEqual({3, 19, 22, 23}, set(mm.getGroupArr(-1, "Grp").getValues())) + + refValues = DataArrayDouble([1.1, 1.2, 1.3, 1.4, 1.1, 1.2, 1.3, 1.4]) + valsToTest = mm.getMeshAtLevel(0).getMeasureField(True).getArray() + delta = valsToTest - refValues + delta.abs() + self.assertTrue(delta.getMaxValue()[0] < 1e-10) + + @WriteInTmpDir + def testBuildInnerBoundaryAlongM1Group4(self): + """Test case where cells touch the M1 group on some nodes only and not + on full edges (triangle mesh for ex) + """ + # fmt: off + coo = DataArrayDouble( + [0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 0.0, 1.0, 1.0, 1.0, + 2.0, 1.0, 3.0, 1.0, 0.0, 2.0, 1.0, 2.0, 2.0, 2.0, 3.0, 2.0], + 12, + 2, + ) + conn = [3, 0, 4, 1, + 3, 1, 4, 5, + 3, 5, 9, 10, + 3, 5, 10, 6, + 3, 2, 6, 7, + 3, 2, 7, 3, + 3, 4, 8, 9, + 3, 4, 9, 5, + 3, 1, 5, 6, + 3, 1, 6, 2, + 3, 6, 10, 11, + 3, 6, 11, 7] + # fmt: on + # Only TRI3: + connI = DataArrayInt() + connI.alloc(13, 1) + connI.iota() + connI *= 4 + m2 = MEDCouplingUMesh("2D", 2) + m2.setCoords(coo) + m2.setConnectivity(DataArrayInt(conn), connI) + m2.checkConsistency() + m1, _, _, _, _ = m2.buildDescendingConnectivity() + grpIds = DataArrayInt([9, 11]) + grpIds.setName("group") + grpIds2 = DataArrayInt([0, 1]) + grpIds2.setName("group2") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m2) + mfu.setMeshAtLevel(-1, m1) + mfu.setGroupsAtLevel(-1, [grpIds, grpIds2]) + nNod = m2.getNumberOfNodes() + + ref0 = [3, 5, 10, 6, 3, 2, 7, 3, 3, 6, 10, 11] + ref1 = [3, 2, 6, 7, 3, 1, 5, 6, 3, 1, 6, 2, 3, 6, 11, 7] + ref2 = [3, 2, 13, 7, 3, 1, 5, 13, 3, 1, 13, 2, 3, 6, 11, 12] + + self.assertEqual(ref0, m2[[3, 5, 10]].getNodalConnectivity().getValues()) + self.assertEqual(ref1, m2[[4, 8, 9, 11]].getNodalConnectivity().getValues()) + + c2o2nN = mfu.crackAlong("group") + self.assertEqual({9: {6: 13}, 8: {6: 13}, 4: {6: 13}, 11: {7: 12}}, c2o2nN) + + self.assertEqual(ref0, m2[[3, 5, 10]].getNodalConnectivity().getValues()) + self.assertEqual(ref2, m2[[4, 8, 9, 11]].getNodalConnectivity().getValues()) + + m2_bis = mfu.getMeshAtLevel(0) + m2_bis.checkConsistency() + m1_bis = mfu.getMeshAtLevel(-1) + m1_bis.checkConsistency() + self.assertEqual(nNod + 2, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 2, m2_bis.getNumberOfNodes()) + self.assertEqual(nNod + 2, m1_bis.getNumberOfNodes()) + self.assertEqual([9, 11, 23, 24], mfu.getGroupArr(-1, "group").getValues()) + self.assertEqual([0, 1], mfu.getGroupArr(-1, "group2").getValues()) + + m_bis0 = mfu.getMeshAtLevel(-1) + m_desc, _, _, _, _ = m_bis0.buildDescendingConnectivity() + m_bis0.checkDeepEquivalOnSameNodesWith(mfu.getMeshAtLevel(-1), 2, 9.9999999) + + @WriteInTmpDir + def testBuildInnerBoundary5(self): + """Full 3D test with tetras only. In this case a tri from the group is not duplicated because it is made only + of non duplicated nodes. The tri in question is hence not part of the final new "dup" group.""" + # fmt: off + coo = DataArrayDouble( + [200.0, 200.0, 0.0, 200.0, 200.0, 200.0, 200.0, 0.0, 200.0, 200.0, 0.0, 0.0, 0.0, 200.0, 0.0, 0.0, 200.0, 200.0, 0.0, 0.0, 0.0, 0.0, 0.0, 200.0, 400.0, 200.0, 0.0, 400.0, 200.0, 200.0, 400.0, 0.0, 0.0, 400.0, 0.0, 200.0, 0.0, 100.00000000000016, 200.0, 63.15203310314546, 200.0, 200.0, 134.45205700643342, 200.0, 200.0, 200.0, 100.00000000000016, 200.0, 63.15203310314546, 0.0, 200.0, 134.45205700643342, 0.0, 200.0, 0.0, 100.00000000000016, 0.0, 63.15203310314546, 200.0, 0.0, 134.45205700643342, 200.0, 0.0, 200.0, 100.00000000000016, 0.0, 63.15203310314546, 0.0, 0.0, 134.45205700643342, 0.0, 0.0, 200.0, 200.0, 100.02130053568538, 0.0, 200.0, 100.00938163175135, 200.0, 0.0, 100.02130053568538, 0.0, 0.0, 100.00938163175135, 299.3058739933347, 200.0, 200.0, 400.0, 98.68100542924483, 200.0, 302.8923433403344, 0.0, 200.0, 302.8923433403344, 200.0, 0.0, 400.0, 100.00000000000016, 0.0, 302.8923433403344, 0.0, 0.0, 400.0, 200.0, 98.55126825835082, 400.0, 0.0, 100.02162286181577, 99.31624553977466, 99.99999998882231, 200.0, 99.31624576683302, 100.00000010178034, 0.0, 99.31624560596512, 200.0, 100.0050761312483, 99.31624560612883, 0.0, 100.00507613125338, 200.0, 99.99999995813045, 100.00950673487786, 0.0, 99.99999989928207, 100.0041870621175, 301.29063354383015, 100.0000000093269, 0.0, 301.29063360689975, 0.0, 100.00957769061164, 140.52853868782435, 99.99999963972768, 100.00509135751312, 297.87779091770784, 97.16750463405486, 97.18018457127863], + 46, + 3, + ) + c0 = [14, 45, 31, 21, 42, 14, 37, 38, 20, 44, 14, 39, 36, 41, 44, 14, 5, 25, 12, 13, 14, 38, 36, 44, 41, 14, 21, 20, 24, 44, 14, 38, 25, 41, 19, 14, 37, 38, 44, 41, 14, 16, 27, 39, 41, 14, 21, 45, 26, 40, 14, 39, 37, 44, 41, 14, 14, 15, 24, 44, 14, 25, 38, 41, 13, 14, 27, 18, 6, 22, 14, 38, 36, 41, 13, 14, 44, 14, 15, 36, 14, 44, 23, 39, 26, 14, 21, 26, 23, 44, 14, 38, 44, 14, 24, 14, 39, 37, 41, 22, 14, 21, 33, 45, 42, 14, 27, 22, 39, 41, 14, 23, 26, 21, 3, 14, 27, 18, 22, 41, 14, 39, 36, 44, 17, 14, 21, 26, 44, 40, 14, 39, 37, 22, 23, 14, 37, 38, 41, 19, 14, 25, 12, 13, 41, 14, 30, 26, 43, 45, 14, 38, 36, 13, 14, 14, 12, 36, 13, 41, 14, 20, 44, 21, 37, 14, 16, 36, 12, 41, 14, 39, 36, 17, 16, 14, 44, 20, 24, 38, 14, 27, 16, 12, 41, 14, 26, 15, 17, 44, 14, 19, 18, 41, 37, 14, 40, 45, 26, 15, 14, 37, 38, 19, 20, 14, 17, 15, 26, 2, 14, 39, 36, 16, 41, 14, 24, 21, 44, 40, 14, 16, 7, 27, 12, 14, 22, 18, 37, 41, 14, 21, 31, 45, 24, 14, 44, 40, 15, 24, 14, 24, 45, 15, 28, 14, 44, 40, 26, 15, 14, 24, 20, 21, 0, 14, 38, 36, 14, 44, 14, 39, 37, 23, 44, 14, 45, 31, 42, 32, 14, 25, 18, 19, 4, 14, 36, 44, 17, 15, 14, 25, 19, 18, 41, 14, 24, 15, 14, 1, 14, 45, 24, 34, 28, 14, 35, 45, 30, 43, 14, 17, 44, 39, 26, 14, 44, 23, 21, 37, 14, 30, 45, 29, 15, 14, 45, 35, 33, 43, 14, 30, 15, 26, 45, 14, 31, 21, 0, 24, 14, 33, 35, 32, 10, 14, 29, 45, 34, 28, 14, 32, 45, 34, 29, 14, 45, 31, 32, 34, 14, 33, 26, 45, 43, 14, 45, 31, 34, 24, 14, 33, 26, 21, 45, 14, 11, 30, 35, 29, 14, 33, 35, 45, 32, 14, 33, 45, 42, 32, 14, 32, 8, 34, 31, 14, 21, 26, 33, 3, 14, 35, 45, 32, 29, 14, 29, 34, 9, 28, 14, 15, 45, 24, 40, 14, 29, 45, 28, 15, 14, 21, 24, 45, 40, 14, 24, 15, 1, 28, 14, 35, 45, 29, 30, 14, 26, 15, 30, 2] + cI0 = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170, 175, 180, 185, 190, 195, 200, 205, 210, 215, 220, 225, 230, 235, 240, 245, 250, 255, 260, 265, 270, 275, 280, 285, 290, 295, 300, 305, 310, 315, 320, 325, 330, 335, 340, 345, 350, 355, 360, 365, 370, 375, 380, 385, 390, 395, 400, 405, 410, 415, 420, 425, 430] + # fmt: on + m3 = MEDCouplingUMesh("3D", 3) + m3.setCoords(coo) + m3.setConnectivity(DataArrayInt(c0), DataArrayInt(cI0)) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([36, 74]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + grpIds3D = DataArrayInt([0, 1]) + grpIds3D.setName("group_3d") + mfu.setGroupsAtLevel(0, [grpIds3D]) # just to check preservation of 3D group + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m3.getNumberOfNodes() + + c2o2nN = mfu.crackAlong("group", grpMustBeFullyDup=False) + self.assertEqual({77: {3: 46}}, c2o2nN) + + m3_bis = mfu.getMeshAtLevel(0) + m3_bis.checkConsistency() + m2_bis = mfu.getMeshAtLevel(-1) + m2_bis.checkConsistency() + self.assertEqual(nNod + 1, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 1, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod + 1, m2_bis.getNumberOfNodes()) + self.assertEqual( + m3_bis.getCoords()[3].getValues(), m3_bis.getCoords()[nNod:].getValues() + ) + self.assertEqual([36, 74, 213, 214], mfu.getGroupArr(-1, "group").getValues()) + self.assertEqual([0, 1], mfu.getGroupArr(0, "group_3d").getValues()) + m_bis0 = mfu.getMeshAtLevel(-1) + m_desc, _, _, _, _ = m_bis0.buildDescendingConnectivity() + m_bis0.checkDeepEquivalOnSameNodesWith(mfu.getMeshAtLevel(-1), 2, 9.9999999) + + @WriteInTmpDir + def testBuildInnerBoundary6(self): + """3D test where the crack has a funny shape with a singular point (i.e. two faces of the M1 group are only connected by one point, not a full segment) + The singular point was wrongly duplicated. + """ + # fmt: off + coo = DataArrayDouble( + [(-1.38778e-17, 0.226, 0), (-1.38778e-17, -1.38778e-17, 0), (0.226, 0.226, 0), (0.226, -1.38778e-17, 0), (0.452, 0.226, 0), (0.452, -1.38778e-17, 0), (-1.38778e-17, 0.452, 0), (0.226, 0.452, 0), (0.452, 0.452, 0), (-1.38778e-17, 0.226, 0.25), (0.226, 0.226, 0.25), (0.226, -1.38778e-17, 0.25), (-1.38778e-17, -1.38778e-17, 0.25), (-1.38778e-17, 0.226, 0.779375), (0.226, 0.226, 0.779375), (0.226, -1.38778e-17, 0.779375), (-1.38778e-17, -1.38778e-17, 0.779375), (-1.38778e-17, 0.226, 1.30875), (0.226, 0.226, 1.30875), (0.226, -1.38778e-17, 1.30875), (-1.38778e-17, -1.38778e-17, 1.30875), (0.452, 0.226, 0.25), (0.452, -1.38778e-17, 0.25), (0.452, 0.226, 0.779375), (0.452, -1.38778e-17, 0.779375), (0.452, 0.226, 1.30875), (0.452, -1.38778e-17, 1.30875), (-1.38778e-17, 0.452, 0.25), (0.226, 0.452, 0.25), (-1.38778e-17, 0.452, 0.779375), (0.226, 0.452, 0.779375), (-1.38778e-17, 0.452, 1.30875), (0.226, 0.452, 1.30875), (0.452, 0.452, 0.25), (0.452, 0.452, 0.779375), (0.452, 0.452, 1.30875), (0.146, 0.226, 0.779375), (0.146, -1.38778e-17, 0.779375), (0.146, 0.226, 1.30875), (0.146, -1.38778e-17, 1.30875), (0.146, 0.452, 0.779375), (0.146, 0.452, 1.30875)] + ) + c0 = [18, 0, 2, 3, 1, 9, 10, 11, 12, 18, 9, 10, 11, 12, 13, 36, 37, 16, 18, 13, 36, 37, 16, 17, 38, 39, 20, 18, 2, 4, 5, 3, 10, 21, 22, 11, 18, 10, 21, 22, 11, 14, 23, 24, 15, 18, 14, 23, 24, 15, 18, 25, 26, 19, 18, 6, 7, 2, 0, 27, 28, 10, 9, 18, 27, 28, 10, 9, 29, 40, 36, 13, 18, 29, 40, 36, 13, 31, 41, 38, 17, 18, 7, 8, 4, 2, 28, 33, 21, 10, 18, 28, 33, 21, 10, 30, 34, 23, 14, 18, 30, 34, 23, 14, 32, 35, 25, 18] + # fmt: on + cI0 = [0, 9, 18, 27, 36, 45, 54, 63, 72, 81, 90, 99, 108] + m3 = MEDCouplingUMesh("3D", 3) + m3.setCoords(coo) + m3.setConnectivity(DataArrayInt(c0), DataArrayInt(cI0)) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([7, 12, 22, 27]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m3.getNumberOfNodes() + c2o2nN = mfu.crackAlong("group") + self.assertEqual( + { + 7: {13: 49, 36: 48}, + 8: {13: 49, 36: 48, 17: 46, 38: 45}, + 10: {23: 47, 14: 43}, + 11: {23: 47, 25: 44, 14: 43, 18: 42}, + }, + c2o2nN, + ) + m3_bis = mfu.getMeshAtLevel(0) + m3_bis.checkConsistency() + m2_bis = mfu.getMeshAtLevel(-1) + m2_bis.checkConsistency() + self.assertEqual(nNod + 8, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 8, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod + 8, m2_bis.getNumberOfNodes()) + # self.assertEqual([13, 14, 17, 18, 23, 25, 36, 38], nodesDup.getValues()) + # self.assertEqual( + # m3_bis.getCoords()[nodesDup].getValues(), + # m3_bis.getCoords()[nNod:].getValues(), + # ) + # self.assertEqual(set([1, 2, 4, 5]), set(cells1.getValues())) + # self.assertEqual(set([7, 8, 10, 11]), set(cells2.getValues())) + self.assertEqual( + [7, 12, 22, 27, 56, 57, 58, 59], mfu.getGroupArr(-1, "group").getValues() + ) + m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity() + m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) + + @WriteInTmpDir + def testBuildInnerBoundary7(self): + """3D test where the crack has another funny shape with another singular point (i.e. two faces of the M1 group are only connected by one point, not a full segment) + Once the crack is inserted, the cells on either side of the crack do not necessarily form a connex spread zone. This was not properly handled either. + """ + # fmt: off + coo = DataArrayDouble( + [(5, 17, 0), (0, 17, 0), (0, 12, 0), (5, 12, 0), (15, 17, 0), (15, 12, 0), (20, 12, 0), (20, 17, 0), (20, 2, 0), (15, 2, 0), (15, -3, 0), (20, -3, 0), (5, -3, 0), (5, 2, 0), (0, -3, 0), (0, 2, 0), (5, 17, 10), (5, 17, 20), (5, 17, 30), (5, 17, 40), (0, 17, 10), (0, 17, 20), (0, 17, 30), (0, 17, 40), (0, 12, 10), (0, 12, 20), (0, 12, 30), (0, 12, 40), (5, 12, 10), (5, 12, 20), (5, 12, 30), (5, 12, 40), (15, 17, 10), (15, 17, 20), (15, 17, 30), (15, 17, 40), (15, 12, 10), (15, 12, 20), (15, 12, 30), (15, 12, 40), (20, 12, 10), (20, 12, 20), (20, 12, 30), (20, 12, 40), (20, 17, 10), (20, 17, 20), (20, 17, 30), (20, 17, 40), (20, 2, 10), (20, 2, 20), (20, 2, 30), (20, 2, 40), (15, 2, 10), (15, 2, 20), (15, 2, 30), (15, 2, 40), (15, -3, 10), (15, -3, 20), (15, -3, 30), (15, -3, 40), (20, -3, 10), (20, -3, 20), (20, -3, 30), (20, -3, 40), (5, -3, 10), (5, -3, 20), (5, -3, 30), (5, -3, 40), (5, 2, 10), (5, 2, 20), (5, 2, 30), (5, 2, 40), (0, -3, 10), (0, -3, 20), (0, -3, 30), (0, -3, 40), (0, 2, 10), (0, 2, 20), (0, 2, 30), (0, 2, 40), (20, 8, 0), (0, 8, 0), (20, 8, 10), (20, 8, 20), (20, 8, 30), (20, 8, 40), (15, 8, 30), (15, 8, 40), (5, 8, 30), (5, 8, 40), (0, 8, 10), (0, 8, 20), (0, 8, 30), (0, 8, 40)] + ) + c = DataArrayInt( + [31, 0, 3, 2, 1, -1, 16, 20, 24, 28, -1, 0, 16, 28, 3, -1, 3, 28, 24, 2, -1, 2, 24, 20, 1, -1, 1, 20, 16, 0, 31, 16, 28, 24, 20, -1, 17, 21, 25, 29, -1, 16, 17, 29, 28, -1, 28, 29, 25, 24, -1, 24, 25, 21, 20, -1, 20, 21, 17, 16, 31, 17, 29, 25, 21, -1, 18, 22, 26, 30, -1, 17, 18, 30, 29, -1, 29, 30, 26, 25, -1, 25, 26, 22, 21, -1, 21, 22, 18, 17, 31, 18, 30, 26, 22, -1, 19, 23, 27, 31, -1, 18, 19, 31, 30, -1, 30, 31, 27, 26, -1, 26, 27, 23, 22, -1, 22, 23, 19, 18, 31, 4, 5, 3, 0, -1, 32, 16, 28, 36, -1, 4, 32, 36, 5, -1, 5, 36, 28, 3, -1, 3, 28, 16, 0, -1, 0, 16, 32, 4, 31, 32, 36, 28, 16, -1, 33, 17, 29, 37, -1, 32, 33, 37, 36, -1, 36, 37, 29, 28, -1, 28, 29, 17, 16, -1, 16, 17, 33, 32, 31, 33, 37, 29, 17, -1, 34, 18, 30, 38, -1, 33, 34, 38, 37, -1, 37, 38, 30, 29, -1, 29, 30, 18, 17, -1, 17, 18, 34, 33, 31, 34, 38, 30, 18, -1, 35, 19, 31, 39, -1, 34, 35, 39, 38, -1, 38, 39, 31, 30, -1, 30, 31, 19, 18, -1, 18, 19, 35, 34, 31, 6, 5, 4, 7, -1, 40, 44, 32, 36, -1, 6, 40, 36, 5, -1, 5, 36, 32, 4, -1, 4, 32, 44, 7, -1, 7, 44, 40, 6, 31, 40, 36, 32, 44, -1, 41, 45, 33, 37, -1, 40, 41, 37, 36, -1, 36, 37, 33, 32, -1, 32, 33, 45, 44, -1, 44, 45, 41, 40, 31, 41, 37, 33, 45, -1, 42, 46, 34, 38, -1, 41, 42, 38, 37, -1, 37, 38, 34, 33, -1, 33, 34, 46, 45, -1, 45, 46, 42, 41, 31, 42, 38, 34, 46, -1, 43, 47, 35, 39, -1, 42, 43, 39, 38, -1, 38, 39, 35, 34, -1, 34, 35, 47, 46, -1, 46, 47, 43, 42, 31, 80, 9, 5, 6, -1, 82, 40, 36, 52, -1, 80, 82, 52, 9, -1, 9, 52, 36, 5, -1, 5, 36, 40, 6, -1, 6, 40, 82, 80, 31, 82, 52, 36, 40, -1, 83, 41, 37, 53, -1, 82, 83, 53, 52, -1, 52, 53, 37, 36, -1, 36, 37, 41, 40, -1, 40, 41, 83, 82, 31, 83, 53, 37, 41, -1, 84, 42, 38, 86, -1, 83, 84, 86, 53, -1, 53, 86, 38, 37, -1, 37, 38, 42, 41, -1, 41, 42, 84, 83, 31, 84, 86, 38, 42, -1, 85, 43, 39, 87, -1, 84, 85, 87, 86, -1, 86, 87, 39, 38, -1, 38, 39, 43, 42, -1, 42, 43, 85, 84, 31, 10, 9, 8, 11, -1, 56, 60, 48, 52, -1, 10, 56, 52, 9, -1, 9, 52, 48, 8, -1, 8, 48, 60, 11, -1, 11, 60, 56, 10, 31, 56, 52, 48, 60, -1, 57, 61, 49, 53, -1, 56, 57, 53, 52, -1, 52, 53, 49, 48, -1, 48, 49, 61, 60, -1, 60, 61, 57, 56, 31, 57, 53, 49, 61, -1, 58, 62, 50, 54, -1, 57, 58, 54, 53, -1, 53, 54, 50, 49, -1, 49, 50, 62, 61, -1, 61, 62, 58, 57, 31, 58, 54, 50, 62, -1, 59, 63, 51, 55, -1, 58, 59, 55, 54, -1, 54, 55, 51, 50, -1, 50, 51, 63, 62, -1, 62, 63, 59, 58, 31, 12, 13, 9, 10, -1, 64, 56, 52, 68, -1, 12, 64, 68, 13, -1, 13, 68, 52, 9, -1, 9, 52, 56, 10, -1, 10, 56, 64, 12, 31, 64, 68, 52, 56, -1, 65, 57, 53, 69, -1, 64, 65, 69, 68, -1, 68, 69, 53, 52, -1, 52, 53, 57, 56, -1, 56, 57, 65, 64, 31, 65, 69, 53, 57, -1, 66, 58, 54, 70, -1, 65, 66, 70, 69, -1, 69, 70, 54, 53, -1, 53, 54, 58, 57, -1, 57, 58, 66, 65, 31, 66, 70, 54, 58, -1, 67, 59, 55, 71, -1, 66, 67, 71, 70, -1, 70, 71, 55, 54, -1, 54, 55, 59, 58, -1, 58, 59, 67, 66, 31, 14, 15, 13, 12, -1, 72, 64, 68, 76, -1, 14, 72, 76, 15, -1, 15, 76, 68, 13, -1, 13, 68, 64, 12, -1, 12, 64, 72, 14, 31, 72, 76, 68, 64, -1, 73, 65, 69, 77, -1, 72, 73, 77, 76, -1, 76, 77, 69, 68, -1, 68, 69, 65, 64, -1, 64, 65, 73, 72, 31, 73, 77, 69, 65, -1, 74, 66, 70, 78, -1, 73, 74, 78, 77, -1, 77, 78, 70, 69, -1, 69, 70, 66, 65, -1, 65, 66, 74, 73, 31, 74, 78, 70, 66, -1, 75, 67, 71, 79, -1, 74, 75, 79, 78, -1, 78, 79, 71, 70, -1, 70, 71, 67, 66, -1, 66, 67, 75, 74, 31, 2, 3, 13, 81, -1, 24, 90, 68, 28, -1, 2, 24, 28, 3, -1, 3, 28, 68, 13, -1, 13, 68, 90, 81, -1, 81, 90, 24, 2, 31, 24, 28, 68, 90, -1, 25, 91, 69, 29, -1, 24, 25, 29, 28, -1, 28, 29, 69, 68, -1, 68, 69, 91, 90, -1, 90, 91, 25, 24, 31, 25, 29, 69, 91, -1, 26, 92, 88, 30, -1, 25, 26, 30, 29, -1, 29, 30, 88, 69, -1, 69, 88, 92, 91, -1, 91, 92, 26, 25, 31, 26, 30, 88, 92, -1, 27, 93, 89, 31, -1, 26, 27, 31, 30, -1, 30, 31, 89, 88, -1, 88, 89, 93, 92, -1, 92, 93, 27, 26, 31, 13, 3, 5, 9, -1, 68, 52, 36, 28, -1, 13, 68, 28, 3, -1, 3, 28, 36, 5, -1, 5, 36, 52, 9, -1, 9, 52, 68, 13, 31, 68, 28, 36, 52, -1, 69, 53, 37, 29, -1, 68, 69, 29, 28, -1, 28, 29, 37, 36, -1, 36, 37, 53, 52, -1, 52, 53, 69, 68, 31, 69, 29, 37, 53, -1, 88, 86, 38, 30, -1, 69, 88, 30, 29, -1, 29, 30, 38, 37, -1, 37, 38, 86, 53, -1, 53, 86, 88, 69, 31, 88, 30, 38, 86, -1, 89, 87, 39, 31, -1, 88, 89, 31, 30, -1, 30, 31, 39, 38, -1, 38, 39, 87, 86, -1, 86, 87, 89, 88] + ) + cI = DataArrayInt( + [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600, 630, 660, 690, 720, 750, 780, 810, 840, 870, 900, 930, 960, 990, 1020, 1050, 1080] + ) + # fmt: on + m3 = MEDCouplingUMesh("box", 3) + m3.setCoords(coo) + m3.setConnectivity(c, cI) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([2, 7, 12, 17, 95, 99, 103, 107, 129, 133, 137, 141]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m3.getNumberOfNodes() + + c2o2nN = mfu.crackAlong("group") + + self.assertEqual( + { + 7: {19: 115, 31: 114, 18: 112, 30: 111}, + 27: {70: 116, 71: 99, 66: 95, 67: 94}, + 26: {70: 116, 65: 98, 69: 96, 66: 95}, + 25: {64: 118, 68: 106, 65: 98, 69: 96}, + 29: {68: 107, 69: 97}, + 30: {69: 97}, + 24: {64: 118, 12: 117, 68: 106, 13: 100}, + 28: {68: 107, 13: 101}, + 32: {28: 108, 3: 103}, + 35: {31: 114, 89: 113, 30: 111, 88: 110}, + 6: {18: 112, 30: 111, 17: 109, 29: 105}, + 4: {28: 108, 0: 104, 3: 103, 16: 102}, + 33: {28: 108, 29: 105}, + 5: {17: 109, 28: 108, 29: 105, 16: 102}, + 34: {30: 111, 88: 110, 29: 105}, + }, + c2o2nN, + ) + m3_bis = mfu.getMeshAtLevel(0) + m3_bis.checkConsistency() + m2_bis = mfu.getMeshAtLevel(-1) + m2_bis.checkConsistency() + self.assertEqual(nNod + 25, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 25, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod + 25, m2_bis.getNumberOfNodes()) + # self.assertEqual( + # [0, 3, 12, 13, 16, 17, 18, 19, 28, 29, 30, 31, 64, 65, 66, 67, 68, 69, 70, 71, 88, 89], + # nodesDup.getValues(), + # ) # fmt: skip + # self.assertEqual( + # m3_bis.getCoords()[nodesDup].getValues(), + # m3_bis.getCoords()[nNod:].getValues(), + # ) + # self.assertEqual( + # set([0, 1, 2, 3, 24, 25, 26, 27, 28, 29, 30, 31]), set(cells1.getValues()) + # ) + # self.assertEqual( + # set([4, 5, 6, 7, 20, 21, 22, 23, 32, 33, 34, 35]), set(cells2.getValues()) + # ) + self.assertEqual( + [2 , 7 , 12 , 17 , 95 , 99 , 103, 107, 129, 133, 137, 141, + 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162], + mfu.getGroupArr(-1, "group").getValues(), + ) # fmt: skip + m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity() + m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) + + # def testBuildInnerBoundary8(self): + # """3D test where the crack leaves 'naked' cells. If we call a 'close-to-crack cell' a cell which shares a face with the M1 group, + # a 'naked cell' is a cell that has some node duplicated, but which do not share any face with a 'close-to-crack cell'. In this case + # it is tricky to decide whether this cell should be renumbered or not ... + # Warning: on the mesh below some points have already been doubled by a previous cut. + # """ + # # fmt: off + # coo = DataArrayDouble( + # [(0, 15, 0), (0, 5, 0), (3, 5, 0), (5, 5, 0), (5, 15, 0), (5, 20, 0), (0, 20, 0), (15, 20, 0), (15, 15, 0), (20, 15, 0), (20, 20, 0), (20, 5, 0), (15, 5, 0), (15, 0, 0), (20, 0, 0), (5, -1.60551e-25, 0), (5, 3, 0), (3, 0, 0), (3, 3, 0), (0, 0, 0), (0, 3, 0), (0, 15, 10), (0, 15, 20), (0, 15, 30), (0, 15, 40), (0, 5, 10), (0, 5, 20), (0, 5, 30), (0, 5, 40), (3, 5, 10), (3, 5, 20), (3, 5, 30), (3, 5, 40), (5, 5, 10), (5, 5, 20), (5, 5, 30), (5, 5, 40), (5, 15, 10), (5, 15, 20), (5, 15, 30), (5, 15, 40), (5, 20, 10), (5, 20, 20), (5, 20, 30), (5, 20, 40), (0, 20, 10), (0, 20, 20), (0, 20, 30), (0, 20, 40), (15, 20, 10), (15, 20, 20), (15, 20, 30), (15, 20, 40), (15, 15, 10), (15, 15, 20), (15, 15, 30), (15, 15, 40), (20, 15, 10), (20, 15, 20), (20, 15, 30), (20, 15, 40), (20, 20, 10), (20, 20, 20), (20, 20, 30), (20, 20, 40), (20, 5, 10), (20, 5, 20), (20, 5, 30), (20, 5, 40), (15, 5, 10), (15, 5, 20), (15, 5, 30), (15, 5, 40), (15, 0, 10), (15, 0, 20), (15, 0, 30), (15, 0, 40), (20, 0, 10), (20, 0, 20), (20, 0, 30), (20, 0, 40), (5, -1.60551e-25, 10), (5, -1.60551e-25, 20), (5, -1.60551e-25, 30), (5, -1.60551e-25, 40), (5, 3, 10), (5, 3, 20), (5, 3, 30), (5, 3, 40), (3, 0, 10), (3, 0, 20), (3, 0, 30), (3, 0, 40), (3, 3, 10), (3, 3, 20), (3, 3, 30), (3, 3, 40), (0, 0, 10), (0, 0, 20), (0, 0, 30), (0, 0, 40), (0, 3, 10), (0, 3, 20), (0, 3, 30), (0, 3, 40), (0, 9, 0), (3, 9, 0), (20, 9, 0), (0, 9, 10), (0, 9, 20), (0, 9, 30), (0, 9, 40), (3, 9, 10), (3, 9, 20), (3, 9, 30), (3, 9, 40), (5, 9, 30), (5, 9, 40), (20, 9, 10), (20, 9, 20), (20, 9, 30), (20, 9, 40), (15, 9, 30), (15, 9, 40), (0, 15, 0), (20, 15, 0), (0, 15, 10), (0, 15, 20), (0, 15, 30), (0, 15, 40), (5, 15, 30), (5, 15, 40), (15, 15, 30), (15, 15, 40), (20, 15, 10), (20, 15, 20), (20, 15, 30), (20, 15, 40)] + # ) + # c = DataArrayInt( + # [31, 5, 4, 124, 6, -1, 41, 45, 126, 37, -1, 5, 41, 37, 4, -1, 4, 37, 126, 124, -1, 124, 126, 45, 6, -1, 6, 45, 41, 5, 31, 41, 37, 126, 45, -1, 42, 46, 127, 38, -1, 41, 42, 38, 37, -1, 37, 38, 127, 126, -1, 126, 127, 46, 45, -1, 45, 46, 42, 41, 31, 42, 38, 127, 46, -1, 43, 47, 128, 130, -1, 42, 43, 130, 38, -1, 38, 130, 128, 127, -1, 127, 128, 47, 46, -1, 46, 47, 43, 42, 31, 43, 130, 128, 47, -1, 44, 48, 129, 131, -1, 43, 44, 131, 130, -1, 130, 131, 129, 128, -1, 128, 129, 48, 47, -1, 47, 48, 44, 43, 31, 7, 8, 4, 5, -1, 49, 41, 37, 53, -1, 7, 49, 53, 8, -1, 8, 53, 37, 4, -1, 4, 37, 41, 5, -1, 5, 41, 49, 7, 31, 49, 53, 37, 41, -1, 50, 42, 38, 54, -1, 49, 50, 54, 53, -1, 53, 54, 38, 37, -1, 37, 38, 42, 41, -1, 41, 42, 50, 49, 31, 50, 54, 38, 42, -1, 51, 43, 130, 132, -1, 50, 51, 132, 54, -1, 54, 132, 130, 38, -1, 38, 130, 43, 42, -1, 42, 43, 51, 50, 31, 51, 132, 130, 43, -1, 52, 44, 131, 133, -1, 51, 52, 133, 132, -1, 132, 133, 131, 130, -1, 130, 131, 44, 43, -1, 43, 44, 52, 51, 31, 125, 8, 7, 10, -1, 134, 61, 49, 53, -1, 125, 134, 53, 8, -1, 8, 53, 49, 7, -1, 7, 49, 61, 10, -1, 10, 61, 134, 125, 31, 134, 53, 49, 61, -1, 135, 62, 50, 54, -1, 134, 135, 54, 53, -1, 53, 54, 50, 49, -1, 49, 50, 62, 61, -1, 61, 62, 135, 134, 31, 135, 54, 50, 62, -1, 136, 63, 51, 132, -1, 135, 136, 132, 54, -1, 54, 132, 51, 50, -1, 50, 51, 63, 62, -1, 62, 63, 136, 135, 31, 136, 132, 51, 63, -1, 137, 64, 52, 133, -1, 136, 137, 133, 132, -1, 132, 133, 52, 51, -1, 51, 52, 64, 63, -1, 63, 64, 137, 136, 31, 107, 12, 8, 9, -1, 118, 57, 53, 69, -1, 107, 118, 69, 12, -1, 12, 69, 53, 8, -1, 8, 53, 57, 9, -1, 9, 57, 118, 107, 31, 118, 69, 53, 57, -1, 119, 58, 54, 70, -1, 118, 119, 70, 69, -1, 69, 70, 54, 53, -1, 53, 54, 58, 57, -1, 57, 58, 119, 118, 31, 119, 70, 54, 58, -1, 120, 59, 55, 122, -1, 119, 120, 122, 70, -1, 70, 122, 55, 54, -1, 54, 55, 59, 58, -1, 58, 59, 120, 119, 31, 120, 122, 55, 59, -1, 121, 60, 56, 123, -1, 120, 121, 123, 122, -1, 122, 123, 56, 55, -1, 55, 56, 60, 59, -1, 59, 60, 121, 120, 31, 13, 12, 11, 14, -1, 73, 77, 65, 69, -1, 13, 73, 69, 12, -1, 12, 69, 65, 11, -1, 11, 65, 77, 14, -1, 14, 77, 73, 13, 31, 73, 69, 65, 77, -1, 74, 78, 66, 70, -1, 73, 74, 70, 69, -1, 69, 70, 66, 65, -1, 65, 66, 78, 77, -1, 77, 78, 74, 73, 31, 74, 70, 66, 78, -1, 75, 79, 67, 71, -1, 74, 75, 71, 70, -1, 70, 71, 67, 66, -1, 66, 67, 79, 78, -1, 78, 79, 75, 74, 31, 75, 71, 67, 79, -1, 76, 80, 68, 72, -1, 75, 76, 72, 71, -1, 71, 72, 68, 67, -1, 67, 68, 80, 79, -1, 79, 80, 76, 75, 31, 17, 18, 16, 15, -1, 89, 81, 85, 93, -1, 17, 89, 93, 18, -1, 18, 93, 85, 16, -1, 16, 85, 81, 15, -1, 15, 81, 89, 17, 31, 89, 93, 85, 81, -1, 90, 82, 86, 94, -1, 89, 90, 94, 93, -1, 93, 94, 86, 85, -1, 85, 86, 82, 81, -1, 81, 82, 90, 89, 31, 90, 94, 86, 82, -1, 91, 83, 87, 95, -1, 90, 91, 95, 94, -1, 94, 95, 87, 86, -1, 86, 87, 83, 82, -1, 82, 83, 91, 90, 31, 91, 95, 87, 83, -1, 92, 84, 88, 96, -1, 91, 92, 96, 95, -1, 95, 96, 88, 87, -1, 87, 88, 84, 83, -1, 83, 84, 92, 91, 31, 19, 20, 18, 17, -1, 97, 89, 93, 101, -1, 19, 97, 101, 20, -1, 20, 101, 93, 18, -1, 18, 93, 89, 17, -1, 17, 89, 97, 19, 31, 97, 101, 93, 89, -1, 98, 90, 94, 102, -1, 97, 98, 102, 101, -1, 101, 102, 94, 93, -1, 93, 94, 90, 89, -1, 89, 90, 98, 97, 31, 98, 102, 94, 90, -1, 99, 91, 95, 103, -1, 98, 99, 103, 102, -1, 102, 103, 95, 94, -1, 94, 95, 91, 90, -1, 90, 91, 99, 98, 31, 99, 103, 95, 91, -1, 100, 92, 96, 104, -1, 99, 100, 104, 103, -1, 103, 104, 96, 95, -1, 95, 96, 92, 91, -1, 91, 92, 100, 99, 31, 1, 2, 18, 20, -1, 25, 101, 93, 29, -1, 1, 25, 29, 2, -1, 2, 29, 93, 18, -1, 18, 93, 101, 20, -1, 20, 101, 25, 1, 31, 25, 29, 93, 101, -1, 26, 102, 94, 30, -1, 25, 26, 30, 29, -1, 29, 30, 94, 93, -1, 93, 94, 102, 101, -1, 101, 102, 26, 25, 31, 26, 30, 94, 102, -1, 27, 103, 95, 31, -1, 26, 27, 31, 30, -1, 30, 31, 95, 94, -1, 94, 95, 103, 102, -1, 102, 103, 27, 26, 31, 27, 31, 95, 103, -1, 28, 104, 96, 32, -1, 27, 28, 32, 31, -1, 31, 32, 96, 95, -1, 95, 96, 104, 103, -1, 103, 104, 28, 27, 31, 3, 4, 8, 12, -1, 33, 69, 53, 37, -1, 3, 33, 37, 4, -1, 4, 37, 53, 8, -1, 8, 53, 69, 12, -1, 12, 69, 33, 3, 31, 33, 37, 53, 69, -1, 34, 70, 54, 38, -1, 33, 34, 38, 37, -1, 37, 38, 54, 53, -1, 53, 54, 70, 69, -1, 69, 70, 34, 33, 31, 34, 38, 54, 70, -1, 116, 122, 55, 39, -1, 34, 116, 39, 38, -1, 38, 39, 55, 54, -1, 54, 55, 122, 70, -1, 70, 122, 116, 34, 31, 116, 39, 55, 122, -1, 117, 123, 56, 40, -1, 116, 117, 40, 39, -1, 39, 40, 56, 55, -1, 55, 56, 123, 122, -1, 122, 123, 117, 116, 31, 16, 18, 2, 3, -1, 85, 33, 29, 93, -1, 16, 85, 93, 18, -1, 18, 93, 29, 2, -1, 2, 29, 33, 3, -1, 3, 33, 85, 16, 31, 85, 93, 29, 33, -1, 86, 34, 30, 94, -1, 85, 86, 94, 93, -1, 93, 94, 30, 29, -1, 29, 30, 34, 33, -1, 33, 34, 86, 85, 31, 86, 94, 30, 34, -1, 87, 35, 31, 95, -1, 86, 87, 95, 94, -1, 94, 95, 31, 30, -1, 30, 31, 35, 34, -1, 34, 35, 87, 86, 31, 87, 95, 31, 35, -1, 88, 36, 32, 96, -1, 87, 88, 96, 95, -1, 95, 96, 32, 31, -1, 31, 32, 36, 35, -1, 35, 36, 88, 87, 31, 4, 3, 106, 105, 0, -1, 37, 21, 108, 112, 33, -1, 3, 4, 37, 33, -1, 106, 3, 33, 112, -1, 105, 106, 112, 108, -1, 0, 105, 108, 21, -1, 4, 0, 21, 37, 31, 37, 33, 112, 108, 21, -1, 38, 22, 109, 113, 34, -1, 33, 37, 38, 34, -1, 112, 33, 34, 113, -1, 108, 112, 113, 109, -1, 21, 108, 109, 22, -1, 37, 21, 22, 38, 31, 38, 34, 113, 109, 22, -1, 39, 23, 110, 114, 116, -1, 34, 38, 39, 116, -1, 113, 34, 116, 114, -1, 109, 113, 114, 110, -1, 22, 109, 110, 23, -1, 38, 22, 23, 39, 31, 39, 116, 114, 110, 23, -1, 40, 24, 111, 115, 117, -1, 116, 39, 40, 117, -1, 114, 116, 117, 115, -1, 110, 114, 115, 111, -1, 23, 110, 111, 24, -1, 39, 23, 24, 40, 31, 16, 3, 12, 13, 15, -1, 85, 81, 73, 69, 33, -1, 3, 16, 85, 33, -1, 12, 3, 33, 69, -1, 13, 12, 69, 73, -1, 15, 13, 73, 81, -1, 16, 15, 81, 85, 31, 85, 33, 69, 73, 81, -1, 86, 82, 74, 70, 34, -1, 33, 85, 86, 34, -1, 69, 33, 34, 70, -1, 73, 69, 70, 74, -1, 81, 73, 74, 82, -1, 85, 81, 82, 86, 31, 86, 34, 70, 74, 82, -1, 87, 83, 75, 71, 35, -1, 34, 86, 87, 35, -1, 70, 34, 35, 71, -1, 74, 70, 71, 75, -1, 82, 74, 75, 83, -1, 86, 82, 83, 87, 31, 87, 35, 71, 75, 83, -1, 88, 84, 76, 72, 36, -1, 35, 87, 88, 36, -1, 71, 35, 36, 72, -1, 75, 71, 72, 76, -1, 83, 75, 76, 84, -1, 87, 83, 84, 88] + # ) + # cI = DataArrayInt( + # [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450, 480, 510, 540, 570, 600, 630, 660, 690, 720, 750, 780, 810, 840, 870, 900, 930, 960, 990, 1020, 1050, 1080, 1110, 1140, 1170, 1200, 1237, 1274, 1311, 1348, 1385, 1422, 1459, 1496] + # ) + # # fmt: on + # m3 = MEDCouplingUMesh("box", 3) + # m3.setCoords(coo) + # m3.setConnectivity(c, cI) + # m3.mergeNodes(1e-12) + # m3.checkConsistency() + # m2, _, _, _, _ = m3.buildDescendingConnectivity() + # # grpIds = DataArrayInt( + # # [2, 7, 12, 17, 101, 106, 111, 116, 160, 164, 170, 173, 176, 179] + # # ) + # grpIds = DataArrayInt( + # [16, 149, 12, 146, 8, 143, 137, 108, 134, 103, 131, 98, 128, 93] + # ) + # grpIds.setName("group") + # mfu = MEDFileUMesh() + # mfu.setMeshAtLevel(0, m3) + # mfu.setMeshAtLevel(-1, m2) + # mfu.setGroupsAtLevel(-1, [grpIds]) + # nNod = m3.getNumberOfNodes() + + # mfu.write("case8dupnode_in.med", 2) + + # c2o2nN = mfu.crackAlong("group", False) + + # mfu.write("case8dupnode_out.med", 2) + + # mfu.OpenCrack(c2o2nN) + + # mfu.write("case8dupnode_cracked.med", 2) + + # self.assertEqual( + # { + # 44: {85: 166, 15: 163, 81: 159, 16: 158}, + # 6: {130: 162, 43: 161, 38: 160, 42: 157}, + # 41: {38: 160, 37: 154}, + # 47: {83: 146, 84: 143, 87: 148, 35: 139, 88: 145, 36: 138}, + # 46: {86: 150, 82: 149, 83: 146, 87: 148, 35: 139}, + # 7: {131: 165, 44: 164, 130: 162, 43: 161}, + # 36: {33: 155, 3: 140}, + # 43: {116: 147, 117: 144, 39: 142, 40: 141}, + # 33: {38: 160, 37: 154}, + # 4: {4: 156, 37: 154, 41: 153, 5: 152}, + # 42: {38: 160, 116: 147, 39: 142}, + # 45: {85: 166, 81: 159, 86: 150, 82: 149}, + # 32: {4: 156, 37: 154}, + # 37: {33: 155, 34: 151}, + # 38: {34: 151}, + # 34: {38: 160}, + # 5: {38: 160, 42: 157, 37: 154, 41: 153}, + # 40: {4: 156, 37: 154}, + # }, + # c2o2nN, + # ) + # m3_bis = mfu.getMeshAtLevel(0) + # m3_bis.checkConsistency() + # m2_bis = mfu.getMeshAtLevel(-1) + # m2_bis.checkConsistency() + # self.assertEqual(nNod + 23, mfu.getNumberOfNodes()) + # self.assertEqual(nNod + 23, m3_bis.getNumberOfNodes()) + # self.assertEqual(nNod + 23, m2_bis.getNumberOfNodes()) + # self.assertEqual( + # [5, 15, 16, 35, 36, 39, 40, 41, 42, 43, 44, 81, 82, 83, 84, 85, 86, 87, 88, 116, 117, 130, 131], + # nodesDup.getValues(), + # ) # fmt: skip + # self.assertEqual( + # m3_bis.getCoords()[nodesDup].getValues(), + # m3_bis.getCoords()[nNod:].getValues(), + # ) + # self.assertEqual( + # set([0, 1, 2, 3, 20, 21, 22, 23, 34, 35, 36, 37, 38, 39]), + # set(cells1.getValues()), + # ) + # self.assertEqual( + # set([4, 5, 6, 7, 42, 43, 44, 45, 46, 47]), set(cells2.getValues()) + # ) + # self.assertEqual( + # [2, 7, 12, 17, 101, 106, 111, 116, 160, 164, 170, 173, 176, 179], + # mfu.getGroupArr(-1, "group").getValues(), + # ) + # self.assertEqual( + # [212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225], + # mfu.getGroupArr(-1, "group_dup").getValues(), + # ) # fmt: skip + # m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity() + # m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) + + def testBuildInnerBoundary9(self): + """3D test where the crack is performed so that two non-connex parts are found facing one single connex part on the other side + of the crack. + """ + # fmt: off + coo = DataArrayDouble( + [(0, 4.6, 0), (3, 4.6, 0), (5, 4.6, 0), (15, 4.6, 0), (15, 0, 0), (5, -1.60551e-25, 0), (5, 3, 0), (3, 0, 0), (3, 3.8, 0), (0, 0, 0), (0, 3.8, 0), (0, 4.6, 10), (0, 4.6, 20), (3, 4.6, 10), (3, 4.6, 20), (5, 4.6, 10), (5, 4.6, 20), (15, 4.6, 10), (15, 4.6, 20), (15, 0, 10), (15, 0, 20), (5, -1.60551e-25, 10), (5, -1.60551e-25, 20), (5, 3, 10), (5, 3, 20), (3, 0, 10), (3, 0, 20), (3, 3.8, 10), (3, 3.8, 20), (0, 0, 10), (0, 0, 20), (0, 3.8, 10), (0, 3.8, 20), (3, 3, 0), (0, 3, 0), (3, 3, 10), (3, 3, 20), (0, 3, 10), (0, 3, 20), ] + ) + c = DataArrayInt( + [31, 7, 33, 6, 5, -1, 25, 21, 23, 35, -1, 7, 25, 35, 33, -1, 33, 35, 23, 6, -1, 6, 23, 21, 5, -1, 5, 21, 25, 7, 31, 25, 35, 23, 21, -1, 26, 22, 24, 36, -1, 25, 26, 36, 35, -1, 35, 36, 24, 23, -1, 23, 24, 22, 21, -1, 21, 22, 26, 25, 31, 9, 34, 33, 7, -1, 29, 25, 35, 37, -1, 9, 29, 37, 34, -1, 34, 37, 35, 33, -1, 33, 35, 25, 7, -1, 7, 25, 29, 9, 31, 29, 37, 35, 25, -1, 30, 26, 36, 38, -1, 29, 30, 38, 37, -1, 37, 38, 36, 35, -1, 35, 36, 26, 25, -1, 25, 26, 30, 29, 31, 0, 1, 8, 10, -1, 11, 31, 27, 13, -1, 0, 11, 13, 1, -1, 1, 13, 27, 8, -1, 8, 27, 31, 10, -1, 10, 31, 11, 0, 31, 11, 13, 27, 31, -1, 12, 32, 28, 14, -1, 11, 12, 14, 13, -1, 13, 14, 28, 27, -1, 27, 28, 32, 31, -1, 31, 32, 12, 11, 31, 6, 8, 1, 2, -1, 23, 15, 13, 27, -1, 6, 23, 27, 8, -1, 8, 27, 13, 1, -1, 1, 13, 15, 2, -1, 2, 15, 23, 6, 31, 23, 27, 13, 15, -1, 24, 16, 14, 28, -1, 23, 24, 28, 27, -1, 27, 28, 14, 13, -1, 13, 14, 16, 15, -1, 15, 16, 24, 23, 31, 6, 2, 3, 4, 5, -1, 23, 21, 19, 17, 15, -1, 2, 6, 23, 15, -1, 3, 2, 15, 17, -1, 4, 3, 17, 19, -1, 5, 4, 19, 21, -1, 6, 5, 21, 23, 31, 23, 15, 17, 19, 21, -1, 24, 22, 20, 18, 16, -1, 15, 23, 24, 16, -1, 17, 15, 16, 18, -1, 19, 17, 18, 20, -1, 21, 19, 20, 22, -1, 23, 21, 22, 24] + ) + # fmt: on + cI = DataArrayInt([0, 30, 60, 90, 120, 150, 180, 210, 240, 277, 314]) + m3 = MEDCouplingUMesh("box", 3) + m3.setCoords(coo) + m3.setConnectivity(c, cI) + m3.checkConsistency() + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([4, 9, 35, 39]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + mfu.setGroupsAtLevel(-1, [grpIds]) + m2, _, _, _, _ = m3.buildDescendingConnectivity() + grpIds = DataArrayInt([4, 9, 35, 39]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m3) + mfu.setMeshAtLevel(-1, m2) + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m3.getNumberOfNodes() + + c2o2nN = mfu.crackAlong("group") + + self.assertEqual({6: {6: 49, 23: 47}, 7: {23: 47, 24: 43}, 8: {6: 50, 23: 48, 21: 46, 5: 45, 2: 41, 15: 40}, 9: {23: 48, 21: 46, 24: 44, 22: 42, 15: 40, 16: 39}}, c2o2nN) + + m3_bis = mfu.getMeshAtLevel(0) + m3_bis.checkConsistency() + m2_bis = mfu.getMeshAtLevel(-1) + m2_bis.checkConsistency() + self.assertEqual(nNod + 12, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 12, m3_bis.getNumberOfNodes()) + self.assertEqual(nNod + 12, m2_bis.getNumberOfNodes()) + # self.assertEqual([2, 5, 6, 15, 16, 21, 22, 23, 24], nodesDup.getValues()) + # self.assertEqual( + # m3_bis.getCoords()[nodesDup].getValues(), + # m3_bis.getCoords()[nNod:].getValues(), + # ) + # self.assertEqual(set([0, 1, 6, 7]), set(cells1.getValues())) + # self.assertEqual(set([8, 9]), set(cells2.getValues())) + self.assertEqual([4, 9, 35, 39, 49, 50, 51, 52], mfu.getGroupArr(-1, "group").getValues()) + m_desc, _, _, _, _ = m3_bis.buildDescendingConnectivity() + m_desc.checkDeepEquivalOnSameNodesWith(m2_bis, 2, 9.9999) + + def testBuildInnerBoundary10(self): + """2D tests where some cells are touching the M1 group with just a node, and are **not** neighbor + with any cells touching the M1 group by a face. + """ + # fmt: off + coo = DataArrayDouble( + [0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 1.0, 3.0, 1.0, 0.0, 2.0, 1.0, 2.0, 2.0, 2.0, 3.0, 2.0, 2.5, 0.0, 3.0, 0.5, ], + 14, + 2, + ) + c = DataArrayInt( + [3, 12, 2, 6, 3, 3, 12, 6, 3, 13, 3, 6, 3, 7, 13, 6, 4, 1, 0, 4, 5, 4, 2, 1, 5, 6, 4, 5, 4, 8, 9, 4, 6, 5, 9, 10, 4, 7, 6, 10, 11] + ) + # fmt: on + cI = DataArrayInt([0, 4, 8, 12, 16, 21, 26, 31, 36, 41]) + m2 = MEDCouplingUMesh("mesh", 2) + m2.setCoords(coo) + m2.setConnectivity(c, cI) + m2.checkConsistency() + m1, _, _, _, _ = m2.buildDescendingConnectivity() + grpIds = DataArrayInt([8, 14]) + grpIds.setName("group") + mfu = MEDFileUMesh() + mfu.setMeshAtLevel(0, m2) + mfu.setMeshAtLevel(-1, m1) + mfu.setGroupsAtLevel(-1, [grpIds]) + nNod = m2.getNumberOfNodes() + + c2o2nN = mfu.crackAlong("group") + + self.assertEqual({7: {6: 15}, 8: {6: 15, 7: 14}}, c2o2nN) + + m2_bis = mfu.getMeshAtLevel(0) + m2_bis.checkConsistency() + m1_bis = mfu.getMeshAtLevel(-1) + m1_bis.checkConsistency() + self.assertEqual(nNod + 2, mfu.getNumberOfNodes()) + self.assertEqual(nNod + 2, m2_bis.getNumberOfNodes()) + self.assertEqual(nNod + 2, m1_bis.getNumberOfNodes()) + # self.assertEqual([6, 7], nodesDup.getValues()) + # self.assertEqual([2.0, 1.0, 3.0, 1.0], m2_bis.getCoords()[nNod:].getValues()) + # self.assertEqual(set([0, 1, 2, 3, 5]), set(cells1.getValues())) + # self.assertEqual(set([7, 8]), set(cells2.getValues())) + self.assertEqual([8, 14, 22, 23], mfu.getGroupArr(-1, "group").getValues()) + + +if __name__ == "__main__": + unittest.main() diff --git a/src/MEDLoader/Swig/MEDLoaderCommon.i b/src/MEDLoader/Swig/MEDLoaderCommon.i index aebbe75c3..6f97b005f 100644 --- a/src/MEDLoader/Swig/MEDLoaderCommon.i +++ b/src/MEDLoader/Swig/MEDLoaderCommon.i @@ -1729,6 +1729,60 @@ namespace MEDCoupling return ret; } + PyObject * crackAlong(const std::string& grpNameM1, bool grpMustBeFullyDup=true) + { + const std::unordered_map> cellOld2NewNode = self->crackAlong(grpNameM1, grpMustBeFullyDup); + PyObject * cellOld2NewNodePy = PyDict_New(); + for (const auto & cellIt: cellOld2NewNode) { + PyObject * old2NewPy = PyDict_New(); + for (const auto & keyVal: cellIt.second) { + PyDict_SetItem(old2NewPy, PyInt_FromLong(keyVal.first), PyInt_FromLong(keyVal.second)); + } + PyDict_SetItem(cellOld2NewNodePy, PyInt_FromLong(cellIt.first), old2NewPy); + } + return cellOld2NewNodePy; + } + + // void OpenCrack(const std::unordered_map>& c2o2nN, double factor=0.9) + // { + // self->OpenCrack(c2o2nN, factor); + // } + + void openCrack(PyObject* c2o2nNPy, double factor=0.9) + { + if (!PyDict_Check(c2o2nNPy)) + throw INTERP_KERNEL::Exception("c2o2nNPy is not a python dict"); + + std::unordered_map> c2o2nN; + Py_ssize_t i_cells = 0; + PyObject *cell, *val; + + while (PyDict_Next(c2o2nNPy, &i_cells, &cell, &val)) { + if (!PyInt_Check(cell)) + throw INTERP_KERNEL::Exception("One of c2o2nNPy cell id is not an int."); + + if (!PyDict_Check(val)) + throw INTERP_KERNEL::Exception("One of c2o2nNPy val is not a python dict"); + std::unordered_map o2nN {}; + Py_ssize_t i_node = 0; + PyObject *oldN, *newN; + + while (PyDict_Next(val, &i_node, &oldN, &newN)) { + if (!PyInt_Check(oldN)) + throw INTERP_KERNEL::Exception( + "One of c2o2nNPy old node id is not an int."); + if (!PyInt_Check(newN)) + throw INTERP_KERNEL::Exception( + "One of c2o2nNPy new node id is not an int."); + o2nN[PyInt_AsLong(oldN)] = PyInt_AsLong(newN); + } + + c2o2nN[PyInt_AsLong(cell)] = o2nN; + } + + self->openCrack(c2o2nN, factor); + } + MEDCoupling1GTUMesh *getDirectUndergroundSingleGeoTypeMesh(INTERP_KERNEL::NormalizedCellType gt) const { MEDCoupling1GTUMesh *ret(self->getDirectUndergroundSingleGeoTypeMesh(gt)); diff --git a/src/MEDLoader/Swig/tests.set b/src/MEDLoader/Swig/tests.set index 2e58e3739..f57a3e06b 100644 --- a/src/MEDLoader/Swig/tests.set +++ b/src/MEDLoader/Swig/tests.set @@ -23,6 +23,7 @@ SET(BASE_TESTS MEDLoaderExamplesTest.py UsersGuideExamplesTest_ML.py MeshFormatWriterTest1.py + CrackAlongTest.py ) # if numpy is used diff --git a/src/MEDLoader/Test/CMakeLists.txt b/src/MEDLoader/Test/CMakeLists.txt index b1b1b907e..274ce2fec 100644 --- a/src/MEDLoader/Test/CMakeLists.txt +++ b/src/MEDLoader/Test/CMakeLists.txt @@ -31,6 +31,7 @@ INCLUDE_DIRECTORIES( SET(TestMEDLoader_SOURCES TestMEDLoader.cxx + CrackAlgoTest.cxx MEDLoaderTest.cxx ) diff --git a/src/MEDLoader/Test/CrackAlgoTest.cxx b/src/MEDLoader/Test/CrackAlgoTest.cxx new file mode 100644 index 000000000..7fdb906d8 --- /dev/null +++ b/src/MEDLoader/Test/CrackAlgoTest.cxx @@ -0,0 +1,657 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Aymeric SONOLET (CEA/DES) + +#include "CrackAlgoTest.hxx" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using std::vector; +using std::string; +using std::pair; +using std::copy; + +using MEDCoupling::CrackAlgoTest; +using MEDCoupling::MEDFileUMesh; +using MEDCoupling::MEDCouplingCMesh; +using MEDCoupling::MEDCouplingUMesh; +using MEDCoupling::DataArrayIdType; +using MEDCoupling::MCAuto; + +using MFU = MCAuto; +using MCC = MCAuto; +using MCU = MCAuto; +using DAI = MCAuto; + + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ +--------+--------+ +--------+--------+ + * | | | | | | | o | + * | | | | | | | o | + * | | | | | | | +-+ | + * 1 +--------+--------+ +--------+--------+ +ooooooo|o|ooooooo+ + * | | | | | | | +-+ | + * | | | | | | | o | + * | | | | | | | o | + * 0 +--------+--------+ +--------+--------+ +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DHalfCrossCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[6] {1.0, 1.0, 1.0, 1.0, 0.1, 0.9}; + MCAuto g1 = m1->getCellsInBoundingBox(box, 0.01); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 4); + const auto res = TestCrack(mf, "M1", "halfCrossCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ ooooooooooooooooooo +--------+--------+ + * | o | ooooooooooooooooooo | o | + * | o | ooooooooooooooooooo | o | + * | +-+ | oooooooo+-+oooooooo | +-+ | + * 1 +ooooooo|o|ooooooo+ oooooooo|o|oooooooo +ooooooo|o|ooooooo+ + * | +-+ | oooooooo+-+oooooooo | +-+ | + * | o | ooooooooooooooooooo | o | + * | o | ooooooooooooooooooo | o | + * 0 +--------+--------+ ooooooooooooooooooo +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DFullFullCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[6] {1.0, 1.0, 1.0, 1.0, 0.1, 1.9}; + DAI g1 = m1->getCellsInBoundingBox(box, 0.01); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + const auto res = TestCrack(mf, "M1", "fullFullCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 12); +} + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ +--------+--------+ +--------+--------+ + * | o | | o | | o | + * | o | | o | | o | + * | +-+ | | o | | +-+ | + * 1 +ooooooo|o|ooooooo+ +ooooooooooooooooo+ +ooooooo|o|ooooooo+ + * | +-+ | | o | | +-+ | + * | o | | o | | o | + * | o | | o | | o | + * 0 +--------+--------+ +--------+--------+ +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DFullCrossCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[6] {1.0, 1.0, 1.0, 1.0, 0.1, 0.9}; + MCAuto g1 = m1->getCellsInBoundingBox(box, 0.01); + double box2[6] {1.0, 1.0, 1.0, 1.0, 1.1, 1.9}; + MCAuto g2 = m1->getCellsInBoundingBox(box, 0.01); + g1->aggregate(g2); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 8); + + const auto res = TestCrack(mf, "M1", "fullCrossCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ +--------+--------+ +--------+--------+ + * | | | | | | | | | + * | | | | | | | | | + * | | | | | | | +-+ | +-+ | + * 1 +--------+--------+ +--------+--------+ +oo|o|ooooooo|o|oo+ + * | | | | | | | +-+ | +-+ | + * | | | | | | | | | + * | | | | | | | | | + * 0 +--------+--------+ +--------+--------+ +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DHalfCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box1[6] {0.5, 0.5, 1.0, 1.0, 0.1, 0.9}; + DAI g1 = m1->getCellsInBoundingBox(box1, 0.01); + double box2[6] {1.5, 1.5, 1.0, 1.0, 0.1, 0.9}; + DAI g2 = m1->getCellsInBoundingBox(box2, 0.01); + + g1->aggregate(g2); + + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 2); + const auto res = TestCrack(mf, "M1", "halfCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ +--------+--------+ +--------+--------+ + * | | | | | | | | | + * | | | | | | | | | + * | +-+ | +-+ | | | | | +-+ | +-+ | + * 1 +oo|o|ooooooo|o|oo+ +ooooooooooooooooo+ +oo|o|ooooooo|o|oo+ + * | +-+ | +-+ | | | | | +-+ | +-+ | + * | | | | | | | | | + * | | | | | | | | | + * 0 +--------+--------+ +--------+--------+ +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DFullCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box1[6] {0.5, 0.5, 1.0, 1.0, 0.1, 0.9}; + DAI g1 = m1->getCellsInBoundingBox(box1, 0.01); + double box2[6] {1.5, 1.5, 1.0, 1.0, 0.1, 0.9}; + DAI g2 = m1->getCellsInBoundingBox(box2, 0.01); + double box3[6] {0.5, 0.5, 1.0, 1.0, 1.1, 1.9}; + DAI g3 = m1->getCellsInBoundingBox(box1, 0.01); + double box4[6] {1.5, 1.5, 1.0, 1.0, 1.1, 1.9}; + DAI g4 = m1->getCellsInBoundingBox(box2, 0.01); + + g1->aggregate(g2); + g1->aggregate(g3); + g1->aggregate(g4); + + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 4); + const auto res = TestCrack(mf, "M1", "fullCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +/* + * TOP z=2 MIDDLE z=1 BOTTOM z=0 + * 2 +--------+--------+ +--------+--------+ +--------+--------+ + * | | | | | | | | | + * | | | | | | | | | + * | +---+ | | | | | | | | + * 1 oo|ooo|oo+--------+ +--------+--------+ +--------+--------+ + * | +---+ | | | | | | | | + * | | | | | | | | | + * | | | | | | | | | + * 0 +--------+--------+ +--------+--------+ +--------+--------+ + * + * 0 1 2 + */ +void CrackAlgoTest::test3DAngleCut( +) { + MFU mf = make2x2Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[6] {0.1, 0.9, 1.0, 1.0, 1.1, 1.9}; + MCAuto g1 = m1->getCellsInBoundingBox(box, 0.01); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 1); + const auto res = TestCrack(mf, "M1", "angleCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +void CrackAlgoTest::testInnerCrossCut( +) { + MFU mf = make4x4Voxel(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[6] {2.0, 2.0, 2.0, 2.0, 2.0, 2.0}; + MCAuto g1 = m1->getCellsInBoundingBox(box, 0.01); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 12); + const auto res = TestCrack(mf, "M1", "innerCrossCut"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +void CrackAlgoTest::test2DGrp( +) { + MCAuto mf = make2DMesh(); + const auto res = TestCrack(mf, "Grp", "2dGrp"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +void CrackAlgoTest::test1DMesh( +) { + MFU mf = make1DMesh(); + + const MCU m1 = mf->getMeshAtLevel(-1); + double box[2] {1.9, 2.1}; + MCAuto g1 = m1->getCellsInBoundingBox(box, 0.01); + g1->setName("M1"); + mf->setGroupsAtLevel(-1, vector{g1}); + + CPPUNIT_ASSERT(g1->getNumberOfTuples() == 1); + const auto res = TestCrack(mf, "M1", "1DMesh"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +void CrackAlgoTest::test2DMeshNonConnexCut( +) { + MCAuto mf = make2DMesh2(); + const auto res = TestCrack(mf, "Grp", "2dGrpNonConnex"); + CPPUNIT_ASSERT(res.first); + CPPUNIT_ASSERT(res.second); +} + +MEDFileUMesh * CrackAlgoTest::make2x2Voxel( +) { + double coords[3] {0., 1., 2.}; + MCAuto x_coords = DataArrayDouble::New(); + x_coords->alloc(3); + copy(coords, coords + 3, x_coords->getPointer()); + + MCC m_cmesh = MEDCouplingCMesh::New(); + m_cmesh->setCoords(x_coords, x_coords, x_coords); + MCU m0 = m_cmesh->buildUnstructured(); + m0->setName("UMesh"); + + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU m1 = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + + MFU mf = MEDFileUMesh::New(); + mf->setCoords(m0->getCoords()); + mf->setMeshAtLevel(0, m0); + mf->setMeshAtLevel(-1, m1); + + return mf.retn(); +} + +MEDFileUMesh * CrackAlgoTest::make1DMesh( +) { + double coords[5] {0., 1., 2., 3., 4.}; + MCAuto x_coords = DataArrayDouble::New(); + x_coords->alloc(5); + copy(coords, coords + 5, x_coords->getPointer()); + + MCC m_cmesh = MEDCouplingCMesh::New(); + m_cmesh->setCoords(x_coords); + MCU m0 = m_cmesh->buildUnstructured(); + m0->setName("UMesh"); + + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU m1 = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + + MFU mf = MEDFileUMesh::New(); + mf->setCoords(m0->getCoords()); + mf->setMeshAtLevel(0, m0); + mf->setMeshAtLevel(-1, m1); + + return mf.retn(); +} + +MEDFileUMesh * CrackAlgoTest::make4x4Voxel( +) { + double coords[5] {0., 1., 2., 3., 4.}; + MCAuto x_coords = DataArrayDouble::New(); + x_coords->alloc(5); + copy(coords, coords + 5, x_coords->getPointer()); + + MCAuto m_cmesh = MEDCouplingCMesh::New(); + m_cmesh->setCoords(x_coords, x_coords, x_coords); + MCAuto m0 = m_cmesh->buildUnstructured(); + m0->setName("UMesh"); + + DAI desc(DataArrayIdType::New()), descIdx(DataArrayIdType::New()), revDesc(DataArrayIdType::New()), revDescIdx(DataArrayIdType::New()); + MCU m1 = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + + MEDFileUMesh * mf = MEDFileUMesh::New(); + mf->setCoords(m0->getCoords()); + mf->setMeshAtLevel(0, m0); + mf->setMeshAtLevel(-1, m1); + + return mf; +} + +MEDFileUMesh * CrackAlgoTest::make2DMesh( +) { + using DAD = MCAuto; + + double c0[6] {0.0, 1.1, 2.3, 3.6, 5.0, 6.5}; + DAD coords0 = DataArrayDouble::New(); + coords0->alloc(6); + copy(c0, c0 + 6, coords0->rwBegin()); + + double c1[5] {0.0, 1.1, 2.3, 3.6, 5.0}; + DAD coords1 = DataArrayDouble::New(); + coords1->alloc(5); + copy(c1, c1+ 5, coords1->rwBegin()); + + MCC m = MEDCouplingCMesh::New(); + m->setCoordsAt(0, coords0); + m->setCoordsAt(1, coords1); + + MCU m0 = m->buildUnstructured(); + m0->setName("duplicate"); + + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU m2 = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + m2->setName(m0->getName()); + + mcIdType ids_[17] {8,11,14,20,21,22,23,24,25,26,31,32,33,34,35,36,37}; + MCU m3 = m2->buildPartOfMySelf(ids_, ids_+17); + + mcIdType grp_[3] {4, 6, 8}; + DAI grp = DataArrayIdType::New(); + grp->alloc(3); + copy(grp_, grp_+3, grp->rwBegin()); + grp->setName("Grp"); + + mcIdType grp2_[2] {9, 16}; + DAI grp2 = DataArrayIdType::New(); + grp2->alloc(2); + copy(grp2_, grp2_+2, grp2->rwBegin()); + grp2->setName("Grp2"); + + MEDFileUMesh * mm = MEDFileUMesh::New(); + mm->setMeshAtLevel(0, m0); + mm->setMeshAtLevel(-1, m3); + mm->setGroupsAtLevel(-1, {grp, grp2}); + + mcIdType grpNode_[3] {4, 21, 23}; + DAI grpNode = DataArrayIdType::New(); + grpNode->alloc(3); + copy(grpNode_, grpNode_+3, grpNode->rwBegin()); + grpNode->setName("GrpNode"); + mm->setGroupsAtLevel(1, {grpNode}); + + return mm; +} + +MEDFileUMesh * CrackAlgoTest::make2DMesh2( +) { + using DAD = MCAuto; + + double c0[5] {0.0, 1.1, 2.3, 3.6, 5.0}; + DAD coords0 = DataArrayDouble::New(); + coords0->alloc(5); + copy(c0, c0 + 5, coords0->rwBegin()); + + double c1[3] {0.0, 1.0, 2.0}; + DAD coords1 = DataArrayDouble::New(); + coords1->alloc(3); + copy(c1, c1 + 3, coords1->rwBegin()); + + MCC m = MEDCouplingCMesh::New(); + m->setCoordsAt(0, coords0); + m->setCoordsAt(1, coords1); + + MCU m0 = m->buildUnstructured(); + m0->setName("simple"); + + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU m2 = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + m2->setName(m0->getName()); + + mcIdType grp_[2] {3, 19}; + DAI grp = DataArrayIdType::New(); + grp->alloc(2); + copy(grp_, grp_+2, grp->rwBegin()); + grp->setName("Grp"); + + MEDFileUMesh * mm = MEDFileUMesh::New(); + mm->setMeshAtLevel(0, m0); + mm->setMeshAtLevel(-1, m2); + mm->setGroupsAtLevel(-1, {grp}); + + return mm; +} + + +pair CrackAlgoTest::TestCrack( + MEDFileUMesh * mm_init, + const string & grp_name, + const string & test_name +) { + const auto c2cBroken = GetC2CBroken(mm_init, grp_name); + const auto c2cPreserved = GetC2CPreserved(mm_init, grp_name); + + const MCU f2dup = mm_init->getGroup(-1, grp_name); + + mm_init->write(test_name + "_in.med", 2); + const auto cellOld2NewNode = CrackAlgo::Compute(mm_init, grp_name); + mm_init->write(test_name + "_out.med", 2); + + const MCU f2dup_b = mm_init->getGroup(-1, grp_name); + + const auto res = pair { + CheckM0Mesh(mm_init, c2cBroken, c2cPreserved), + CheckM1Mesh(f2dup, f2dup_b) + }; + + CrackAlgo::OpenCrack(mm_init, cellOld2NewNode); + mm_init->write(test_name + "_cracked.med", 2); + + return res; +} + +bool CrackAlgoTest::CheckM1Mesh( + const MEDCouplingUMesh * f2dup_before, + const MEDCouplingUMesh * f2dup_after +) { + const mcIdType nbFaces_0 = f2dup_before->getNumberOfCells(); + const mcIdType nbFaces_1 = f2dup_after->getNumberOfCells(); + + MCU f2dup_before_copy = f2dup_after->deepCopy(); + f2dup_before_copy->zipCoords(); + + MCU f2dup_after_copy = f2dup_after->deepCopy(); + bool nodesAreMerged = false; + mcIdType newNbOfNodes = 0; + DAI o2n = f2dup_after_copy->mergeNodes(1e-12, nodesAreMerged, newNbOfNodes); + f2dup_after_copy->zipCoords(); + + DataArrayIdType * cellCor; + DataArrayIdType * nodeCor; + f2dup_after_copy->checkGeoEquivalWith(f2dup_before_copy, 12, 1e-12, cellCor, nodeCor); + DAI cellCorrI(cellCor); + DAI nodeCorrI(nodeCor); + + return (2 * nbFaces_0 == nbFaces_1); +} + +vector> +CrackAlgoTest::GetC2CBroken( + const MEDFileUMesh * mm, + const string & grp_name +) { + vector> res; + MCU m0 = mm->getMeshAtLevel(0); + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU mf = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + MCU f2dupMesh = mm->getGroup(-1, grp_name); + DataArrayIdType * f2dup_p; + mf->areCellsIncludedIn(f2dupMesh, 2, f2dup_p); + DAI f2dup(f2dup_p); + + for (const auto & f : *f2dup) { + const auto & cellsId_start = revDescIdx->begin()[f]; + const auto & cellsId_end = revDescIdx->begin()[f+1]; + if (cellsId_end - cellsId_start != 2) + throw INTERP_KERNEL::Exception( + "crackAlong: A face of group M1 does not have two neighbors."); + const auto & cell0 = revDesc->begin()[cellsId_start]; + const auto & cell1 = revDesc->begin()[cellsId_start + 1]; + res.push_back(pair(cell0, cell1)); + } + return res; +} + +vector> +CrackAlgoTest::GetC2CPreserved( + const MEDFileUMesh * mm, + const string & grp_name +) { + vector> res; + MCU m0 = mm->getMeshAtLevel(0); + DAI desc(DataArrayIdType::New()), + descIdx(DataArrayIdType::New()), + revDesc(DataArrayIdType::New()), + revDescIdx(DataArrayIdType::New()); + MCU mf = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + MCU f2dupMesh = mm->getGroup(-1, grp_name); + + DataArrayIdType * f2dup_p; + mf->areCellsIncludedIn(f2dupMesh, 2, f2dup_p); + DAI f2dup(f2dup_p); + + mcIdType nbOfFaces = mf->getNumberOfCells(); + for (mcIdType f = 0; f < nbOfFaces; f++) { + if (std::find(f2dup->begin(), f2dup->end(), f) != f2dup->end()) + continue; + + const auto & cellsId_start = revDescIdx->begin()[f]; + const auto & cellsId_end = revDescIdx->begin()[f+1]; + if (cellsId_end - cellsId_start == 1) + continue; + if (cellsId_end - cellsId_start > 2) + throw INTERP_KERNEL::Exception( + "crackAlong: A face of the mf has more than two neighbors."); + + const auto & cell0 = revDesc->begin()[cellsId_start]; + const auto & cell1 = revDesc->begin()[cellsId_start + 1]; + res.push_back(pair(cell0, cell1)); + } + return res; +} + +bool +CrackAlgoTest::CheckM0Mesh( + const MEDFileUMesh * mm, + const Connections& c2cBrokenConnection, + const Connections& c2cPreservedConnection +) { + MCU m0 = mm->getMeshAtLevel(0); + DAI desc(DataArrayIdType::New()), descIdx(DataArrayIdType::New()), revDesc(DataArrayIdType::New()), revDescIdx(DataArrayIdType::New()); + MCU mf = m0->buildDescendingConnectivity(desc, descIdx, revDesc, revDescIdx); + + CrackAlgo::Map2Set c2c; + + const mcIdType * descIdx_ptr = descIdx->begin(); + const mcIdType * desc_ptr = desc->begin(); + const mcIdType * revDescIdx_ptr = revDescIdx->begin(); + const mcIdType * revDesc_ptr = revDesc->begin(); + for (mcIdType cell = 0; cell < m0->getNumberOfCells(); cell++) { + auto & neighbors = c2c[cell]; + for ( + mcIdType descId = descIdx_ptr[cell]; + descId < descIdx_ptr[cell+1]; + descId++ + ) { + const auto & face = desc_ptr[descId]; + for ( + mcIdType revDescId = revDescIdx_ptr[face]; + revDescId < revDescIdx_ptr[face+1]; + revDescId++ + ) { + const auto & otherCell = revDesc_ptr[revDescId]; + if (otherCell != cell) + neighbors.insert(otherCell); + } + } + } + + bool res = true; + for (const auto & conn : c2cBrokenConnection) { + const auto & cell0 = conn.first; + const auto & cell1 = conn.second; + // NOTE: all cells are in c2c so cell0 and cell1 are in c2c + if (c2c.at(cell0).find(cell1) != c2c.at(cell0).end()) + res = false; + if (c2c.at(cell1).find(cell0) != c2c.at(cell1).end()) + res = false; + } + + for (const auto & conn : c2cPreservedConnection) { + const auto & cell0 = conn.first; + const auto & cell1 = conn.second; + // NOTE: all cells are in c2c so cell0 and cell1 are in c2c + if (c2c.at(cell0).find(cell1) == c2c.at(cell0).end()) + res = false; + if (c2c.at(cell1).find(cell0) == c2c.at(cell1).end()) + res = false; + } + return res; +} diff --git a/src/MEDLoader/Test/CrackAlgoTest.hxx b/src/MEDLoader/Test/CrackAlgoTest.hxx new file mode 100644 index 000000000..549942d03 --- /dev/null +++ b/src/MEDLoader/Test/CrackAlgoTest.hxx @@ -0,0 +1,108 @@ +// Copyright (C) 2007-2024 CEA, EDF +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author : Aymeric SONOLET (CEA/DES) + +#ifndef __CRACKALGOTEST_HXX__ +#define __CRACKALGOTEST_HXX__ + +#include +#include +#include +#include + +#include + +namespace MEDCoupling { + +class MEDFileUMesh; +class MEDCouplingUMesh; + +class CrackAlgoTest : public CppUnit::TestFixture { + CPPUNIT_TEST_SUITE(CrackAlgoTest); + CPPUNIT_TEST(test3DHalfCrossCut); + CPPUNIT_TEST(test3DFullCrossCut); + CPPUNIT_TEST(test3DHalfCut); + CPPUNIT_TEST(test3DFullCut); + CPPUNIT_TEST(test3DFullFullCut); + CPPUNIT_TEST(test3DAngleCut); + CPPUNIT_TEST(test2DGrp); + CPPUNIT_TEST(test2DMeshNonConnexCut); + CPPUNIT_TEST(test1DMesh); + CPPUNIT_TEST(testInnerCrossCut); + CPPUNIT_TEST_SUITE_END(); + + public: + + void test3DHalfCrossCut(); + void test3DFullCrossCut(); + void test3DAngleCut(); + void test3DFullCut(); + void test3DFullFullCut(); + void test3DHalfCut(); + void test2DGrp(); + void test2DMeshNonConnexCut(); + void test1DMesh(); + void testInnerCrossCut(); + + private: + using Connections = std::vector>; + static std::pair + TestCrack( + MEDFileUMesh * mm_init, + const std::string & grp_name, + const std::string & test_name); + + static bool + CheckM0Mesh( + const MEDFileUMesh * mm, + const Connections& c2cBrokenConnection, + const Connections& c2cPreservedConnection); + + static bool + CheckM1Mesh( + const MEDCouplingUMesh * f2dup_before, + const MEDCouplingUMesh * f2dup_after); + + static Connections + GetC2CBroken( + const MEDFileUMesh * mm, + const std::string & grp_name); + + static std::vector> + GetC2CPreserved( + const MEDFileUMesh * mm, + const std::string & grp_name); + + static MEDFileUMesh * + make2x2Voxel(); + + static MEDFileUMesh * + make4x4Voxel(); + + static MEDFileUMesh * + make2DMesh(); + + static MEDFileUMesh * + make2DMesh2(); + + static MEDFileUMesh * + make1DMesh(); +}; +} // namespace MEDCoupling +#endif // __CRACKALGOTEST_HXX__ diff --git a/src/MEDLoader/Test/MEDLoaderTest.cxx b/src/MEDLoader/Test/MEDLoaderTest.cxx index 0d1add2c1..5d1c77ddb 100644 --- a/src/MEDLoader/Test/MEDLoaderTest.cxx +++ b/src/MEDLoader/Test/MEDLoaderTest.cxx @@ -19,6 +19,7 @@ // Author : Anthony Geay (CEA/DEN) #include "MEDLoaderTest.hxx" +#include "MEDCouplingCMesh.hxx" #include "MEDLoader.hxx" #include "MEDLoaderBase.hxx" #include "MEDCouplingUMesh.hxx" @@ -27,8 +28,9 @@ #include "MEDCouplingFieldInt64.hxx" #include "MEDCouplingMemArray.hxx" #include "TestInterpKernelUtils.hxx" // getResourceFile() +#include "MEDFileMesh.hxx" -#include +#include #include using namespace MEDCoupling; @@ -1569,6 +1571,3 @@ MEDCouplingFieldDouble *MEDLoaderTest::buildVecFieldOnGaussNE_1() m->decrRef(); return f; } - - - diff --git a/src/MEDLoader/Test/MEDLoaderTest.hxx b/src/MEDLoader/Test/MEDLoaderTest.hxx index baccfc282..dc3332af2 100644 --- a/src/MEDLoader/Test/MEDLoaderTest.hxx +++ b/src/MEDLoader/Test/MEDLoaderTest.hxx @@ -61,6 +61,7 @@ namespace MEDCoupling CPPUNIT_TEST(testMEDLoaderPolygonRead); CPPUNIT_TEST(testMEDLoaderPolyhedronRead); + CPPUNIT_TEST_SUITE_END(); public: void testMesh1DRW(); @@ -89,6 +90,7 @@ namespace MEDCoupling void testMEDLoaderRead1(); void testMEDLoaderPolygonRead(); void testMEDLoaderPolyhedronRead(); + private: MEDCouplingUMesh *build1DMesh_1(); MEDCouplingUMesh *build2DCurveMesh_1(); diff --git a/src/MEDLoader/Test/TestMEDLoader.cxx b/src/MEDLoader/Test/TestMEDLoader.cxx index 68387f1f2..9e73982d1 100644 --- a/src/MEDLoader/Test/TestMEDLoader.cxx +++ b/src/MEDLoader/Test/TestMEDLoader.cxx @@ -20,7 +20,9 @@ #include "CppUnitTest.hxx" #include "MEDLoaderTest.hxx" +#include "CrackAlgoTest.hxx" CPPUNIT_TEST_SUITE_REGISTRATION( MEDCoupling::MEDLoaderTest ); +CPPUNIT_TEST_SUITE_REGISTRATION( MEDCoupling::CrackAlgoTest ); #include "BasicMainTest.hxx" -- 2.39.2