-// Copyright (C) 2007-2019 CEA/DEN, EDF R&D
+// 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
}
}
-MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<mcIdType> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<mcIdType,mcIdType>& mergedNodes, const std::vector< std::vector<mcIdType> >& colinear2, const std::vector< std::vector<mcIdType> >& 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<mcIdType> >& intersectEdge2,
+ const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<mcIdType,mcIdType>& mergedNodes,
+ const std::vector< std::vector<mcIdType> >& colinear2, const std::vector< std::vector<mcIdType> >& intersectEdge1,
MCAuto<DataArrayIdType>& idsInRetColinear, MCAuto<DataArrayIdType>& idsInMesh1DForIdsInRetColinear)
{
idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1);
mcIdType nbSubEdge=ToIdType(subEdges.size()/2);
for(mcIdType j=0;j<nbSubEdge;j++,kk++)
{
- MCAuto<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
+ MCAuto<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),
+ n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
INTERP_KERNEL::Edge *e2Ptr(e2);
std::map<mcIdType,mcIdType>::const_iterator itm;
}
e->decrRef();
}
+
MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
ret->setConnectivity(cOut,ciOut,true);
MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
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();
* 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,
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<double>(),1./normm));
- std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
+ std::transform(v1,v1+3,v1,std::bind(std::multiplies<double>(),std::placeholders::_1,1./normm));
+ std::transform(v2,v2+3,v2,std::bind(std::multiplies<double>(),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]);
* 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
* 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
{
if(!mesh2D || !mesh1D)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
+
mesh2D->checkFullyDefined();
mesh1D->checkFullyDefined();
+
+ { // (In a scope to free mem early)
+ // Check if any node within mesh1D is connected to more than two edges
+ DataArrayIdType *mesh1DDesc1(DataArrayIdType::New()),*mesh1DDescIndx1(DataArrayIdType::New()),*mesh1DRevDesc1(DataArrayIdType::New()),*mesh1DRevDescIndx1(DataArrayIdType::New());
+ MCAuto<DataArrayIdType> m1dd1(mesh1DDesc1),m1dd2(mesh1DDescIndx1),m1dd3(mesh1DRevDesc1),m1dd4(mesh1DRevDescIndx1);
+ MCAuto<MEDCouplingUMesh> mesh1Desc(mesh1D->buildDescendingConnectivity(mesh1DDesc1,mesh1DDescIndx1,mesh1DRevDesc1,mesh1DRevDescIndx1));
+ MCAuto<DataArrayIdType> dsi=mesh1DRevDescIndx1->deltaShiftIndex();
+ MCAuto<DataArrayIdType> tpls = dsi->findIdsGreaterOrEqualTo(3);
+ if(tpls->getNumberOfTuples() != 0)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : Certain nodes of the 1D mesh are connected to more than 2 edges !\
+ If the 1D mesh consists of mutiple lines, make sure those lines are not overlapping. ");
+ }
+
const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
MCAuto<DataArrayIdType> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
MCAuto<DataArrayIdType> ret2(DataArrayIdType::New()); ret2->alloc(0,1);
MCAuto<MEDCouplingUMesh> 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<DataArrayDouble> centerOfMassRet1(ret1->computeCellCenterOfMass());
+ MCAuto<DataArrayIdType> cells,cellsIndex;
+ mesh2D->getCellsContainingPoints(centerOfMassRet1->begin(),centerOfMassRet1->getNumberOfTuples(),eps,cells,cellsIndex);
+ MCAuto<DataArrayIdType> cellsIndex2(DataArrayIdType::New()); cellsIndex2->alloc(0,1);
+ if (cellsIndex->getNumberOfTuples() > 1)
+ cellsIndex2 = cellsIndex->deltaShiftIndex();
+ MCAuto<DataArrayIdType> idsInRet1With2Contacts(cellsIndex2->findIdsEqual(2));
+
+ MCAuto<DataArrayIdType> 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<mcIdType> s1(realIdsInDesc2D->begin()+dd2->begin()[*cellId1], realIdsInDesc2D->begin()+dd2->begin()[*cellId1+1]);
+ std::set<mcIdType> s2(realIdsInDesc2D->begin()+dd2->begin()[*cellId2], realIdsInDesc2D->begin()+dd2->begin()[*cellId2+1]);
+
+ std::vector<mcIdType> 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<DataArrayIdType> ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<mcIdType>::max()); ret3->rearrange(2);
MCAuto<DataArrayIdType> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
// deal with cells in mesh2D that are not cut but only some of their edges are
MCAuto<DataArrayIdType> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
+
MCAuto<DataArrayIdType> 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())
{
}
//
MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
- MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
- MCAuto<DataArrayIdType> elts,eltsIndex;
- mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
+ MCAuto<DataArrayDouble> baryRet1(centerOfMassRet1->selectByTupleId(idsInRet1NotColinear->begin(), idsInRet1NotColinear->end()));
+ DataArrayIdType *out(0),*outi(0);
+ DataArrayIdType::ExtractFromIndexedArrays(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end(),cells,cellsIndex,out,outi);
+ MCAuto<DataArrayIdType> elts(out),eltsIndex(outi);
+
MCAuto<DataArrayIdType> eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1);
if (eltsIndex->getNumberOfTuples() > 1)
eltsIndex2 = eltsIndex->deltaShiftIndex();
* 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)
MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
const double * coo(_coords->begin());
MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+ const double carMeshSz = getCaracteristicDimension();
{
/*************************
// Keep only candidates whose normal matches the normal of current face
for(vector<mcIdType>::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]),
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<mcIdType>::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)
INTERP_KERNEL::TranslationRotationMatrix rotation;
INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
- MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
+ MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is NOT called
MCAuto<DataArrayIdType> nodeMap(mPartCand->zipCoordsTraducer());
mcIdType nbElemsNotM1;
{
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<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
vector<std::size_t> compo; compo.push_back(1);
MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);