X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh_intersection.cxx;h=d1a9a974c0199350143827d9cbcfedf0a7e1424e;hb=e7a9d4f59978fd384ee98db1dfdd5ec2118331ca;hp=e72c9b5e2f1b909f27a4b123f361f38432127ae8;hpb=e0c843a1fe827a90af91ada8d2033ffb3a7dd1d8;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index e72c9b5e2..d1a9a974c 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -556,7 +556,15 @@ void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MED } } -MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector >& intersectEdge2, const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, +/* + * Build the final 1D mesh resulting from the newly created points after intersection with the segments of the descending 2D mesh. + * @param[out] idsInRetColinear IDs of edges in the result (ret) that are colinears to one of the segment of the descending 2D mesh. Indexing scheme + * is the one of the ret 1D mesh. + * @param[out] idsInMesh1DForIdsInRetColinear same IDs as above in the descending 2D mesh + */ +MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector >& intersectEdge2, + const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, + const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, MCAuto& idsInRetColinear, MCAuto& idsInMesh1DForIdsInRetColinear) { idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1); @@ -581,7 +589,8 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: mcIdType nbSubEdge=ToIdType(subEdges.size()/2); for(mcIdType j=0;j n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); + MCAuto n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)), + n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); MCAuto e2(e->buildEdgeLyingOnMe(n1,n2)); INTERP_KERNEL::Edge *e2Ptr(e2); std::map::const_iterator itm; @@ -619,6 +628,7 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: } e->decrRef(); } + MCAuto ret(MEDCouplingUMesh::New(mesh1D->getName(),1)); ret->setConnectivity(cOut,ciOut,true); MCAuto arr3(DataArrayDouble::New()); @@ -1372,8 +1382,10 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the * (newly created) nodes corresponding to the edge intersections. * Output params: - * @param[out] cr, crI connectivity of the resulting mesh - * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 + * @param[out] cr connectivity of the resulting mesh + * @param[out] crI connectivity of the resulting mesh + * @param[out] cNb1 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 + * @param[out] cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 * TODO: describe input parameters */ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const mcIdType *desc1, const mcIdType *descIndx1, @@ -1477,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]); @@ -1670,6 +1682,9 @@ mcIdType MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayIdType *desc, co * The meshes should be in 2D space. In * addition, returns two arrays mapping cells of the result mesh to cells of the input * meshes. + * \b WARNING: the two meshes should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and + * orientCorrectly2DCells() can be used for this. + * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells) * \param [in] m1 - the first input mesh which is a partitioned object. The mesh must be so that each point in the space covered by \a m1 * must be covered exactly by one entity, \b no \b more. If it is not the case, some tools are available to heal the mesh (conformize2D, mergeNodes) * \param [in] m2 - the second input mesh which is a partition tool. The mesh must be so that each point in the space covered by \a m2 @@ -1753,6 +1768,10 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 * and finally, in case of quadratic polygon the centers of edges new nodes. * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input. * + * \b WARNING: the 2D mesh should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and + * orientCorrectly2DCells() can be used for this. + * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells) + * * \param [in] mesh2D - the 2D mesh (spacedim=meshdim=2) to be intersected using \a mesh1D tool. The mesh must be so that each point in the space covered by \a mesh2D * must be covered exactly by one entity, \b no \b more. If it is not the case, some tools are available to heal the mesh (conformize2D, mergeNodes) * \param [in] mesh1D - the 1D mesh (spacedim=2 meshdim=1) the is the tool that will be used to intersect \a mesh2D. \a mesh1D must be ordered consecutively. If it is not the case @@ -1812,13 +1831,57 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MCAuto idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear; MCAuto ret2(DataArrayIdType::New()); ret2->alloc(0,1); MCAuto ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1, - idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); + idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); + + // ### Colinearity fix : + // if a node in ret1 has been merged with an already existing node in mesh2D, + // we might end up with edges in ret1 that are colinear with some edges in mesh2D + // even if none of the edges of the two original meshes were colinear. + // this procedure detects such edges and adds them in idsInRet1Colinear/idsInDescMesh2DForIdsInRetColinear + // a- find edges in ret1 that are in contact with 2 cells + MCAuto centerOfMassRet1(ret1->computeCellCenterOfMass()); + MCAuto cells,cellsIndex; + mesh2D->getCellsContainingPoints(centerOfMassRet1->begin(),centerOfMassRet1->getNumberOfTuples(),eps,cells,cellsIndex); + MCAuto cellsIndex2(DataArrayIdType::New()); cellsIndex2->alloc(0,1); + if (cellsIndex->getNumberOfTuples() > 1) + cellsIndex2 = cellsIndex->deltaShiftIndex(); + MCAuto idsInRet1With2Contacts(cellsIndex2->findIdsEqual(2)); + + MCAuto realIdsInDesc2D(desc1->deepCopy()); + realIdsInDesc2D->abs(); realIdsInDesc2D->applyLin(1,-1); + const mcIdType *cRet1(ret1->getNodalConnectivity()->begin()),*ciRet1(ret1->getNodalConnectivityIndex()->begin()); + for(const mcIdType *it=idsInRet1With2Contacts->begin();it!=idsInRet1With2Contacts->end();it++) + { + // b- find the edge that the 2 cells in m1Desc have in common: + // this is the edge which is colinear with the one in ret1 + const mcIdType* cellId1 = cells->begin() + cellsIndex->begin()[*it]; + const mcIdType* cellId2 = cells->begin() + cellsIndex->begin()[*it]+1; + + 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]); + + 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; + mcIdType commonEdgeIdCorrectlyOriented; + if(IsColinearOfACellOf(intersectEdge1, commonEdgeId, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1, commonEdgeIdCorrectlyOriented)) + { + idsInRet1Colinear->pushBackSilent(*it); + idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeIdCorrectlyOriented); + } + } + // ### End colinearity fix + MCAuto ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits::max()); ret3->rearrange(2); MCAuto idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells())); // deal with cells in mesh2D that are not cut but only some of their edges are MCAuto idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy()); idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1); idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique(); + MCAuto out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells if(!idsInDesc2DToBeRefined->empty()) { @@ -1831,9 +1894,11 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, } // MCAuto ret1NonCol(static_cast(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end()))); - MCAuto baryRet1(ret1NonCol->computeCellCenterOfMass()); - MCAuto elts,eltsIndex; - mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex); + MCAuto baryRet1(centerOfMassRet1->selectByTupleId(idsInRet1NotColinear->begin(), idsInRet1NotColinear->end())); + DataArrayIdType *out(0),*outi(0); + DataArrayIdType::ExtractFromIndexedArrays(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end(),cells,cellsIndex,out,outi); + MCAuto elts(out),eltsIndex(outi); + MCAuto eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1); if (eltsIndex->getNumberOfTuples() > 1) eltsIndex2 = eltsIndex->deltaShiftIndex();