#include "MEDCouplingUMesh.hxx"
#include "MEDCoupling1GTUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "CellModel.hxx"
#include "VolSurfUser.txx"
#include "InterpolationUtils.hxx"
#include "InterpKernelGeo2DEdgeLin.hxx"
#include "InterpKernelGeo2DEdgeArcCircle.hxx"
#include "InterpKernelGeo2DQuadraticPolygon.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "VectorUtils.hxx"
+#include "MEDCouplingSkyLineArray.hxx"
#include <sstream>
#include <fstream>
#include <cstring>
#include <limits>
#include <list>
+#include <assert.h>
using namespace MEDCoupling;
}
/**
- * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI).
+ * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
*/
void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
// This initializes posBaseElt.
if(nbOfTurn==0)
{
- for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell wih one single edge
+ for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
{
cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
}
// Now move forward:
const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
- for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell wih one single edge
+ for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell with one single edge
{
cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #j's connectivity
INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
break;
}
//push [posBaseElt,posEndElt) in newConnOfCell using e
- // The if clauses below are (volontary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its begining!
+ // The if clauses below are (voluntary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its beginning!
if(nbOfTurn==0)
- // at the begining of the connectivity (insert type)
+ // at the beginning of the connectivity (insert type)
EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
else if((nbOfHit+nbOfTurn) != (nbs-1))
// in the middle
VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
std::size_t size() const { return _pool.size(); }
int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
- void setMeshAt(int pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs);
+ void setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs);
const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
return zeMesh->getCellContainingPoint(barys->begin(),eps);
}
-void VectorOfCellInfo::setMeshAt(int pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
+void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
{
get(pos);//to check pos
bool isFast(pos==0 && _pool.size()==1);
_edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
//
std::vector<CellInfo> pool(_pool.size()-1+sz);
+ for(
std::size_t i=0;i<pos;i++)
pool[i]=_pool[i];
for(std::size_t j=0;j<sz;j++)
pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2, std::vector<double>& addCoo, std::map<int,int>& mergedNodes)
{
static const int SPACEDIM=2;
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
// Build BB tree of all edges in the tool mesh (second mesh)
- MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree());
+ MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
intersectEdge1.resize(nDescCell1);
const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
int offset2(offset1+m2->getNumberOfNodes());
int offset3(offset2+((int)addCoords.size())/2);
- MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
+ MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
// Here a BBTree on 2D-cells, not on segments:
BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
for(int i=0;i<nbOfCells;i++)
{
- int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset),start(-1),end(-1);
+ int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
for(int j=0;j<nbOfFaces;j++)
{
const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
}
/*!
- * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
+ * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additional nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
*
* \sa MEDCouplingUMesh::split2DCells
*/
/*!
- * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
+ * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
*
* \return int - the number of new nodes created.
* \sa MEDCouplingUMesh::split2DCells
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
m1->checkFullyDefined();
m2->checkFullyDefined();
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
/*!
* Partitions the first given 2D mesh using the second given 1D mesh as a tool.
* Thus the final result contains the aggregation of nodes of \a mesh2D, then nodes of \a mesh1D, then new nodes that are the result of the intersection
- * and finaly, in case of quadratic polygon the centers of edges new nodes.
+ * 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.
*
* \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
// Step 1: compute all edge intersections (new nodes)
std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
//
// Build desc connectivity
DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
MCAuto<DataArrayInt> elts,eltsIndex;
mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
- MCAuto<DataArrayInt> eltsIndex2(eltsIndex->deltaShiftIndex());
+ MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
+ if (eltsIndex->getNumberOfTuples() > 1)
+ eltsIndex2 = eltsIndex->deltaShiftIndex();
MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
MCAuto<MEDCouplingUMesh> splitOfOneCell(BuildMesh2DCutFrom(eps,*it,m1Desc,partOfMesh1CuttingCur2DCell,dd1->begin()+dd2->getIJ(*it,0),dd1->begin()+dd2->getIJ((*it)+1,0),intersectEdge1,ret2->getNumberOfTuples(),partOfRet3));
ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
outMesh2DSplit.push_back(splitOfOneCell);
- for(int i=0;i<splitOfOneCell->getNumberOfCells();i++)
+ for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
ret2->pushBackSilent(*it);
}
//
MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
// To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
ret3->rearrange(1);
- MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStricltyNegative());
+ MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
{
int old2DCellId(-ret3->getIJ(*it,0)-1);
MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
- MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
+ MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
std::vector<double> addCoo;
BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
for(int i=0;i<nDescCell;i++)
{
std::vector<int> candidates;
checkConsistencyLight();
if(getSpaceDimension()!=2 || getMeshDimension()!=2)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
return ret.retn();
}
+///@cond INTERNAL
+/**
+ * c, cI describe a wire mesh in 3D space, in local numbering
+ * startNode, endNode in global numbering
+ *\return true if the segment is indeed split
+ */
+bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
+ const int * c, const int * cI, const int *idsBg, const int *endBg,
+ std::vector<int> & pointIds, std::vector<int> & hitSegs)
+{
+ using namespace std;
+
+ const int SPACEDIM=3;
+ typedef pair<double, int> PairDI;
+ set< PairDI > x;
+ for (const int * it = idsBg; it != endBg; ++it)
+ {
+ assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
+ int start = c[cI[*it]+1], end = c[cI[*it]+2];
+ x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
+ x.insert(make_pair(coo[end*SPACEDIM], end));
+ }
+
+ vector<PairDI> xx(x.begin(), x.end());
+ sort(xx.begin(),xx.end());
+ pointIds.reserve(xx.size());
+
+ // Keep what is inside [startNode, endNode]:
+ int go = 0;
+ for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
+ {
+ const int idx = (*it).second;
+ if (!go)
+ {
+ if (idx == startNode) go = 1;
+ if (idx == endNode) go = 2;
+ if (go) pointIds.push_back(idx);
+ continue;
+ }
+ pointIds.push_back(idx);
+ if (idx == endNode || idx == startNode)
+ break;
+ }
+
+// vector<int> pointIds2(pointIds.size()+2);
+// copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
+// pointIds2[0] = startNode;
+// pointIds2[pointIds2.size()-1] = endNode;
+
+ if (go == 2)
+ reverse(pointIds.begin(), pointIds.end());
+
+ // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
+ for (const int * it = idsBg; it != endBg; ++it)
+ {
+ int start = c[cI[*it]+1], end = c[cI[*it]+2];
+ vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
+ if (itStart == pointIds.end()) continue;
+ vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
+ if (itEnd == pointIds.end()) continue;
+ if (abs(distance(itEnd, itStart)) != 1) continue;
+ hitSegs.push_back(*it); // segment is undivided.
+ }
+
+ return (pointIds.size() > 2); // something else apart start and end node
+}
+
+void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
+ const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
+{
+ using namespace std;
+ int dst = distance(sIdxConn, sIdxConnE);
+ modifiedFace.reserve(dst + insidePoints.size()-2);
+ modifiedFace.resize(dst);
+ copy(sIdxConn, sIdxConnE, modifiedFace.data());
+
+ vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
+ vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
+ if (startPos == shortEnd)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
+ vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
+ if (endPos == shortEnd)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
+ int d = distance(startPos, endPos);
+ if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
+ modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
+ else
+ modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
+}
+
+///@endcond
+
+
+/*!
+ * \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.
+ *
+ * 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 DataArrayInt * - 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.
+ * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
+ */
+DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
+{
+ using namespace std;
+
+ static const int SPACEDIM=3;
+ checkConsistencyLight();
+ if(getSpaceDimension()!=3 || getMeshDimension()!=3)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : 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::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
+
+ MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
+ const double * coo(_coords->begin());
+ MCAuto<DataArrayInt> ret(DataArrayInt::New());
+
+ {
+ /*************************
+ * STEP 1 -- faces
+ *************************/
+ MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
+ const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+ const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
+
+ // Build BBTree
+ MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
+ const double *bbox(bboxArr->begin()); getCoords()->begin();
+ int nDescCell(mDesc->getNumberOfCells());
+ BBTree<SPACEDIM,int> 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,int> > S;
+ for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
+ {
+ pair<double,int> p = make_pair(surfs->begin()[i], i);
+ S.push_back(p);
+ }
+ sort(S.rbegin(),S.rend()); // reverse sort
+ vector<bool> hit(nDescCell);
+ fill(hit.begin(), hit.end(), false);
+ vector<int> hitPoly; // the final result: which 3D cells have been modified.
+
+ for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
+ {
+ int faceIdx = (*it).second;
+ if (hit[faceIdx]) continue;
+
+ vector<int> candidates, cands2;
+ myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
+ // Keep only candidates whose normal matches the normal of current face
+ for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
+ {
+ bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
+ 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
+ INTERP_KERNEL::TranslationRotationMatrix rotation;
+ INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
+ coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
+ coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
+
+ 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
+ double * cooPartRef(mPartRef->_coords->getPointer());
+ double * cooPartCand(mPartCand->_coords->getPointer());
+ for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
+ rotation.transform_vector(cooPartRef+SPACEDIM*ii);
+ for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
+ rotation.transform_vector(cooPartCand+SPACEDIM*ii);
+
+ // Localize faces in 2D thanks to barycenters
+ MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
+ vector<int> compo; compo.push_back(2);
+ MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
+ MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
+ if (!idsGoodPlane->getNumberOfTuples())
+ continue;
+
+ baryPart = baryPart->selectByTupleId(*idsGoodPlane);
+
+ compo[0] = 0; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
+ mPartRef->changeSpaceDimension(2,0.0);
+ MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
+ mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
+
+ if (!cc->getNumberOfTuples())
+ continue;
+ MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
+
+ {
+ MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
+ if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
+ {
+ ostringstream oss;
+ oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+ }
+
+ MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
+ // Boundary face:
+ if (!ids->getNumberOfTuples())
+ continue;
+
+ double checkSurf=0.0;
+ const int * idsGoodPlaneP(idsGoodPlane->begin());
+ for (const int * ii = ids->begin(); ii != ids->end(); ii++)
+ {
+ int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
+ hit[faceIdx2] = true;
+ checkSurf += surfs->begin()[faceIdx2];
+ }
+ if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
+ {
+ ostringstream oss;
+ oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+
+ // For all polyhedrons using this face, replace connectivity:
+ vector<int> polyIndices, packsIds, facePack;
+ for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
+ polyIndices.push_back(revDescP[ii]);
+ ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
+
+ // Current face connectivity
+ const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
+ const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
+ connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
+ // Deletion of old faces
+ int jj=0;
+ for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
+ {
+ if (packsIds[jj] == -1)
+ // The below should never happen - if a face is used several times, with a different layout of the nodes
+ // it means that is is already conform, so it is *not* hit by the algorithm. The algorithm only hits
+ // faces which are actually used only once, by a single cell. This is different for edges below.
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
+ else
+ connSla->deletePack(*it2, packsIds[jj]);
+ }
+ // Insertion of new faces:
+ for (const int * ii = ids->begin(); ii != ids->end(); ii++)
+ {
+ // Build pack from the face to insert:
+ int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
+ int facePack2Sz;
+ const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
+ // Insert it in all hit polyhedrons:
+ for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
+ connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
+ }
+ }
+ } // end step1
+
+ // Set back modified connectivity
+ MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
+ MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
+ connSla->convertToPolyhedronConn(cAuto, cIAuto);
+
+ {
+ /************************
+ * STEP 2 -- edges
+ ************************/
+ // Now we have a face-conform mesh.
+
+ // Recompute descending
+ MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
+ // Rebuild desc connectivity with orientation this time!!
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
+ const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+ const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
+ const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
+ MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
+ MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
+ MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::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 int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
+ const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
+ MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
+ const double *bbox2(bboxArr->begin());
+ int nDesc2Cell=mDesc2->getNumberOfCells();
+ BBTree<SPACEDIM,int> 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,int> > S;
+ for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
+ {
+ pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
+ S.push_back(p);
+ }
+ sort(S.rbegin(),S.rend()); // reverse sort
+
+ vector<bool> hit(nDesc2Cell);
+ fill(hit.begin(), hit.end(), false);
+
+ for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
+ {
+ int eIdx = (*it).second;
+ if (hit[eIdx])
+ continue;
+
+ vector<int> candidates, cands2;
+ myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
+ // Keep only candidates colinear with current edge
+ double vCurr[3];
+ unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
+ for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
+ vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
+ for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
+ {
+ double vOther[3];
+ unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
+ for (int 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
+ int startNode = cDesc2[cIDesc2[eIdx]+1];
+ int 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<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
+ int nbElemsNotM1;
+ {
+ MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
+ nbElemsNotM1 = tmp->getNbOfElems();
+ }
+ MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
+ double * cooPartRef(mPartRef->_coords->getPointer());
+ double * cooPartCand(mPartCand->_coords->getPointer());
+ for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
+ rotation.transform_vector(cooPartRef+SPACEDIM*ii);
+ for (std::size_t 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();
+ vector<int> compo; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
+ compo[0] = 2;
+ MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
+ MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
+ if (!idsGoodLine->getNumberOfTuples())
+ continue;
+
+ // Now the ordering along the Ox axis:
+ std::vector<int> 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<int>::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<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
+ *iit = nodeMapInv->begin()[*iit];
+
+ vector<int> polyIndices, packsIds, facePack;
+ // For each face implying this edge
+ for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
+ {
+ int 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 int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
+ const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
+
+ vector<int> 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<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
+ MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
+ MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
+ newConn->set3(superIdx, idx, vals);
+ for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
+ for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
+ {
+ int sz, faceIdx = abs(descP[jj])-1;
+ bool orient = descP[jj]>0;
+ const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
+ if (orient)
+ newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
+ else
+ {
+ vector<int> rev(sz-1);
+ for (int 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(cAuto, cIAuto);
+ } // end step2
+
+ ret = ret->buildUniqueNotSorted();
+ return ret.retn();
+}
+
