[EDF27988] : Implementation of MEDCouplingUMesh.explodeMeshTo for MEDFileUMesh.reduce...

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh_internal.cxx

diff --git a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx
index 9bf681736475ebd16957fc2e34381e25cd8ad03a..4acfeb5601cd7c7c3f6a835da791c3d8df95d523 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingUMesh_internal.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh_internal.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2023  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -547,7 +547,7 @@ void MEDCouplingUMesh::fillCellIdsToKeepFromNodeIds(const mcIdType *begin, const
  * This method has one in/out parameter : 'cut3DCurve'.
  * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
  * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
- * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
+ * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully (or partially) part of plane defined by ('origin','vec').
  * This method will throw an exception if \a this contains a non linear segment.
  */
 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<mcIdType>& cut3DCurve)
@@ -1313,6 +1313,74 @@ DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCoupling
   return ret.retn();
 }
 
+/*!
+ * \param [in] arrIn array part of the indexed array pair to be treated
+ * \param [in] arrIndxIn array index part of the indexed array pair to be treated
+ * \param [out] arrOut array part of the indexed array pair containing result
+ * \param [out] arrIndxOut array index part of the indexed array pair containing result
+ */
+void MEDCouplingUMesh::DeleteCellTypeInIndexedArray(const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, MCAuto<DataArrayIdType>& arrOut, MCAuto<DataArrayIdType>& arrIndxOut)
+{
+  if( !arrIn || !arrIn->isAllocated() )
+    THROW_IK_EXCEPTION("input array is null or not allocated !");
+  if( !arrIndxIn || !arrIndxIn->isAllocated() )
+    THROW_IK_EXCEPTION("input indexed array is null or not allocated !");
+  arrIn->checkNbOfComps(1,"input array"); arrIndxIn->checkNbOfComps(1,"input indexed array");
+  mcIdType arrNbTuples(arrIn->getNumberOfTuples()),arrIndxNbTuples(arrIndxIn->getNumberOfTuples());
+  if( arrIndxNbTuples < 1 )
+    THROW_IK_EXCEPTION("Indexed array is supposed to have length >= 1 !");
+  if( arrNbTuples < arrIndxNbTuples )
+    THROW_IK_EXCEPTION("Number of tuples of input array is to low compared to indexed array one !");
+  const mcIdType *inArrPtr(arrIn->begin()),*inArrIndxPtr(arrIndxIn->begin());
+  if( *inArrIndxPtr != 0 )
+    THROW_IK_EXCEPTION("First value of indexed array must be 0 !");
+  arrOut = DataArrayIdType::New(); arrOut->alloc(arrNbTuples-arrIndxNbTuples+1,1);
+  arrIndxOut = DataArrayIdType::New(); arrIndxOut->alloc(arrIndxNbTuples,1);
+  bool presenceOfPolyh = false;
+  {
+    mcIdType *outArrPtr(arrOut->getPointer()),*outArrIndxPtr(arrIndxOut->getPointer());
+    *outArrIndxPtr++ = 0;
+    for( mcIdType i = 0 ; i < arrIndxNbTuples - 1 ; ++i )
+    {
+      mcIdType startPos(*inArrIndxPtr++);
+      mcIdType endPos(*inArrIndxPtr);
+      if(inArrPtr[startPos] == INTERP_KERNEL::NORM_POLYHED)
+        {
+          presenceOfPolyh = true;
+          break;
+        }
+      outArrPtr = std::copy(inArrPtr+startPos+1,inArrPtr+endPos,outArrPtr);
+      *outArrIndxPtr++ = endPos - i - 1;
+    }
+  }
+  if(!presenceOfPolyh)
+    return ;
+  //
+  {
+    arrOut = DataArrayIdType::New(); arrOut->alloc(0,1);
+    inArrIndxPtr = arrIndxIn->begin();
+    mcIdType *outArrIndxPtr = arrIndxOut->getPointer();
+    *outArrIndxPtr++ = 0;
+    for( mcIdType i = 0 ; i < arrIndxNbTuples - 1 ; ++i,++outArrIndxPtr )
+    {
+      mcIdType startPos(*inArrIndxPtr++);
+      mcIdType endPos(*inArrIndxPtr);
+      if(inArrPtr[startPos] != INTERP_KERNEL::NORM_POLYHED)
+      {
+        arrOut->insertAtTheEnd(inArrPtr+startPos+1,inArrPtr+endPos);
+        outArrIndxPtr[0] = outArrIndxPtr[-1] + (endPos-startPos-1);
+      }
+      else
+      {
+        std::set<mcIdType> s(inArrPtr+startPos+1,inArrPtr+endPos);
+        s.erase( -1 );
+        arrOut->insertAtTheEnd(s.begin(),s.end());
+        outArrIndxPtr[0] = outArrIndxPtr[-1] + s.size();
+      }
+    }
+  }
+}
+
 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
 {
   if(a.empty())
@@ -1624,11 +1692,20 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>&
                 res.push_back(status);
               else
                 {
-                  res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
-                  res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
+                  const mcIdType& node1 = nodal3DCurve[nodalIndx3DCurve[edgeId]+1];
+                  const mcIdType& node2 = nodal3DCurve[nodalIndx3DCurve[edgeId]+2];
+                  // Here, we have an edge that has either one or both of its nodes intersecting the plane
+                  // we're only adding the nodes from the edges fully contained in the plane
+                  if(std::find(nodesOnP.begin(), nodesOnP.end(),node1) != nodesOnP.end()
+                  && std::find(nodesOnP.begin(), nodesOnP.end(),node2) != nodesOnP.end())
+                    {
+                      res.push_back(node1);
+                      res.push_back(node2);
+                    }
                 }
             }
         }
+
       switch(res.size())
       {
         case 2: