-// 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
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();
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]);
* 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);
