X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh_intersection.cxx;h=f938bb103dd89d0f5b34d36151a7d98d6c0e6221;hb=0cc7b8ec634eb17df8b0b70c26e7473378e16a5b;hp=b480bec30b67740d31480fafa6b1662deac5ea0b;hpb=40904f8238992f0f20e02194b0816b03df977690;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index b480bec30..f938bb103 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.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 @@ -1127,7 +1127,7 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh if (i < nbCellsInSplitMesh1D-1) { // Build circular permutation to shift consecutive edges together - renumb->iota(i+1); + renumb->iota(nbCellsInSplitMesh1D-1-i); renumb->applyModulus(nbCellsInSplitMesh1D); splitMesh1D->renumberCells(renumbP, false); cSplitPtr = splitMesh1D->getNodalConnectivity()->begin(); @@ -1489,8 +1489,8 @@ void InsertNodeInConnIfNecessary(mcIdType nodeIdToInsert, std::vector& mcIdType pt0(conn[i]),pt1(conn[(i+1)%sz]); double v1[3]={coords[3*pt1+0]-coords[3*pt0+0],coords[3*pt1+1]-coords[3*pt0+1],coords[3*pt1+2]-coords[3*pt0+2]},v2[3]={coords[3*nodeIdToInsert+0]-coords[3*pt0+0],coords[3*nodeIdToInsert+1]-coords[3*pt0+1],coords[3*nodeIdToInsert+2]-coords[3*pt0+2]}; double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])); - std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies(),1./normm)); - std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies(),1./normm)); + std::transform(v1,v1+3,v1,std::bind(std::multiplies(),std::placeholders::_1,1./normm)); + std::transform(v2,v2+3,v2,std::bind(std::multiplies(),std::placeholders::_1,1./normm)); double v3[3]; v3[0]=v1[1]*v2[2]-v1[2]*v2[1]; v3[1]=v1[2]*v2[0]-v1[0]*v2[2]; v3[2]=v1[0]*v2[1]-v1[1]*v2[0]; double normm2(sqrt(v3[0]*v3[0]+v3[1]*v3[1]+v3[2]*v3[2])),dotTest(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]); @@ -1859,17 +1859,18 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, std::set s1(realIdsInDesc2D->begin()+dd2->begin()[*cellId1], realIdsInDesc2D->begin()+dd2->begin()[*cellId1+1]); std::set s2(realIdsInDesc2D->begin()+dd2->begin()[*cellId2], realIdsInDesc2D->begin()+dd2->begin()[*cellId2+1]); - mcIdType commonEdgeId; - std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(), &commonEdgeId); + + std::vector commonEdgeId; + std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(), std::back_inserter(commonEdgeId)); // c- find correct orientation for commonEdgeId const mcIdType* firstNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+1; const mcIdType* secondNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+2; - std::vector v; v.push_back(commonEdgeId); - if(IsColinearOfACellOf(intersectEdge1, v, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1,commonEdgeId)) + mcIdType commonEdgeIdCorrectlyOriented; + if(IsColinearOfACellOf(intersectEdge1, commonEdgeId, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1, commonEdgeIdCorrectlyOriented)) { idsInRet1Colinear->pushBackSilent(*it); - idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeId); + idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeIdCorrectlyOriented); } } // ### End colinearity fix @@ -2335,13 +2336,17 @@ void MEDCouplingUMesh::ReplaceEdgeInFace(const mcIdType * sIdxConn, const mcIdTy * This method expects that all nodes in \a this are not closer than \a eps. * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method. * - * \param [in] eps the relative error to detect merged edges. + * \b WARNING this method only works for 'partition-like' non-conformities. When joining adjacent faces, the set of all small faces must fit exactly into the + * neighboring bigger face. No real face intersection is actually computed. + * + * \param [in] eps the absolute (geometric) error to detect merged edges. * \return DataArrayIdType * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array * that the user is expected to deal with. * * \throw If \a this is not coherent. * \throw If \a this has not spaceDim equal to 3. * \throw If \a this has not meshDim equal to 3. + * \throw If the 'partition-like' condition described above is not respected. * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly() */ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) @@ -2358,6 +2363,7 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); const double * coo(_coords->begin()); MCAuto ret(DataArrayIdType::New()); ret->alloc(0,1); + const double carMeshSz = getCaracteristicDimension(); { /************************* @@ -2404,14 +2410,14 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) // Keep only candidates whose normal matches the normal of current face for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) { - bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps); + bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps/carMeshSz); // using rough relative epsilon if (*it2 != faceIdx && col) cands2.push_back(*it2); } if (!cands2.size()) continue; - // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later + // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes one day INTERP_KERNEL::TranslationRotationMatrix rotation; INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]), coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]), @@ -2574,13 +2580,17 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3]; + double carSz = INTERP_KERNEL::caracteristicDimVector(vCurr), eps2 = eps*carSz*carSz*carSz; for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) { double vOther[3]; mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; for (mcIdType i3=0; i3 < 3; i3++) vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3]; - bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps); + // isColinear() expect unit vecotr, and relative (non geometrical) precision. + // relative prec means going from eps -> eps/curSz + // normalizing vCurr and vOther means multiplying by v^4, knowing that inside isColinear3D() the square (^2) of a dot product (^2) is computed + bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps2); // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine()) if (col) @@ -2595,7 +2605,7 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) INTERP_KERNEL::TranslationRotationMatrix rotation; INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation); MCAuto mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called - MCAuto mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called + MCAuto mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is NOT called MCAuto nodeMap(mPartCand->zipCoordsTraducer()); mcIdType nbElemsNotM1; { @@ -2610,8 +2620,7 @@ DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) rotation.transform_vector(cooPartCand+SPACEDIM*ii); - - // Eliminate all edges for which y or z is not null + // Eliminate all edges for which y or z is not null - those are colinear edges which are simply parallels, but not on the same line MCAuto baryPart = mPartCand->computeCellCenterOfMass(); vector compo; compo.push_back(1); MCAuto baryPartY = baryPart->keepSelectedComponents(compo);