Salome HOME
2D PointLocator remapping now works properly on non-convex polygons.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh.cxx
index a890398e76df2b7f2822927b956d4c1d0cc6c2a0..c20a00b273257223bf222b714f56ef62e6871086 100755 (executable)
@@ -618,7 +618,7 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayIdType *revNodal, Da
       const mcIdType *endNdlConnOfCurCell=conn+connIndex[eltId+1];
       for(const mcIdType *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
         if(*iter>=0)//for polyhedrons
-          *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind2nd(std::equal_to<mcIdType>(),-1))=eltId;
+          *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind(std::equal_to<mcIdType>(),std::placeholders::_1,-1))=eltId;
     }
 }
 
@@ -1733,7 +1733,7 @@ void MEDCouplingUMesh::FindCommonCellsAlg(int compType, mcIdType startCellId, co
         {
           if(!isFetched[i])
             {
-              const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1));
+              const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1));
               std::vector<mcIdType> v,v2;
               if(connOfNode!=connPtr+connIPtr[i+1])
                 {
@@ -1768,7 +1768,7 @@ void MEDCouplingUMesh::FindCommonCellsAlg(int compType, mcIdType startCellId, co
         {
           if(!isFetched[i])
             {
-              const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1));
+              const mcIdType *connOfNode=std::find_if(connPtr+connIPtr[i]+1,connPtr+connIPtr[i+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1));
               // v2 contains the result of successive intersections using rev nodal on on each node of cell #i
               std::vector<mcIdType> v,v2;
               if(connOfNode!=connPtr+connIPtr[i+1])
@@ -2322,7 +2322,6 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const
  * \param [out] cellIdsNeededToBeRenum cell ids in \b this in which the renumber of nodes should be performed.
  * \param [out] cellIdsNotModified cell ids mcIdType \b this that lies on \b otherDimM1OnSameCoords mesh whose connectivity do \b not need to be modified as it is the case for \b cellIdsNeededToBeRenum.
  *
- * \warning This method modifies param \b otherDimM1OnSameCoords (for speed reasons).
  */
 void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayIdType *& nodeIdsToDuplicate,
                                             DataArrayIdType *& cellIdsNeededToBeRenum, DataArrayIdType *& cellIdsNotModified) const
@@ -2339,14 +2338,15 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
 
   // Checking star-shaped M1 group:
   DAInt dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New();
-  MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0);
+  MCUMesh meshM2 = otherDimM1OnSameCoords.buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); // 2D: a mesh of points, 3D: a mesh of segs
   DAInt dsi = rdit0->deltaShiftIndex();
-  DAInt idsTmp0 = dsi->findIdsNotInRange(-1, 3);
+  DAInt idsTmp0 = dsi->findIdsNotInRange(-1, 3);  // for 2D: if a point is connected to more than 2 segs. For 3D: if a seg is connected to more than two faces.
   if(idsTmp0->getNumberOfTuples())
     throw INTERP_KERNEL::Exception("MEDFileUMesh::buildInnerBoundaryAlongM1Group: group is too complex: some points (or edges) have more than two connected segments (or faces)!");
   dt0=0; dit0=0; rdt0=0; rdit0=0; idsTmp0=0;
 
-  // Get extreme nodes from the group (they won't be duplicated), ie nodes belonging to boundary cells of M1
+  // Get extreme nodes from the group (they won't be duplicated except if they also lie on bound of M0 -- see below),
+  // ie nodes belonging to the boundary "cells" (might be points) of M1
   DAInt xtremIdsM2 = dsi->findIdsEqual(1); dsi = 0;
   MCUMesh meshM2Part = static_cast<MEDCouplingUMesh *>(meshM2->buildPartOfMySelf(xtremIdsM2->begin(), xtremIdsM2->end(),true));
   DAInt xtrem = meshM2Part->computeFetchedNodeIds();
@@ -2354,15 +2354,17 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
   dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New();
   MCUMesh m0desc = buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0;
   dsi = rdit0->deltaShiftIndex();
-  DAInt boundSegs = dsi->findIdsEqual(1);   // boundary segs/faces of the M0 mesh
+  DAInt boundSegs = dsi->findIdsEqual(1);  dsi = 0; // boundary segs/faces of the M0 mesh
   MCUMesh m0descSkin = static_cast<MEDCouplingUMesh *>(m0desc->buildPartOfMySelf(boundSegs->begin(),boundSegs->end(), true));
   DAInt fNodes = m0descSkin->computeFetchedNodeIds();
