MEDCOUPLING_EXPORT DataArrayIdType *colinearize2D(double eps);
MEDCOUPLING_EXPORT DataArrayIdType *colinearizeKeepingConform2D(double eps);
MEDCOUPLING_EXPORT DataArrayIdType *conformize3D(double eps);
+ MEDCOUPLING_EXPORT DataArrayIdType *conformize3DIJK(double eps);
MEDCOUPLING_EXPORT mcIdType split2DCells(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI, const DataArrayIdType *midOpt=0, const DataArrayIdType *midOptI=0);
MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da);
MEDCOUPLING_EXPORT static MCAuto<MEDCouplingUMesh> Build1DMeshFromCoords(DataArrayDouble *da);
DataArrayIdType *internalColinearize2D(double eps, bool stayConform);
template<class MAPCLS>
void renumberNodesInConnT(const MAPCLS& newNodeNumbersO2N);
+ void conformize3DEdges(const double * coo, double eps, MCAuto<DataArrayIdType>& c, MCAuto<DataArrayIdType>& cI, MCAuto<DataArrayIdType>& ret);
+
public:
MEDCOUPLING_EXPORT static DataArrayIdType *ComputeRangesFromTypeDistribution(const std::vector<mcIdType>& code);
MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=25;
#include "InterpKernelGeo2DEdgeArcCircle.hxx"
#include "InterpKernelGeo2DQuadraticPolygon.hxx"
#include "MEDCouplingUMesh_internal.hxx"
+#include "TranslationRotationMatrix.hxx"
#include <sstream>
#include <fstream>
}
}
}
+
+
+void MEDCouplingUMesh::conformize3DEdges(const double * coo, double eps, MCAuto<DataArrayIdType>& c, MCAuto<DataArrayIdType>& cI, MCAuto<DataArrayIdType>& ret)
+{
+
+ static const int SPACEDIM=3;
+ // Now we have a face-conform mesh.
+
+ // Recompute descending
+ MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
+ // Rebuild desc connectivity with orientation this time!!
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
+ const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+ const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
+ const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
+ MCAuto<DataArrayIdType> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
+ MCAuto<DataArrayIdType> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
+ MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
+ MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
+ // std::cout << "writing!\n";
+ // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
+ // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
+ const mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
+ const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
+ MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
+ const double *bbox2(bboxArr->begin());
+ mcIdType nDesc2Cell=mDesc2->getNumberOfCells();
+ BBTree<SPACEDIM,mcIdType> myTree2(bbox2,0,0,nDesc2Cell,-eps);
+
+ // Edges - handle longest first
+ MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
+ DataArrayDouble * lens = lenF->getArray();
+
+ // Sort edges by decreasing length:
+ std::vector<std::pair<double,mcIdType> > S;
+ for(mcIdType i=0;i < lens->getNumberOfTuples();i++)
+ {
+ std::pair<double,mcIdType> p = std::make_pair(lens->getIJ(i, 0), i);
+ S.push_back(p);
+ }
+ sort(S.rbegin(),S.rend()); // reverse sort
+
+ std::vector<bool> hit(nDesc2Cell);
+ fill(hit.begin(), hit.end(), false);
+
+ for( std::vector<std::pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
+ {
+ mcIdType eIdx = (*it).second;
+ if (hit[eIdx])
+ continue;
+
+ std::vector<mcIdType> candidates, cands2;
+ myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
+ // Keep only candidates colinear with current edge
+ double vCurr[3];
+ 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];
+ for(std::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);
+ // 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)
+ cands2.push_back(*it2);
+ }
+ if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
+ continue;
+
+ // Now rotate edges to bring them on Ox
+ mcIdType startNode = cDesc2[cIDesc2[eIdx]+1];
+ mcIdType endNode = cDesc2[cIDesc2[eIdx]+2];
+ 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<DataArrayIdType> nodeMap(mPartCand->zipCoordsTraducer());
+ mcIdType nbElemsNotM1;
+ {
+ MCAuto<DataArrayIdType> tmp(nodeMap->findIdsNotEqual(-1));
+ nbElemsNotM1 = tmp->getNbOfElems();
+ }
+ MCAuto<DataArrayIdType> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
+ double * cooPartRef(mPartRef->_coords->getPointer());
+ double * cooPartCand(mPartCand->_coords->getPointer());
+ for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
+ rotation.transform_vector(cooPartRef+SPACEDIM*ii);
+ 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
+ MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
+ std::vector<std::size_t> compo; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
+ compo[0] = 2;
+ MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
+ MCAuto<DataArrayIdType> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayIdType> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayIdType> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
+ if (!idsGoodLine->getNumberOfTuples())
+ continue;
+
+ // Now the ordering along the Ox axis:
+ std::vector<mcIdType> insidePoints, hitSegs;
+ bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
+ mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
+ idsGoodLine->begin(), idsGoodLine->end(),
+ /*out*/insidePoints, hitSegs);
+ // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
+ for (std::vector<mcIdType>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
+ hit[cands2[*its]] = true;
+
+ if (!isSplit) // current segment remains in one piece
+ continue;
+
+ // Get original node IDs in global coords array
+ for (std::vector<mcIdType>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
+ *iit = nodeMapInv->begin()[*iit];
+
+ std::vector<mcIdType> polyIndices, packsIds, facePack;
+ // For each face implying this edge
+ for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
+ {
+ mcIdType faceIdx = revDescP2[ii];
+ // each cell where this face is involved connectivity will be modified:
+ ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
+
+ // Current face connectivity
+ const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
+ const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
+
+ std::vector<mcIdType> modifiedFace;
+ ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
+ modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
+ connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
+ }
+ }
+
+ // Rebuild 3D connectivity from descending:
+ MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
+ MCAuto<DataArrayIdType> superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
+ MCAuto<DataArrayIdType> idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0);
+ MCAuto<DataArrayIdType> vals(DataArrayIdType::New()); vals->alloc(0);
+ newConn->set3(superIdx, idx, vals);
+ mcIdType nbCells=getNumberOfCells();
+ for(mcIdType ii = 0; ii < nbCells; ii++)
+ for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++)
+ {
+ mcIdType sz, faceIdx = abs(descP[jj])-1;
+ bool orient = descP[jj]>0;
+ const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
+ if (orient)
+ newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
+ else
+ {
+ std::vector<mcIdType> rev(sz-1);
+ for (mcIdType kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
+ newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
+ }
+ }
+ // And finally:
+ newConn->convertToPolyhedronConn(c, cI);
+
+}
+
+
+
+
// Build circular permutation to shift consecutive edges together
renumb->iota(i+1);
renumb->applyModulus(nbCellsInSplitMesh1D);
+
splitMesh1D->renumberCells(renumbP, false);
cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
for(mcIdType mm=0;mm<nbCellsInSplitMesh1D;mm++)
pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
+
return pool.getZeMesh().retn();
}
MCAuto<DataArrayIdType> cIAuto; cIAuto.takeRef(_nodal_connec_index);
connSla->convertToPolyhedronConn(cAuto, cIAuto);
+
+ /************************
+ * STEP 2 -- edges
+ ************************/
+ conformize3DEdges(coo, eps, cAuto, cIAuto, ret);
+
+
+ ret = ret->buildUniqueNotSorted();
+ return ret.retn();
+}
+
+
+/*!
+ * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
+ * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
+ * This method performs a conformization of \b this for meshes aligned on Oxyz.
+ *
+ * Only polyhedron cells are supported. You can call convertAllToPoly()
+ *
+ * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
+ * 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.
+ * \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 mesh is not aligned on Oxyz.
+ * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
+ */
+DataArrayIdType *MEDCouplingUMesh::conformize3DIJK(double eps)
+{
+ using namespace std;
+
+ static const int SPACEDIM=3;
+ checkConsistencyLight();
+ if(getSpaceDimension()!=3 || getMeshDimension()!=3)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK : This method only works for meshes with spaceDim=3 and meshDim=3!");
+ if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK : This method only works for polyhedrons! Call convertAllToPoly first.");
+
+ MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
+ const double * coo(_coords->begin());
+ MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+
{
- /************************
- * STEP 2 -- edges
- ************************/
- // Now we have a face-conform mesh.
-
- // Recompute descending
- MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
- // Rebuild desc connectivity with orientation this time!!
- MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
+ /*************************
+ * STEP 1 -- faces
+ *************************/
+ MCAuto<DataArrayIdType> descDNU(DataArrayIdType::New()),descIDNU(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
- const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
- MCAuto<DataArrayIdType> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
- MCAuto<DataArrayIdType> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
- MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
- MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
- MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
-// std::cout << "writing!\n";
-// mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
-// mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
- const mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
- const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
- MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
- const double *bbox2(bboxArr->begin());
- mcIdType nDesc2Cell=mDesc2->getNumberOfCells();
- BBTree<SPACEDIM,mcIdType> myTree2(bbox2,0,0,nDesc2Cell,-eps);
-
- // Edges - handle longest first
- MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
- DataArrayDouble * lens = lenF->getArray();
-
- // Sort edges by decreasing length:
- vector<pair<double,mcIdType> > S;
- for(mcIdType i=0;i < lens->getNumberOfTuples();i++)
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
+
+ // Build BBTree
+ MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
+ const double *bbox(bboxArr->begin()); getCoords()->begin();
+ mcIdType nDescCell=mDesc->getNumberOfCells();
+ BBTree<SPACEDIM,mcIdType> myTree(bbox,0,0,nDescCell,-eps);
+ // Surfaces - handle biggest first
+ MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
+ DataArrayDouble * surfs = surfF->getArray();
+ // Normal field
+ MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
+ DataArrayDouble * normals = normalsF->getArray();
+ const double * normalsP = normals->getConstPointer();
+
+ // Sort faces by decreasing surface:
+ vector< pair<double,mcIdType> > S;
+ for(mcIdType i=0;i < surfs->getNumberOfTuples();i++)
{
- pair<double,mcIdType> p = make_pair(lens->getIJ(i, 0), i);
+ pair<double,mcIdType> p = make_pair(surfs->begin()[i], i);
S.push_back(p);
}
sort(S.rbegin(),S.rend()); // reverse sort
-
- vector<bool> hit(nDesc2Cell);
+ vector<bool> hit(nDescCell);
fill(hit.begin(), hit.end(), false);
+ vector<mcIdType> hitPoly; // the final result: which 3D cells have been modified.
for( vector<pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
{
- mcIdType eIdx = (*it).second;
- if (hit[eIdx])
- continue;
+ mcIdType faceIdx = (*it).second;
+ if (hit[faceIdx]) continue;
vector<mcIdType> candidates, cands2;
- myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
- // Keep only candidates colinear with current edge
- double vCurr[3];
- 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];
+ myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
+ // Keep only candidates whose normal matches the normal of current face
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);
- // 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)
+ bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
+ if (*it2 != faceIdx && col)
cands2.push_back(*it2);
}
- if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
+ if (!cands2.size())
continue;
- // Now rotate edges to bring them on Ox
- mcIdType startNode = cDesc2[cIDesc2[eIdx]+1];
- mcIdType endNode = cDesc2[cIDesc2[eIdx]+2];
+ MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
+ MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
+
+ // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
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<DataArrayIdType> nodeMap(mPartCand->zipCoordsTraducer());
- mcIdType nbElemsNotM1;
- {
- MCAuto<DataArrayIdType> tmp(nodeMap->findIdsNotEqual(-1));
- nbElemsNotM1 = tmp->getNbOfElems();
- }
- MCAuto<DataArrayIdType> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
+ INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
+ coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
+ coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
+
double * cooPartRef(mPartRef->_coords->getPointer());
double * cooPartCand(mPartCand->_coords->getPointer());
for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
rotation.transform_vector(cooPartCand+SPACEDIM*ii);
- // Eliminate all edges for which y or z is not null
+ // Localize faces in 2D thanks to barycenters
MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
- vector<std::size_t> compo; compo.push_back(1);
- MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
- compo[0] = 2;
+ vector<std::size_t> compo; compo.push_back(2);
MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
- MCAuto<DataArrayIdType> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
- MCAuto<DataArrayIdType> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
- MCAuto<DataArrayIdType> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
- if (!idsGoodLine->getNumberOfTuples())
+ MCAuto<DataArrayIdType> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
+ if (!idsGoodPlane->getNumberOfTuples())
continue;
- // Now the ordering along the Ox axis:
- std::vector<mcIdType> insidePoints, hitSegs;
- bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
- mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
- idsGoodLine->begin(), idsGoodLine->end(),
- /*out*/insidePoints, hitSegs);
- // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
- for (vector<mcIdType>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
- hit[cands2[*its]] = true;
-
- if (!isSplit) // current segment remains in one piece
- continue;
+ baryPart = baryPart->selectByTupleId(*idsGoodPlane);
- // Get original node IDs in global coords array
- for (std::vector<mcIdType>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
- *iit = nodeMapInv->begin()[*iit];
+ compo[0] = 0; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
+ mPartRef->changeSpaceDimension(2,0.0);
+ MCAuto<DataArrayIdType> cc(DataArrayIdType::New()), ccI(DataArrayIdType::New());
+ mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
- vector<mcIdType> polyIndices, packsIds, facePack;
- // For each face implying this edge
- for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
- {
- mcIdType faceIdx = revDescP2[ii];
- // each cell where this face is involved connectivity will be modified:
- ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
-
- // Current face connectivity
- const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
- const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
-
- vector<mcIdType> modifiedFace;
- ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
- modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
- connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
- }
- }
+ if (!cc->getNumberOfTuples())
+ continue;
+ MCAuto<DataArrayIdType> dsi = ccI->deltaShiftIndex();
- // Rebuild 3D connectivity from descending:
- MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
- MCAuto<DataArrayIdType> superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
- MCAuto<DataArrayIdType> idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0);
- MCAuto<DataArrayIdType> vals(DataArrayIdType::New()); vals->alloc(0);
- newConn->set3(superIdx, idx, vals);
- mcIdType nbCells=getNumberOfCells();
- for(mcIdType ii = 0; ii < nbCells; ii++)
- for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++)
{
- mcIdType sz, faceIdx = abs(descP[jj])-1;
- bool orient = descP[jj]>0;
- const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
- if (orient)
- newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
- else
+ MCAuto<DataArrayIdType> tmp = dsi->findIdsInRange(0, 2);
+ if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
{
- vector<mcIdType> rev(sz-1);
- for (mcIdType kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
- newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
+ ostringstream oss;
+ oss << "MEDCouplingUMesh::conformize3DIJK: Non expected non-conformity. Only simple non-conformities are handled. Face #" << faceIdx << " violates this condition!";
+ throw INTERP_KERNEL::Exception(oss.str());
}
}
- // And finally:
- newConn->convertToPolyhedronConn(cAuto, cIAuto);
- } // end step2
+
+ MCAuto<DataArrayIdType> ids = dsi->findIdsEqual(1);
+ // Boundary face:
+ if (!ids->getNumberOfTuples())
+ continue;
+
+ const mcIdType * idsGoodPlaneP(idsGoodPlane->begin());
+ vector<mcIdType> goFaces; //Faces that intersect current face
+ for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++)
+ {
+ mcIdType faceIdx2 = cands2[idsGoodPlaneP[*ii]];
+ goFaces.push_back(faceIdx2);
+ hit[faceIdx2] = true;
+ }
+
+ MCAuto<MEDCouplingUMesh> meshGoodPlane(mDesc->buildPartOfMySelf(&goFaces[0], &goFaces[0]+goFaces.size(), false)); // false=zipCoords is called
+
+ double faceIdxNormal = normalsP[faceIdx];
+ int direction = -1;
+ for(int d=0; d<3; d++)
+ {
+ double diff = fabs(normalsP[faceIdx*3+d]) - 1;
+ if(fabs(diff) < eps)
+ direction = d;
+ }
+ if(direction == -1)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK: Not a mesh aligned on Oxyz !");
+
+ // Projecting the plane in a 2D space
+ const DataArrayDouble *cooGoodPlane(meshGoodPlane->getCoords());
+ double extraDim = cooGoodPlane->begin()[direction];
+ map<int, pair<int,int> > projection = {{0, make_pair(1,2)}, {1, make_pair(0,2)}, {2, make_pair(0,1)}};
+ int dir1 = projection.at(direction).first, dir2 = projection.at(direction).second;
+ vector<std::size_t> compo2D; compo2D.push_back(dir1); compo2D.push_back(dir2);
+ MCAuto<DataArrayDouble> coo2D = cooGoodPlane->keepSelectedComponents(compo2D);
+ meshGoodPlane->setCoords(coo2D);
+
+ // Intersection between the plane and the current face
+ const double vec[3] = {0.,0.,1.};
+ mPartRef->changeSpaceDimension(3,0.0); meshGoodPlane->changeSpaceDimension(3,0.0);
+ mPartRef->orientCorrectly2DCells(vec, false); meshGoodPlane->orientCorrectly2DCells(vec, false);
+ mPartRef->changeSpaceDimension(2,0.0); meshGoodPlane->changeSpaceDimension(2,0.0);
+ DataArrayIdType *cellIdInPlane(0),*cellIdInCurrentFace(0);
+ MCAuto<MEDCouplingUMesh> mInter=MEDCouplingUMesh::Intersect2DMeshes(meshGoodPlane,mPartRef,eps,cellIdInPlane,cellIdInCurrentFace);
+ const DataArrayDouble *cooInter(mInter->getCoords());
+ const double * coo_InterD(cooInter->begin());
+
+ goFaces.push_back(faceIdx);
+ mcIdType index = 0;
+ for( vector<mcIdType>::const_iterator face=goFaces.begin(); face!=goFaces.end(); ++face, ++index )
+ {
+ // For all polyhedrons using this face, replace connectivity:
+ vector<mcIdType> polyIndices, packsIds, facePack;
+ for (mcIdType ii=revDescIP[*face]; ii < revDescIP[*face+1]; ii++)
+ polyIndices.push_back(revDescP[ii]);
+ ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
+
+ // Current face connectivity
+ const mcIdType * sIdxConn = cDesc + cIDesc[*face] + 1;
+ const mcIdType * sIdxConnE = cDesc + cIDesc[*face+1];
+ connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
+
+ // Deletion of old faces
+ mcIdType jj=0;
+ for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
+ {
+ if (packsIds[jj] == -1)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3DIJK: Could not find face in connectivity! Internal error.");
+ else
+ connSla->deletePack(*it2, packsIds[jj]);
+ }
+
+ // Insertion of new faces coming from the intersection
+ MCAuto<DataArrayIdType> newFaces= *face==goFaces.back() ? cellIdInCurrentFace->findIdsEqual(0) : cellIdInPlane->findIdsEqual(index);
+ if(newFaces->getNumberOfTuples())
+ {
+ for (const mcIdType *conformFace = newFaces->begin(); conformFace != newFaces->end(); conformFace++)
+ {
+ MCAuto<MEDCouplingSkyLineArray> connSlaDescIntersectMesh(MEDCouplingSkyLineArray::New(mInter->getNodalConnectivityIndex(), mInter->getNodalConnectivity()));
+ mcIdType facePack2Sz;
+ const mcIdType * facePack2 = connSlaDescIntersectMesh->getSimplePackSafePtr(*conformFace, facePack2Sz); // contains the type!
+
+ vector<mcIdType> origNodes;
+ for(mcIdType nodeI=1; nodeI <facePack2Sz; nodeI++ ) //without type
+ {
+ DataArrayDouble *cooIn3DMesh(DataArrayDouble::New());
+ double cooIn3DMesh_[3];
+ cooIn3DMesh_[direction] = extraDim; cooIn3DMesh_[dir1] = coo_InterD[2*facePack2[nodeI]]; cooIn3DMesh_[dir2] = coo_InterD[2*facePack2[nodeI]+1];
+ cooIn3DMesh->useArray(cooIn3DMesh_,false,DeallocType::C_DEALLOC,1,3);
+ MCAuto<DataArrayIdType> nodeIndexIn3DMesh = _coords->findClosestTupleId(cooIn3DMesh);
+ origNodes.push_back(nodeIndexIn3DMesh->begin()[0]);
+ }
+
+ for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
+ connSla->pushBackPack(*it2, origNodes.data(), origNodes.data()+origNodes.size()); // without the type
+ }
+ }
+ else
+ {
+ mcIdType facePack2Sz;
+ const mcIdType * facePack2 = connSlaDesc->getSimplePackSafePtr(*face, facePack2Sz); // contains the type!
+ // Insert it in all hit polyhedrons:
+ for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
+ connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
+ }
+ }
+ cellIdInCurrentFace->decrRef();
+ cellIdInPlane->decrRef();
+ }
+ } // end step1
+
+ // Set back modified connectivity
+ MCAuto<DataArrayIdType> cAuto; cAuto.takeRef(_nodal_connec);
+ MCAuto<DataArrayIdType> cIAuto; cIAuto.takeRef(_nodal_connec_index);
+ connSla->convertToPolyhedronConn(cAuto, cIAuto);
+
+
+ /************************
+ * STEP 2 -- edges
+ ************************/
+ conformize3DEdges(coo, eps, cAuto, cIAuto, ret);
ret = ret->buildUniqueNotSorted();
return ret.retn();
}
+
+
//tools
DataArrayIdType *conformize2D(double eps);
DataArrayIdType *conformize3D(double eps);
+ DataArrayIdType *conformize3DIJK(double eps);
DataArrayIdType *colinearize2D(double eps);
DataArrayIdType *colinearizeKeepingConform2D(double eps);
void shiftNodeNumbersInConn(int delta);
self.assertEqual(set([18]), set(ret.getValues()))
pass
+ def testSwig2Conformize3DIJK(self):
+ eps = 1.0e-8
+ mesh = MEDCouplingUMesh('merge', 3)
+ coo = DataArrayDouble([(0.0, 0.0, -0.5), (1.0, 0.0, -0.5), (0.0, 1.0, -0.5), (1.0, 1.0, -0.5), (0.0, 0.0, 0.5), (1.0, 0.0, 0.5), (0.0, 1.0, 0.5), (1.0, 1.0, 0.5), (0.0, 0.0, 1.5), (1.0, 0.0, 1.5), (0.0, 1.0, 1.5), (1.0, 1.0, 1.5), (-1.0, 0.0, 0.0), (0.0, 0.0, 0.0), (-1.0, 1.0, 0.0), (0.0, 1.0, 0.0), (-1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-1.0, 1.0, 1.0), (0.0, 1.0, 1.0)])
+ mesh.setCoords(coo)
+ c = DataArrayInt([31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3, -1, 3, 7, 5, 1, 31, 5, 4, 6, 7, -1, 9, 11, 10, 8, -1, 5, 9, 8, 4, -1, 4, 8, 10, 6, -1, 6, 10, 11, 7, -1, 7, 11, 9, 5, 31, 13, 12, 14, 15, -1, 17, 19, 18, 16, -1, 13, 17, 16, 12, -1, 12, 16, 18, 14, -1, 14, 18, 19, 15, -1, 15, 19, 17, 13])
+ cI = DataArrayInt([0, 30, 60, 90])
+ mesh.setConnectivity(c, cI)
+
+ ret = mesh.conformize3DIJK(eps)
+
+ mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity()
+ mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity()
+ c0, cI0 = mesh.getNodalConnectivity().getValues(), mesh.getNodalConnectivityIndex().getValues()
+ c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues()
+ c2, cI2 = mretDesc2.getNodalConnectivity().getValues(), mretDesc2.getNodalConnectivityIndex().getValues()
+
+ cRef0 = [31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 13, 0, -1, 2, 15, 6, 7, 3, -1, 3, 7, 5, 1, -1, 15, 6, 4, 13, -1, 13, 0, 2, 15, 31, 4, 6, 7, 5, -1, 9, 11, 10, 8, -1, 5, 9, 8, 17, 4, -1, 6, 19, 10, 11, 7, -1, 7, 11, 9, 5, -1, 17, 4, 6, 19, -1, 19, 10, 8, 17, 31, 13, 12, 14, 15, -1, 17, 19, 18, 16, -1, 13, 4, 17, 16, 12, -1, 12, 16, 18, 14, -1, 14, 18, 19, 6, 15, -1, 15, 6, 4, 13, -1, 17, 4, 6, 19]
+ cIRef0 = [0, 37, 74, 111]
+ cRef = [5, 1, 0, 2, 3, 5, 5, 7, 6, 4, 5, 1, 5, 4, 13, 0, 5, 2, 15, 6, 7, 3, 5, 3, 7, 5, 1, 5, 15, 6, 4, 13, 5, 13, 0, 2, 15, 5, 9, 11, 10, 8, 5, 5, 9, 8, 17, 4, 5, 6, 19, 10, 11, 7, 5, 7, 11, 9, 5, 5, 17, 4, 6, 19, 5, 19, 10, 8, 17, 5, 13, 12, 14, 15, 5, 17, 19, 18, 16, 5, 13, 4, 17, 16, 12, 5, 12, 16, 18, 14, 5, 14, 18, 19, 6, 15]
+ cIRef = [0, 5, 10, 16, 22, 27, 32, 37, 42, 48, 54, 59, 64, 69, 74, 79, 85, 90, 96]
+ cRef2 = [1, 1, 0, 1, 0, 2, 1, 2, 3, 1, 3, 1, 1, 5, 7, 1, 7, 6, 1, 6, 4, 1, 4, 5, 1, 1, 5, 1, 4, 13, 1, 13, 0, 1, 2, 15, 1, 15, 6, 1, 7, 3, 1, 13, 15, 1, 9, 11, 1, 11, 10, 1, 10, 8, 1, 8, 9, 1, 5, 9, 1, 8, 17, 1, 17, 4, 1, 6, 19, 1, 19, 10, 1, 11, 7, 1, 19, 17, 1, 13, 12, 1, 12, 14, 1, 14, 15, 1, 19, 18, 1, 18, 16, 1, 16, 17, 1, 16, 12, 1, 18, 14]
+ cIRef2 = [0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78, 81, 84, 87, 90, 93, 96, 99, 102]
+ self.assertEqual(18, mretDesc.getNumberOfCells())
+ self.assertEqual(34, mretDesc2.getNumberOfCells())
+ self.assertEqual(cRef0, c0)
+ self.assertEqual(cIRef0, cI0)
+ self.assertEqual(cRef, c)
+ self.assertEqual(cIRef, cI)
+ self.assertEqual(cRef2, c2)
+ self.assertEqual(cIRef2, cI2)
+ self.assertEqual(set([0,1,2]), set(ret.getValues()))
+ pass
+
+
if __name__ == '__main__':
unittest.main()
