-// Copyright (C) 2007-2021 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
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();
{
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 !");
* 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);