-  // In 3D, some points on the boundary of M0 still need duplication:
+  // In 3D, some points on the boundary of M0 will NOT be duplicated (where as in 2D, points on the boundary of M0 are always duplicated)
+  // Think of a partial (plane) crack in a cube: the points at the tip of the crack and not located inside the volume of the cube are not duplicated
+  // although they are technically on the skin of the cube.
   DAInt notDup = 0;
   if (getMeshDimension() == 3)
     {
       DAInt dnu1=DataArrayIdType::New(), dnu2=DataArrayIdType::New(), dnu3=DataArrayIdType::New(), dnu4=DataArrayIdType::New();
-      MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4);
+      MCUMesh m0descSkinDesc = m0descSkin->buildDescendingConnectivity(dnu1, dnu2, dnu3, dnu4); // all segments of the skin of the 3D (M0) mesh
       dnu1=0;dnu2=0;dnu3=0;dnu4=0;
       DataArrayIdType * corresp=0;
       meshM2->areCellsIncludedIn(m0descSkinDesc,2,corresp);
@@ -2370,10 +2372,26 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
       corresp->decrRef();
       if (validIds->getNumberOfTuples())
         {
+          // Build the set of segments which are: in the desc mesh of the skin of the 3D mesh (M0) **and** in the desc mesh of the M1 group:
           MCUMesh m1IntersecSkin = static_cast<MEDCouplingUMesh *>(m0descSkinDesc->buildPartOfMySelf(validIds->begin(), validIds->end(), true));
+          // Its boundary nodes should no be duplicated (this is for example the tip of the crack inside the cube described above)
           DAInt notDuplSkin = m1IntersecSkin->findBoundaryNodes();
           DAInt fNodes1 = fNodes->buildSubstraction(notDuplSkin);
-          notDup = xtrem->buildSubstraction(fNodes1);
+
+          // Also, in this (segment) mesh, nodes connected to more than 3 segs should not be dup either (singular points - see testBuildInnerBoundary6())
+          dt0=DataArrayIdType::New(),dit0=DataArrayIdType::New(),rdt0=DataArrayIdType::New(),rdit0=DataArrayIdType::New();
+          MCUMesh meshM2Desc = meshM2->buildDescendingConnectivity(dt0, dit0, rdt0, rdit0); dt0=0; dit0=0; rdt0=0;  // a mesh made of node cells
+          dsi = rdit0->deltaShiftIndex();
+          DAInt singPoints = dsi->findIdsNotInRange(-1,4);   // points connected to (strictly) more than 3 segments
+          const mcIdType *cc = meshM2Desc->getNodalConnectivity()->begin(), *ccI = meshM2Desc->getNodalConnectivityIndex()->begin();
+          mcIdType * singPointsP = singPoints->rwBegin();
+          for (mcIdType j=0; j < singPoints->getNumberOfTuples(); j++) // replace ids in singPoints by real coordinate index (was index of cells in notDuplSkin)
+            {
+              mcIdType nodeCellIdx = singPointsP[j];
+              singPointsP[j] = cc[ccI[nodeCellIdx]+1];  // +1 to skip type
+            }
+          DAInt fNodes2 = fNodes1->buildSubstraction(singPoints);
+          notDup = xtrem->buildSubstraction(fNodes2);
         }
       else
         notDup = xtrem->buildSubstraction(fNodes);
@@ -3014,7 +3032,7 @@ mcIdType MEDCouplingUMesh::getNumberOfNodesInCell(mcIdType cellId) const
   if(pt[ptI[cellId]]!=INTERP_KERNEL::NORM_POLYHED)
     return ptI[cellId+1]-ptI[cellId]-1;
   else
-    return ToIdType(std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind2nd(std::not_equal_to<mcIdType>(),-1)));
+    return ToIdType(std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind(std::not_equal_to<mcIdType>(),std::placeholders::_1,-1)));
 }
 
 /*!
@@ -3251,7 +3269,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const
           area_vol[iel]=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[iel+1]-ipt-1,coords,dim_space);
         }
       if(isAbs)
-        std::transform(area_vol,area_vol+nbelem,area_vol,std::ptr_fun<double,double>(fabs));
+        std::transform(area_vol,area_vol+nbelem,area_vol,[](double c){return fabs(c);});
     }
   else
     {
@@ -3304,7 +3322,7 @@ DataArrayDouble *MEDCouplingUMesh::getPartMeasureField(bool isAbs, const mcIdTyp
           *area_vol++=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[*iel+1]-ipt-1,coords,dim_space);
         }
       if(isAbs)
-        std::transform(array->getPointer(),area_vol,array->getPointer(),std::ptr_fun<double,double>(fabs));
+        std::transform(array->getPointer(),area_vol,array->getPointer(),[](double c){return fabs(c);});
     }
   else
     {
@@ -3401,7 +3419,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
               mcIdType offset=connI[i];
               INTERP_KERNEL::crossprod<3>(locPtr+3*i,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
               double n=INTERP_KERNEL::norm<3>(vals);
-              std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
+              std::transform(vals,vals+3,vals,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
             }
         }
       else
@@ -3420,7 +3438,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
           mcIdType offset=connI[i];
           std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
           double n=INTERP_KERNEL::norm<2>(tmp);
-          std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
+          std::transform(tmp,tmp+2,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
           *vals++=-tmp[1];
           *vals++=tmp[0];
         }
@@ -3480,7 +3498,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildPartOrthogonalField(const mcIdTyp
               mcIdType offset=connI[*i];
               INTERP_KERNEL::crossprod<3>(locPtr,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
               double n=INTERP_KERNEL::norm<3>(vals);
-              std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
+              std::transform(vals,vals+3,vals,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
             }
         }
       else
@@ -3497,7 +3515,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildPartOrthogonalField(const mcIdTyp
           mcIdType offset=connI[*i];
           std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
           double n=INTERP_KERNEL::norm<2>(tmp);
-          std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
+          std::transform(tmp,tmp+2,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1./n));
           *vals++=-tmp[1];
           *vals++=tmp[0];
         }
@@ -6412,7 +6430,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::keepSpecifiedCells(INTERP_KERNEL::Normalized
         for(mcIdType j=0;j<code[3*i+1];j++)
           *idsPtr++=offset+j;
       else
-        idsPtr=std::transform(idsPerGeoTypeBg,idsPerGeoTypeEnd,idsPtr,std::bind2nd(std::plus<mcIdType>(),offset));
+        idsPtr=std::transform(idsPerGeoTypeBg,idsPerGeoTypeEnd,idsPtr,std::bind(std::plus<mcIdType>(),std::placeholders::_1,offset));
       offset+=code[3*i+1];
     }
   MCAuto<MEDCouplingUMesh> ret=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(idsTokeep->begin(),idsTokeep->end(),true));
@@ -6547,7 +6565,7 @@ DataArrayDouble *MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell() const
             }
           mcIdType nbOfNodesInCell=nodalI[i+1]-nodalI[i]-1;
           if(nbOfNodesInCell>0)
-            std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+            std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)nbOfNodesInCell));
           else
             {
               std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on cell #" << i << " presence of cell with no nodes !";
@@ -6569,7 +6587,7 @@ DataArrayDouble *MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell() const
                 }
             }
           if(!s.empty())
-            std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind2nd(std::multiplies<double>(),1./(double)s.size()));
+            std::transform(ptToFill,ptToFill+spaceDim,ptToFill,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)s.size()));
           else
             {
               std::ostringstream oss; oss << "MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell : on polyhedron cell #" << i << " there are no nodes !";
@@ -6684,7 +6702,7 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
               for(mcIdType offset=nodalI[0]+1;offset<nodalI[1];offset++)
                 std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
               mcIdType nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
-              std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+              std::transform(dd,dd+3,dd,std::bind(std::multiplies<double>(),std::placeholders::_1,1./(double)nbOfNodesInCell));
               std::copy(dd,dd+3,matrix+4*2);
               INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
               retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
@@ -6893,7 +6911,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesOnSameCoords(const std::vector<c
       mcIdType meshLgth2=(*iter)->getNodalConnectivityArrayLen();
       nodalPtr=std::copy(nod,nod+meshLgth2,nodalPtr);
       if(iter!=meshes.begin())
-        nodalIndexPtr=std::transform(index+1,index+nbOfCells+1,nodalIndexPtr,std::bind2nd(std::plus<mcIdType>(),offset));
+        nodalIndexPtr=std::transform(index+1,index+nbOfCells+1,nodalIndexPtr,std::bind(std::plus<mcIdType>(),std::placeholders::_1,offset));
       else
         nodalIndexPtr=std::copy(index,index+nbOfCells+1,nodalIndexPtr);
       offset+=meshLgth2;