1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingUMesh.hxx"
22 #include "MEDCoupling1GTUMesh.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
24 #include "CellModel.hxx"
25 #include "VolSurfUser.txx"
26 #include "InterpolationUtils.hxx"
27 #include "PointLocatorAlgos.txx"
29 #include "BBTreeDst.txx"
30 #include "DirectedBoundingBox.hxx"
31 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
32 #include "InterpKernelAutoPtr.hxx"
33 #include "InterpKernelGeo2DNode.hxx"
34 #include "InterpKernelGeo2DEdgeLin.hxx"
35 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
36 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
37 #include "TranslationRotationMatrix.hxx"
38 #include "VectorUtils.hxx"
39 #include "MEDCouplingSkyLineArray.hxx"
49 using namespace MEDCoupling;
53 int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
59 int ret(nodesCnter++);
61 e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
62 addCoo.insertAtTheEnd(newPt,newPt+2);
67 int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
73 int ret(nodesCnter++);
75 e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt);
76 addCoo.insertAtTheEnd(newPt,newPt+2);
82 void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
85 int trueStart(start>=0?start:nbOfEdges+start);
86 tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
87 newConnOfCell->insertAtTheEnd(tmp,tmp+3);
92 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
93 InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2);
94 middles.push_back(tmp3+offset);
97 middles.push_back(connBg[trueStart+nbOfEdges]);
101 void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
103 int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
104 newConnOfCell->pushBackSilent(tmpEnd);
109 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
110 InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
111 middles.push_back(tmp3+offset);
114 middles.push_back(connBg[start+nbOfEdges]);
118 void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
120 // only the quadratic point to deal with:
123 if(stp-start>1) // if we are covering more than one segment we need to create a new mid point
125 int tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]); // % to handle last seg.
126 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
127 InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
128 middles.push_back(tmp3+offset);
131 middles.push_back(connBg[start+nbOfEdges]);
135 void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
137 MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
138 std::map<MCAuto<INTERP_KERNEL::Node>,int>::const_iterator it(m.find(nTmp));
140 throw INTERP_KERNEL::Exception("Internal error in remapping !");
142 if(v==forbVal0 || v==forbVal1)
144 if(std::find(isect.begin(),isect.end(),v)==isect.end())
148 bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
153 bool presenceOfOn(false);
154 for(int i=0;i<sz;i++)
156 INTERP_KERNEL::ElementaryEdge *e(c[i]);
157 if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
159 IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
160 IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
166 namespace MEDCoupling
169 INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map< MCAuto<INTERP_KERNEL::Node>,int>& m)
171 INTERP_KERNEL::Edge *ret(0);
172 MCAuto<INTERP_KERNEL::Node> n0(new INTERP_KERNEL::Node(coords2D[2*bg[0]],coords2D[2*bg[0]+1])),n1(new INTERP_KERNEL::Node(coords2D[2*bg[1]],coords2D[2*bg[1]+1]));
173 m[n0]=bg[0]; m[n1]=bg[1];
176 case INTERP_KERNEL::NORM_SEG2:
178 ret=new INTERP_KERNEL::EdgeLin(n0,n1);
181 case INTERP_KERNEL::NORM_SEG3:
183 INTERP_KERNEL::Node *n2(new INTERP_KERNEL::Node(coords2D[2*bg[2]],coords2D[2*bg[2]+1])); m[n2]=bg[2];
184 INTERP_KERNEL::EdgeLin *e1(new INTERP_KERNEL::EdgeLin(n0,n2)),*e2(new INTERP_KERNEL::EdgeLin(n2,n1));
185 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
186 // is the SEG3 degenerated, and thus can be reduced to a SEG2?
187 bool colinearity(inters.areColinears());
188 delete e1; delete e2;
190 { ret=new INTERP_KERNEL::EdgeLin(n0,n1); }
192 { ret=new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1); }
196 throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge2 : Expecting a mesh with spaceDim==2 and meshDim==1 !");
201 INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
203 INTERP_KERNEL::Edge *ret=0;
206 case INTERP_KERNEL::NORM_SEG2:
208 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
211 case INTERP_KERNEL::NORM_SEG3:
213 INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
214 INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
215 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
216 // is the SEG3 degenerated, and thus can be reduced to a SEG2?
217 bool colinearity=inters.areColinears();
218 delete e1; delete e2;
220 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
222 ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
223 mapp2[bg[2]].second=false;
227 throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
233 * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from
234 * the global mesh 'mDesc'.
235 * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2.
236 * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
238 INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
239 std::map<INTERP_KERNEL::Node *,int>& mapp)
242 std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3.
243 const double *coo=mDesc->getCoords()->getConstPointer();
244 const int *c=mDesc->getNodalConnectivity()->getConstPointer();
245 const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
247 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
248 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
249 for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
251 INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
252 mapp2[*it2]=std::pair<INTERP_KERNEL::Node *,bool>(n,true);
254 INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
255 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
257 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
258 ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
260 for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
262 if((*it2).second.second)
263 mapp[(*it2).second.first]=(*it2).first;
264 ((*it2).second.first)->decrRef();
269 INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
273 int locId=nodeId-offset2;
274 return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
278 int locId=nodeId-offset1;
279 return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
281 return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
285 * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
287 void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
288 const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
289 /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
291 for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
293 int eltId1=abs(*desc1)-1;
294 for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
296 std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
297 if(it==mappRev.end())
299 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
311 * Returns true if a colinearization has been found in the given cell. If false is returned the content pushed in \a newConnOfCell is equal to [ \a connBg , \a connEnd ) .
312 * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
313 * \param forbiddenPoints the list of points that should not be removed in the process
315 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset,
316 const std::map<int, bool>& forbiddenPoints,
317 DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
319 std::size_t sz(std::distance(connBg,connEnd));
320 if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
321 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
323 INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
324 INTERP_KERNEL::AutoPtr<int> tmpConn2(new int[sz]);
325 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
326 unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
327 unsigned nbOfHit(0); // number of fusions operated
328 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
329 const unsigned int maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3; // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least
330 INTERP_KERNEL::NormalizedCellType typeOfSon;
331 std::vector<int> middles;
333 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
335 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
336 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
337 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
338 posEndElt = posBaseElt+1;
340 // Look backward first: are the final edges of the cells colinear with the first ones?
341 // This initializes posBaseElt.
344 for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
346 cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn2,typeOfSon);
347 // Identify common point:
348 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
349 auto itE(forbiddenPoints.end());
350 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
352 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
353 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
354 bool isColinear=eint->areColinears();
365 // Update last connectivity
366 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
370 const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
371 for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell with one single edge
373 cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn2,typeOfSon); // get edge #j's connectivity
374 // Identify common point:
375 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
376 auto itE(forbiddenPoints.end());
377 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
379 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
380 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
381 bool isColinear(eint->areColinears());
392 // Update last connectivity
393 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
395 //push [posBaseElt,posEndElt) in newConnOfCell using e
396 // 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!
398 // at the beginning of the connectivity (insert type)
399 EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
400 else if((nbOfHit+nbOfTurn) != (nbs-1))
402 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
403 if ((nbOfHit+nbOfTurn) == (nbs-1))
404 // at the end (only quad points to deal with)
405 EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
406 posBaseElt=posEndElt;
410 newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
416 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
418 if(candidates.empty())
420 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
422 const std::vector<int>& pool(intersectEdge1[*it]);
423 int tmp[2]; tmp[0]=start; tmp[1]=stop;
424 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
429 tmp[0]=stop; tmp[1]=start;
430 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
440 * This method performs the 2nd step of Partition of 2D mesh.
441 * This method has 4 inputs :
442 * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
443 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
444 * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
445 * The aim of this method is to sort the splitting nodes, if any, and to put them in 'intersectEdge' output parameter based on edges of mesh 'm2'
446 * Nodes end up lying consecutively on a cutted edge.
447 * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method.
448 * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
449 * \param m2 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method.
450 * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
451 * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
453 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
454 const std::vector<double>& addCoo,
455 const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
457 int offset1=m1->getNumberOfNodes();
458 int ncell2=m2->getNumberOfCells();
459 const int *c=m2->getNodalConnectivity()->begin();
460 const int *cI=m2->getNodalConnectivityIndex()->begin();
461 const double *coo=m2->getCoords()->begin();
462 const double *cooBis=m1->getCoords()->begin();
463 int offset2=offset1+m2->getNumberOfNodes();
464 intersectEdge.resize(ncell2);
465 for(int i=0;i<ncell2;i++,cI++)
467 const std::vector<int>& divs=subDiv[i];
468 int nnode=cI[1]-cI[0]-1;
469 std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
470 std::map<INTERP_KERNEL::Node *, int> mapp22;
471 for(int j=0;j<nnode;j++)
473 INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
474 int nnid=c[(*cI)+j+1];
475 mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
476 mapp22[nn]=nnid+offset1;
478 INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
479 for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
480 ((*it).second.first)->decrRef();
481 std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
482 std::map<INTERP_KERNEL::Node *,int> mapp3;
483 for(std::size_t j=0;j<divs.size();j++)
486 INTERP_KERNEL::Node *tmp=0;
488 tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
490 tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
492 tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
496 e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
497 for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
503 MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<int> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<int,int>& mergedNodes, const std::vector< std::vector<int> >& colinear2, const std::vector< std::vector<int> >& intersectEdge1,
504 MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
506 idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
507 idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
508 int nCells(mesh1D->getNumberOfCells());
509 if(nCells!=(int)intersectEdge2.size())
510 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
511 const DataArrayDouble *coo2(mesh1D->getCoords());
512 const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
513 const double *coo2Ptr(coo2->begin());
514 int offset1(coords1->getNumberOfTuples());
515 int offset2(offset1+coo2->getNumberOfTuples());
516 int offset3(offset2+addCoo.size()/2);
517 std::vector<double> addCooQuad;
518 MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
519 int tmp[4],cicnt(0),kk(0);
520 for(int i=0;i<nCells;i++)
522 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
523 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
524 const std::vector<int>& subEdges(intersectEdge2[i]);
525 int nbSubEdge(subEdges.size()/2);
526 for(int j=0;j<nbSubEdge;j++,kk++)
528 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));
529 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
530 INTERP_KERNEL::Edge *e2Ptr(e2);
531 std::map<int,int>::const_iterator itm;
532 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
534 tmp[0]=INTERP_KERNEL::NORM_SEG3;
535 itm=mergedNodes.find(subEdges[2*j]);
536 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
537 itm=mergedNodes.find(subEdges[2*j+1]);
538 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
539 tmp[3]=offset3+(int)addCooQuad.size()/2;
541 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
543 cOut->insertAtTheEnd(tmp,tmp+4);
544 ciOut->pushBackSilent(cicnt);
548 tmp[0]=INTERP_KERNEL::NORM_SEG2;
549 itm=mergedNodes.find(subEdges[2*j]);
550 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
551 itm=mergedNodes.find(subEdges[2*j+1]);
552 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
554 cOut->insertAtTheEnd(tmp,tmp+3);
555 ciOut->pushBackSilent(cicnt);
558 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
560 idsInRetColinear->pushBackSilent(kk);
561 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
566 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
567 ret->setConnectivity(cOut,ciOut,true);
568 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
569 arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
570 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,(int)addCooQuad.size()/2,2);
571 std::vector<const DataArrayDouble *> coordss(4);
572 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
573 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
578 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
580 std::vector<int> allEdges;
581 for(const int *it2(descBg);it2!=descEnd;it2++)
583 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
585 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
587 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
589 std::size_t nb(allEdges.size());
591 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
592 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
593 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
594 ret->setCoords(coords);
595 ret->allocateCells(1);
596 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
597 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
598 connOut[kk]=allEdges[2*kk];
599 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
603 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
605 const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
606 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
608 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
609 if(sz!=std::distance(descBg,descEnd))
610 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
611 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
612 std::vector<int> allEdges,centers;
613 const double *coordsPtr(coords->begin());
614 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
615 int offset(coords->getNumberOfTuples());
616 for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
618 INTERP_KERNEL::NormalizedCellType typeOfSon;
619 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
620 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
622 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
624 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
626 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
628 {//the current edge has been subsplit -> create corresponding centers.
629 std::size_t nbOfCentersToAppend(edge1.size()/2);
630 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
631 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
632 std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
633 for(std::size_t k=0;k<nbOfCentersToAppend;k++)
636 const double *aa(coordsPtr+2*(*it3++));
637 const double *bb(coordsPtr+2*(*it3++));
638 ee->getMiddleOfPoints(aa,bb,tmpp);
639 addCoo->insertAtTheEnd(tmpp,tmpp+2);
640 centers.push_back(offset+k);
644 std::size_t nb(allEdges.size());
646 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
647 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
648 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
650 ret->setCoords(coords);
653 addCoo->rearrange(2);
654 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
655 ret->setCoords(addCoo);
657 ret->allocateCells(1);
658 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
659 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
660 connOut[kk]=allEdges[2*kk];
661 connOut.insert(connOut.end(),centers.begin(),centers.end());
662 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
667 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
670 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
672 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
674 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
675 if(!cm.isQuadratic())
676 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
678 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
681 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
684 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
686 const INTERP_KERNEL::Edge *ee(*it);
687 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
691 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
694 const double *coo(mesh2D->getCoords()->begin());
695 std::size_t sz(conn.size());
696 std::vector<double> addCoo;
697 std::vector<int> conn2(conn);
698 int offset(mesh2D->getNumberOfNodes());
699 for(std::size_t i=0;i<sz;i++)
702 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
703 addCoo.insert(addCoo.end(),tmp,tmp+2);
704 conn2.push_back(offset+(int)i);
706 mesh2D->getCoords()->rearrange(1);
707 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
708 mesh2D->getCoords()->rearrange(2);
709 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
714 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
716 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
717 * a set of edges defined in \a splitMesh1D.
719 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
720 std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
722 std::size_t nb(edge1Bis.size()/2);
723 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
724 int iEnd(splitMesh1D->getNumberOfCells());
726 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
728 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
729 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
730 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
733 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
734 out0.resize(1); out1.resize(1);
735 std::vector<int>& connOut(out0[0]);
736 connOut.resize(nbOfEdgesOf2DCellSplit);
737 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
738 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
739 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
741 connOut[kk]=edge1Bis[2*kk];
742 edgesPtr[kk]=edge1BisPtr[2*kk];
747 // [i,iEnd[ contains the
748 out0.resize(2); out1.resize(2);
749 std::vector<int>& connOutLeft(out0[0]);
750 std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
751 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
752 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
753 for(std::size_t k=ii;k<jj+1;k++)
754 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
755 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
756 for(int ik=0;ik<iEnd;ik++)
758 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
759 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
762 for(int ik=iEnd-1;ik>=0;ik--)
763 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
764 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
765 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
766 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
767 for(int ik=0;ik<iEnd;ik++)
768 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
769 eright.insert(eright.end(),ees.begin(),ees.end());
777 CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
779 std::vector<int> _edges;
780 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
783 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
785 std::size_t nbe(edges.size());
786 std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
787 for(std::size_t i=0;i<nbe;i++)
789 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
790 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
792 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
793 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
794 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
800 EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
801 EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
802 bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
803 void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
804 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
808 MCAuto<MEDCouplingUMesh> _mesh;
809 MCAuto<INTERP_KERNEL::Edge> _edge;
814 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
816 const MEDCouplingUMesh *mesh(_mesh);
822 { _left++; _right++; return ; }
825 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
826 if((isLeft && isRight) || (!isLeft && !isRight))
827 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
838 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
839 if((isLeft && isRight) || (!isLeft && !isRight))
840 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
855 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
857 const MEDCouplingUMesh *mesh(_mesh);
860 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
863 {// not fully splitting cell case
864 if(mesh2D->getNumberOfCells()==1)
865 {//little optimization. 1 cell no need to find in which cell mesh is !
866 neighbors[0]=offset; neighbors[1]=offset;
871 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
872 int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
874 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
875 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
880 class VectorOfCellInfo
883 VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
884 std::size_t size() const { return _pool.size(); }
885 int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
886 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);
887 const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
888 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
889 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
890 void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
892 int getZePosOfEdgeGivenItsGlobalId(int pos) const;
893 void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
894 const CellInfo& get(int pos) const;
895 CellInfo& get(int pos);
897 std::vector<CellInfo> _pool;
898 MCAuto<MEDCouplingUMesh> _ze_mesh;
899 std::vector<EdgeInfo> _edge_info;
902 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
904 _pool[0]._edges=edges;
905 _pool[0]._edges_ptr=edgesPtr;
908 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
911 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
914 const MEDCouplingUMesh *zeMesh(_ze_mesh);
916 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
917 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
918 return zeMesh->getCellContainingPoint(barys->begin(),eps);
921 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
922 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
923 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
925 get(pos);//to check pos
926 bool isFast(pos==0 && _pool.size()==1);
927 std::size_t sz(edges.size());
928 // dealing with edges
930 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
932 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
934 std::vector<CellInfo> pool(_pool.size()-1+sz);
935 for(std::size_t i=0;i<pos;i++)
937 for(std::size_t j=0;j<sz;j++)
938 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
939 for(int i=pos+1;i<(int)_pool.size();i++)
940 pool[i+sz-1]=_pool[i];
944 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
952 std::vector< MCAuto<MEDCouplingUMesh> > ms;
955 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
959 if(pos<_ze_mesh->getNumberOfCells()-1)
961 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
964 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
965 for(std::size_t j=0;j<ms2.size();j++)
967 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
970 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
972 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
975 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
978 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
980 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
982 if((*it).isInMyRange(pos))
985 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
988 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
990 get(pos);//to perform the sanity check;
991 if(_edge_info.empty())
993 std::size_t sz(_edge_info.size()-1);
994 for(std::size_t i=0;i<sz;i++)
995 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
998 const CellInfo& VectorOfCellInfo::get(int pos) const
1000 if(pos<0 || pos>=(int)_pool.size())
1001 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1005 CellInfo& VectorOfCellInfo::get(int pos)
1007 if(pos<0 || pos>=(int)_pool.size())
1008 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1014 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1015 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1017 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1019 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1021 * \param [in] allEdges a list of pairs (beginNode, endNode). Represents all edges (already cut) in the single 2D cell being handled here. Linked with \a allEdgesPtr to get the equation of edge.
1023 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1024 MCAuto<DataArrayInt>& idsLeftRight)
1026 int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1027 if(nbCellsInSplitMesh1D==0)
1028 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1029 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1030 std::size_t nb(allEdges.size()),jj;
1032 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1033 std::vector<int> edge1Bis(nb*2);
1034 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1035 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1036 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1037 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1038 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1040 idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1041 int *idsLeftRightPtr(idsLeftRight->getPointer());
1042 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1044 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1045 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1046 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1047 // of the connectivity.
1048 MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1049 renumb->alloc(nbCellsInSplitMesh1D,1);
1050 const int * renumbP(renumb->begin());
1052 int i, first=cSplitPtr[1];
1053 // Follow 1D line backward as long as it is connected:
1054 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1055 first=cSplitPtr[ciSplitPtr[i]+1];
1056 if (i < nbCellsInSplitMesh1D-1)
1058 // Build circular permutation to shift consecutive edges together
1060 renumb->applyModulus(nbCellsInSplitMesh1D);
1061 splitMesh1D->renumberCells(renumbP, false);
1062 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1063 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1069 for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1070 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1072 for(;iEnd<nbCellsInSplitMesh1D;)
1074 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1080 if(iEnd<nbCellsInSplitMesh1D)
1083 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1084 int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1086 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1087 retTmp->setCoords(splitMesh1D->getCoords());
1088 retTmp->allocateCells();
1090 std::vector< std::vector<int> > out0;
1091 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1093 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1094 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1095 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1096 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1100 for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1101 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1103 return pool.getZeMesh().retn();
1107 * splitMesh1D is an input parameter but might have its cells renumbered.
1109 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1110 const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1111 MCAuto<DataArrayInt>& idsLeftRight)
1113 const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1115 std::vector<int> allEdges;
1116 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1117 for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1119 int edgeId(std::abs(*it)-1);
1120 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1121 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1122 const std::vector<int>& edge1(intersectEdge1[edgeId]);
1124 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1126 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1127 std::size_t sz(edge1.size());
1128 for(std::size_t cnt=0;cnt<sz;cnt++)
1129 allEdgesPtr.push_back(ee);
1132 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1135 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1137 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1138 {//easy case comparison not
1139 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1141 else if(typ1.isQuadratic() && typ2.isQuadratic())
1143 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1146 if(conn1[2]==conn2[2])
1148 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1149 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1153 {//only one is quadratic
1154 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1157 const double *a(0),*bb(0),*be(0);
1158 if(typ1.isQuadratic())
1160 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1164 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1166 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1167 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1173 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1174 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1176 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1178 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1180 if(candidatesIn2DEnd==candidatesIn2DBg)
1181 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1182 const double *coo(mesh2DSplit->getCoords()->begin());
1183 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1184 return *candidatesIn2DBg;
1185 int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1186 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1187 if(cellIdInMesh1DSplitRelative<0)
1188 cur1D->changeOrientationOfCells();
1189 const int *c1D(cur1D->getNodalConnectivity()->begin());
1190 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1191 for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1193 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1194 const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1195 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1196 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1197 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1198 for(unsigned it2=0;it2<sz;it2++)
1200 INTERP_KERNEL::NormalizedCellType typeOfSon;
1201 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1202 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1203 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1207 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1211 * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of int given resp start and stop.
1212 * So for all edge \a i in \a m1Desc \a intersectEdge1[i] is of length 2*n where n is the number of sub edges.
1213 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1214 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1215 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1216 * \param [out] addCoo - nodes to be append at the end
1217 * \param [out] mergedNodes - gives all pair of nodes of \a m2Desc that have same location than some nodes in \a m1Desc. key is id in \a m2Desc offsetted and value is id in \a m1Desc.
1219 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1220 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)
1222 static const int SPACEDIM=2;
1223 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1224 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1225 // Build BB tree of all edges in the tool mesh (second mesh)
1226 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1227 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1228 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1229 intersectEdge1.resize(nDescCell1);
1230 colinear2.resize(nDescCell2);
1231 subDiv2.resize(nDescCell2);
1232 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1234 std::vector<int> candidates1(1);
1235 int offset1(m1Desc->getNumberOfNodes());
1236 int offset2(offset1+m2Desc->getNumberOfNodes());
1237 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1239 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1240 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1241 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1243 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1244 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1245 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1247 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1248 // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one
1249 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1250 std::set<INTERP_KERNEL::Node *> nodes;
1251 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1252 std::size_t szz(nodes.size());
1253 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1254 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1255 for(std::size_t iii=0;iii<szz;iii++,itt++)
1256 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1257 // end of protection
1258 // Performs edge cutting:
1259 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1264 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1265 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1271 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1272 * It builds the descending connectivity of the two meshes, and then using a binary tree
1273 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1274 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1276 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1277 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1278 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1279 std::vector<double>& addCoo,
1280 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1282 // Build desc connectivity
1283 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1284 desc2=DataArrayInt::New();
1285 descIndx2=DataArrayInt::New();
1286 revDesc2=DataArrayInt::New();
1287 revDescIndx2=DataArrayInt::New();
1288 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1289 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1290 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1291 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1292 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1293 std::map<int,int> notUsedMap;
1294 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1295 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1296 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1300 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1301 * (newly created) nodes corresponding to the edge intersections.
1303 * @param[out] cr, crI connectivity of the resulting mesh
1304 * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1305 * TODO: describe input parameters
1307 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1308 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1309 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1310 const std::vector<double>& addCoords,
1311 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1313 static const int SPACEDIM=2;
1314 const double *coo1(m1->getCoords()->begin());
1315 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1316 int offset1(m1->getNumberOfNodes());
1317 const double *coo2(m2->getCoords()->begin());
1318 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1319 int offset2(offset1+m2->getNumberOfNodes());
1320 int offset3(offset2+((int)addCoords.size())/2);
1321 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1322 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1323 // Here a BBTree on 2D-cells, not on segments:
1324 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1325 int ncell1(m1->getNumberOfCells());
1327 for(int i=0;i<ncell1;i++)
1329 std::vector<int> candidates2;
1330 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1331 std::map<INTERP_KERNEL::Node *,int> mapp;
1332 std::map<int,INTERP_KERNEL::Node *> mappRev;
1333 INTERP_KERNEL::QuadraticPolygon pol1;
1334 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1335 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1336 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1337 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1338 // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1339 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1340 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1342 std::set<INTERP_KERNEL::Edge *> edges1;// store all edges of pol1 that are NOT consumed by intersect cells. If any after iteration over candidates2 -> a part of pol1 should appear in result
1343 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1344 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1345 for(it1.first();!it1.finished();it1.next())
1346 edges1.insert(it1.current()->getPtr());
1348 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1349 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1351 // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1352 // Again all the additional intersecting nodes are there.
1353 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1355 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1356 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1357 // Complete mapping with elements coming from the current cell it2 in mesh2:
1358 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1359 // pol2 is the new QP in the final merged result.
1360 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1361 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1364 // Now rebuild intersected cells from all this:
1365 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1367 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1368 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1369 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1370 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1372 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1373 // by m2 but that we still want to keep in the final result.
1378 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1380 catch(INTERP_KERNEL::Exception& e)
1382 std::ostringstream oss; oss << "Error when computing residual of cell #" << i << " in source/m1 mesh ! Maybe the neighbours of this cell in mesh are not well connected !\n" << "The deep reason is the following : " << e.what();
1383 throw INTERP_KERNEL::Exception(oss.str());
1386 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1387 (*it).second->decrRef();
1391 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1393 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1396 std::size_t sz(conn.size());
1397 std::size_t found(std::numeric_limits<std::size_t>::max());
1398 for(std::size_t i=0;i<sz;i++)
1400 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1401 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]};
1402 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1403 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1404 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1406 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];
1407 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]);
1409 if(dotTest>eps && dotTest<1.-eps)
1415 if(found==std::numeric_limits<std::size_t>::max())
1416 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1417 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1420 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1422 std::size_t sz(conn.size());
1423 std::vector<int> *curPart(&part0);
1424 for(std::size_t i=0;i<sz;i++)
1426 int nextt(conn[(i+1)%sz]);
1427 (*curPart).push_back(nextt);
1428 if(nextt==pt0 || nextt==pt1)
1434 (*curPart).push_back(nextt);
1440 * this method method splits cur cells 3D Surf in sub cells 3DSurf using the previous subsplit. This method is the last one used to clip.
1442 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1443 const int *desc, const int *descIndx, const double *coords, double eps,
1444 std::vector<std::vector<int> >& res) const
1446 checkFullyDefined();
1447 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1448 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1449 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1450 int nbOfCells(getNumberOfCells());
1452 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1453 for(int i=0;i<nbOfCells;i++)
1455 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1456 for(int j=0;j<nbOfFaces;j++)
1458 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1459 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1460 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1461 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1462 INTERP_KERNEL::NormalizedCellType cmsId;
1463 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1464 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1465 if(p.first!=-1 && p.second!=-1)
1469 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1470 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1471 std::vector<int> elt1,elt2;
1472 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1473 res.push_back(elt1);
1474 res.push_back(elt2);
1486 * 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).
1488 * \sa MEDCouplingUMesh::split2DCells
1490 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1492 checkConnectivityFullyDefined();
1493 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1494 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1495 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1496 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1497 int prevPosOfCi(ciPtr[0]);
1498 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1500 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1501 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1502 for(int j=0;j<sz;j++)
1504 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1505 for(int k=0;k<sz2;k++)
1506 *cPtr++=subPtr[offset2+k];
1508 *cPtr++=oldConn[prevPosOfCi+j+2];
1511 prevPosOfCi=ciPtr[1];
1512 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1515 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1516 _nodal_connec->decrRef();
1517 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1522 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1524 * \return int - the number of new nodes created.
1525 * \sa MEDCouplingUMesh::split2DCells
1527 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1529 checkConsistencyLight();
1530 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1531 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1532 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1533 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1534 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1535 const double *oldCoordsPtr(getCoords()->begin());
1536 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1537 int prevPosOfCi(ciPtr[0]);
1538 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1540 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1541 for(int j=0;j<sz;j++)
1542 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1543 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1544 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1546 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1550 cPtr[1]=oldConn[prevPosOfCi+2+j];
1551 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1554 std::vector<INTERP_KERNEL::Node *> ns(3);
1555 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1556 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1557 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1558 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1559 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1561 cPtr[1]=subPtr[offset2+k];
1562 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1564 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1567 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1569 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1570 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1573 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1574 _nodal_connec->decrRef();
1575 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1576 addCoo->rearrange(2);
1577 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1579 return addCoo->getNumberOfTuples();
1586 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1587 * returns a result mesh constituted by polygons.
1588 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1589 * all nodes from m2.
1590 * The meshes should be in 2D space. In
1591 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1593 * \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
1594 * 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)
1595 * \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
1596 * 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)
1597 * \param [in] eps - precision used to detect coincident mesh entities.
1598 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1599 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1600 * this array using decrRef() as it is no more needed.
1601 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1602 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1603 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1604 * any cell of \a m2. The caller is to delete this array using decrRef() as
1605 * it is no more needed.
1606 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1607 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1608 * is no more needed.
1609 * \throw If the coordinates array is not set in any of the meshes.
1610 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1611 * \throw If any of the meshes is not a 2D mesh in 2D space.
1613 * \sa conformize2D, mergeNodes
1615 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1616 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1619 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1620 m1->checkFullyDefined();
1621 m2->checkFullyDefined();
1622 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1623 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1624 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1626 // Step 1: compute all edge intersections (new nodes)
1627 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1628 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1629 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1630 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1631 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1632 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1633 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1634 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1635 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1636 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1638 // Step 2: re-order newly created nodes according to the ordering found in m2
1639 std::vector< std::vector<int> > intersectEdge2;
1640 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1641 subDiv2.clear(); dd5=0; dd6=0;
1644 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1645 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1646 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1647 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1649 // Step 4: Prepare final result:
1650 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1651 addCooDa->alloc((int)(addCoo.size())/2,2);
1652 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1653 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1654 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1655 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1656 std::vector<const DataArrayDouble *> coordss(4);
1657 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1658 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1659 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1660 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1661 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1662 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1663 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1664 ret->setConnectivity(conn,connI,true);
1665 ret->setCoords(coo);
1666 cellNb1=c1.retn(); cellNb2=c2.retn();
1671 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1672 * 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
1673 * and finally, in case of quadratic polygon the centers of edges new nodes.
1674 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1676 * \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
1677 * 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)
1678 * \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
1679 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1680 * \param [in] eps - precision used to perform intersections and localization operations.
1681 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1682 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1683 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1684 * So this array has a number of tuples equal to the number of cells of \a splitMesh2D and a number of component equal to 1.
1685 * \param [out] cellIdInMesh1D - the array of pair that gives for each cell id \a i in \a splitMesh1D the cell in \a splitMesh2D on the left for the 1st component
1686 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1687 * So this array has a number of tuples equal to the number of cells of \a splitMesh1D and a number of components equal to 2.
1689 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1691 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1693 if(!mesh2D || !mesh1D)
1694 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1695 mesh2D->checkFullyDefined();
1696 mesh1D->checkFullyDefined();
1697 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1698 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1699 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1700 // Step 1: compute all edge intersections (new nodes)
1701 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1702 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1703 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1705 // Build desc connectivity
1706 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1707 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1708 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1709 std::map<int,int> mergedNodes;
1710 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1711 // use mergeNodes to fix intersectEdge1
1712 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1714 std::size_t n((*it0).size()/2);
1715 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1716 std::map<int,int>::const_iterator it1;
1717 it1=mergedNodes.find(eltStart);
1718 if(it1!=mergedNodes.end())
1719 (*it0)[0]=(*it1).second;
1720 it1=mergedNodes.find(eltEnd);
1721 if(it1!=mergedNodes.end())
1722 (*it0)[2*n-1]=(*it1).second;
1725 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1726 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
1727 // Step 2: re-order newly created nodes according to the ordering found in m2
1728 std::vector< std::vector<int> > intersectEdge2;
1729 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1731 // Step 3: compute splitMesh1D
1732 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1733 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1734 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1735 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1736 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1737 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1738 // deal with cells in mesh2D that are not cut but only some of their edges are
1739 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1740 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1741 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1742 MCAuto<DataArrayInt> out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells
1743 if(!idsInDesc2DToBeRefined->empty())
1745 DataArrayInt *out0(0),*outi0(0);
1746 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1747 MCAuto<DataArrayInt> outi0s(outi0);
1749 out0s=out0s->buildUnique();
1753 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1754 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1755 MCAuto<DataArrayInt> elts,eltsIndex;
1756 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1757 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1758 if (eltsIndex->getNumberOfTuples() > 1)
1759 eltsIndex2 = eltsIndex->deltaShiftIndex();
1760 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1761 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1762 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1763 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1764 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1765 if((DataArrayInt *)out0s)
1766 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1767 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1768 // OK all is ready to insert in ret2 mesh
1769 if(!untouchedCells->empty())
1770 {// the most easy part, cells in mesh2D not impacted at all
1771 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1772 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1773 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1775 if((DataArrayInt *)out0s)
1776 {// here dealing with cells in out0s but not in cellsToBeModified
1777 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1778 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1779 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1781 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1782 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1784 int offset(ret2->getNumberOfTuples());
1785 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1786 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1787 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1788 int kk(0),*ret3ptr(partOfRet3->getPointer());
1789 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1791 int faceId(std::abs(*it)-1);
1792 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1794 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1797 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1798 ret3ptr[2*kk]=tmp+offset;
1799 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1800 ret3ptr[2*kk+1]=tmp+offset;
1803 {//the current edge is shared by a 2D cell that will be split just after
1804 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1805 ret3ptr[2*kk]=-(*it2+1);
1806 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1807 ret3ptr[2*kk+1]=-(*it2+1);
1811 m1Desc->setCoords(ret1->getCoords());
1812 ret1NonCol->setCoords(ret1->getCoords());
1813 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1814 if(!outMesh2DSplit.empty())
1816 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1817 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1818 (*itt)->setCoords(da);
1821 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1822 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1824 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1825 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1826 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1827 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1828 MCAuto<DataArrayInt> partOfRet3;
1829 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));
1830 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1831 outMesh2DSplit.push_back(splitOfOneCell);
1832 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1833 ret2->pushBackSilent(*it);
1836 std::size_t nbOfMeshes(outMesh2DSplit.size());
1837 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1838 for(std::size_t i=0;i<nbOfMeshes;i++)
1839 tmp[i]=outMesh2DSplit[i];
1841 ret1->getCoords()->setInfoOnComponents(compNames);
1842 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1843 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1845 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1846 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1848 int old2DCellId(-ret3->getIJ(*it,0)-1);
1849 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1850 ret3->setIJ(*it,0,FindRightCandidateAmong(ret2D,candidates->begin(),candidates->end(),ret1,*it%2==0?-((*it)/2+1):(*it)/2+1,eps));// div by 2 because 2 components natively in ret3
1852 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1855 splitMesh1D=ret1.retn();
1856 splitMesh2D=ret2D.retn();
1857 cellIdInMesh2D=ret2.retn();
1858 cellIdInMesh1D=ret3.retn();
1862 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1863 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1864 * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this method
1865 * will suppress such edges to use sub edges in \a this. So this method does not add nodes in \a this if merged edges are both linear (INTERP_KERNEL::NORM_SEG2).
1866 * In the other cases new nodes can be created. If any are created, they will be appended at the end of the coordinates object before the invocation of this method.
1868 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1869 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1871 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1872 * This method expects that all nodes in \a this are not closer than \a eps.
1873 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1875 * \param [in] eps the relative error to detect merged edges.
1876 * \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
1877 * that the user is expected to deal with.
1879 * \throw If \a this is not coherent.
1880 * \throw If \a this has not spaceDim equal to 2.
1881 * \throw If \a this has not meshDim equal to 2.
1882 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1884 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1886 static const int SPACEDIM=2;
1887 checkConsistencyLight();
1888 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1889 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1890 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1891 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1892 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1893 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1894 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1895 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1896 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1897 std::vector<double> addCoo;
1898 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1899 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1900 for(int i=0;i<nDescCell;i++)
1902 std::vector<int> candidates;
1903 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1904 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1905 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1907 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1908 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1909 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1910 INTERP_KERNEL::MergePoints merge;
1911 INTERP_KERNEL::QuadraticPolygon c1,c2;
1912 e1->intersectWith(e2,merge,c1,c2);
1913 e1->decrRef(); e2->decrRef();
1914 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1915 overlapEdge[i].push_back(*it);
1916 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1917 overlapEdge[*it].push_back(i);
1920 // splitting done. sort intersect point in intersectEdge.
1921 std::vector< std::vector<int> > middle(nDescCell);
1922 int nbOf2DCellsToBeSplit(0);
1923 bool middleNeedsToBeUsed(false);
1924 std::vector<bool> cells2DToTreat(nDescCell,false);
1925 for(int i=0;i<nDescCell;i++)
1927 std::vector<int>& isect(intersectEdge[i]);
1928 int sz((int)isect.size());
1931 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1932 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1933 e->sortSubNodesAbs(coords,isect);
1938 int idx0(rdi[i]),idx1(rdi[i+1]);
1940 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1941 if(!cells2DToTreat[rd[idx0]])
1943 cells2DToTreat[rd[idx0]]=true;
1944 nbOf2DCellsToBeSplit++;
1946 // try to reuse at most eventual 'middle' of SEG3
1947 std::vector<int>& mid(middle[i]);
1948 mid.resize(sz+1,-1);
1949 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1951 middleNeedsToBeUsed=true;
1952 const std::vector<int>& candidates(overlapEdge[i]);
1953 std::vector<int> trueCandidates;
1954 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1955 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1956 trueCandidates.push_back(*itc);
1957 int stNode(c[ci[i]+1]),endNode(isect[0]);
1958 for(int j=0;j<sz+1;j++)
1960 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1962 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1963 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1964 { mid[j]=*itc; break; }
1967 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1972 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1973 if(nbOf2DCellsToBeSplit==0)
1976 int *retPtr(ret->getPointer());
1977 for(int i=0;i<nCell;i++)
1978 if(cells2DToTreat[i])
1981 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1982 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1983 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1984 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1985 if(middleNeedsToBeUsed)
1986 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1987 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1988 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1989 setCoords(modif->getCoords());//if nbOfNodesCreated==0 modif and this have the same coordinates pointer so this line has no effect. But for quadratic cases this line is important.
1990 setPartOfMySelf(ret->begin(),ret->end(),*modif);
1992 bool areNodesMerged; int newNbOfNodes;
1993 if(nbOfNodesCreated!=0)
1994 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2000 * This non const method works on 2D mesh. This method scans every cell in \a this and look if each edge constituting this cell is not mergeable with neighbors edges of that cell.
2001 * If yes, the cell is "repaired" to minimize at most its number of edges. So this method do not change the overall shape of cells in \a this (with eps precision).
2002 * This method do not take care of shared edges between cells, so this method can lead to a non conform mesh (\a this). If a conform mesh is required you're expected
2003 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2004 * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automatically. For 2D "repaired" quadratic cells
2005 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2007 * If the returned array is empty it means that nothing has changed in \a this (as if it were a const method). If the array is not empty the connectivity of \a this is modified
2008 * using new instance, idem for coordinates.
2010 * If \a this is constituted by only linear 2D cells, this method is close to the computation of the convex hull of each cells in \a this.
2012 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
2014 * \throw If \a this is not coherent.
2015 * \throw If \a this has not spaceDim equal to 2.
2016 * \throw If \a this has not meshDim equal to 2.
2018 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2020 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2022 return internalColinearize2D(eps, false);
2026 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2027 * In a given 2D cell, if two edges are colinear and the junction point is used by a third edge, the two edges will not be
2028 * merged, contrary to colinearize2D().
2030 * \sa MEDCouplingUMesh::colinearize2D
2032 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2034 return internalColinearize2D(eps, true);
2039 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2041 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2043 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2044 checkConsistencyLight();
2045 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2046 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2047 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2048 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2049 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2050 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2051 std::map<int, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2054 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2055 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2056 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2057 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2058 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2059 MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2060 MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2061 const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2062 for(auto it = ids->begin(); it != ids->end(); it++)
2063 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2066 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2067 const double *coords(_coords->begin());
2068 int *newciptr(newci->getPointer());
2069 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2071 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2072 ret->pushBackSilent(i);
2073 newciptr[1]=newc->getNumberOfTuples();
2078 if(!appendedCoords->empty())
2080 appendedCoords->rearrange(2);
2081 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2083 setCoords(newCoords);
2086 setConnectivity(newc,newci,true);
2092 * c, cI describe a wire mesh in 3D space, in local numbering
2093 * startNode, endNode in global numbering
2094 *\return true if the segment is indeed split
2096 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2097 const int * c, const int * cI, const int *idsBg, const int *endBg,
2098 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2100 using namespace std;
2102 const int SPACEDIM=3;
2103 typedef pair<double, int> PairDI;
2105 for (const int * it = idsBg; it != endBg; ++it)
2107 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2108 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2109 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2110 x.insert(make_pair(coo[end*SPACEDIM], end));
2113 vector<PairDI> xx(x.begin(), x.end());
2114 sort(xx.begin(),xx.end());
2115 pointIds.reserve(xx.size());
2117 // Keep what is inside [startNode, endNode]:
2119 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2121 const int idx = (*it).second;
2124 if (idx == startNode) go = 1;
2125 if (idx == endNode) go = 2;
2126 if (go) pointIds.push_back(idx);
2129 pointIds.push_back(idx);
2130 if (idx == endNode || idx == startNode)
2134 // vector<int> pointIds2(pointIds.size()+2);
2135 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2136 // pointIds2[0] = startNode;
2137 // pointIds2[pointIds2.size()-1] = endNode;
2140 reverse(pointIds.begin(), pointIds.end());
2142 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2143 for (const int * it = idsBg; it != endBg; ++it)
2145 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2146 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2147 if (itStart == pointIds.end()) continue;
2148 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2149 if (itEnd == pointIds.end()) continue;
2150 if (abs(distance(itEnd, itStart)) != 1) continue;
2151 hitSegs.push_back(*it); // segment is undivided.
2154 return (pointIds.size() > 2); // something else apart start and end node
2157 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2158 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2160 using namespace std;
2161 int dst = distance(sIdxConn, sIdxConnE);
2162 modifiedFace.reserve(dst + insidePoints.size()-2);
2163 modifiedFace.resize(dst);
2164 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2166 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2167 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2168 if (startPos == shortEnd)
2169 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2170 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2171 if (endPos == shortEnd)
2172 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2173 int d = distance(startPos, endPos);
2174 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2175 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2177 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2184 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2185 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2186 * This method performs a conformization of \b this.
2188 * Only polyhedron cells are supported. You can call convertAllToPoly()
2190 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2191 * This method expects that all nodes in \a this are not closer than \a eps.
2192 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2194 * \param [in] eps the relative error to detect merged edges.
2195 * \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
2196 * that the user is expected to deal with.
2198 * \throw If \a this is not coherent.
2199 * \throw If \a this has not spaceDim equal to 3.
2200 * \throw If \a this has not meshDim equal to 3.
2201 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2203 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2205 using namespace std;
2207 static const int SPACEDIM=3;
2208 checkConsistencyLight();
2209 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2210 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2211 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2212 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2214 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2215 const double * coo(_coords->begin());
2216 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2219 /*************************
2221 *************************/
2222 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2223 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2224 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2225 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2226 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2229 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2230 const double *bbox(bboxArr->begin()); getCoords()->begin();
2231 int nDescCell(mDesc->getNumberOfCells());
2232 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2233 // Surfaces - handle biggest first
2234 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2235 DataArrayDouble * surfs = surfF->getArray();
2237 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2238 DataArrayDouble * normals = normalsF->getArray();
2239 const double * normalsP = normals->getConstPointer();
2241 // Sort faces by decreasing surface:
2242 vector< pair<double,int> > S;
2243 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2245 pair<double,int> p = make_pair(surfs->begin()[i], i);
2248 sort(S.rbegin(),S.rend()); // reverse sort
2249 vector<bool> hit(nDescCell);
2250 fill(hit.begin(), hit.end(), false);
2251 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2253 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2255 int faceIdx = (*it).second;
2256 if (hit[faceIdx]) continue;
2258 vector<int> candidates, cands2;
2259 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2260 // Keep only candidates whose normal matches the normal of current face
2261 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2263 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2264 if (*it2 != faceIdx && col)
2265 cands2.push_back(*it2);
2270 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2271 INTERP_KERNEL::TranslationRotationMatrix rotation;
2272 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2273 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2274 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2276 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2277 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2278 double * cooPartRef(mPartRef->_coords->getPointer());
2279 double * cooPartCand(mPartCand->_coords->getPointer());
2280 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2281 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2282 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2283 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2285 // Localize faces in 2D thanks to barycenters
2286 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2287 vector<int> compo; compo.push_back(2);
2288 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2289 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2290 if (!idsGoodPlane->getNumberOfTuples())
2293 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2295 compo[0] = 0; compo.push_back(1);
2296 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2297 mPartRef->changeSpaceDimension(2,0.0);
2298 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2299 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2301 if (!cc->getNumberOfTuples())
2303 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2306 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2307 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2310 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2311 throw INTERP_KERNEL::Exception(oss.str());
2315 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2317 if (!ids->getNumberOfTuples())
2320 double checkSurf=0.0;
2321 const int * idsGoodPlaneP(idsGoodPlane->begin());
2322 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2324 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2325 hit[faceIdx2] = true;
2326 checkSurf += surfs->begin()[faceIdx2];
2328 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2331 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2332 throw INTERP_KERNEL::Exception(oss.str());
2335 // For all polyhedrons using this face, replace connectivity:
2336 vector<int> polyIndices, packsIds, facePack;
2337 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2338 polyIndices.push_back(revDescP[ii]);
2339 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2341 // Current face connectivity
2342 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2343 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2344 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2345 // Deletion of old faces
2347 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2349 if (packsIds[jj] == -1)
2350 // The below should never happen - if a face is used several times, with a different layout of the nodes
2351 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2352 // faces which are actually used only once, by a single cell. This is different for edges below.
2353 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2355 connSla->deletePack(*it2, packsIds[jj]);
2357 // Insertion of new faces:
2358 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2360 // Build pack from the face to insert:
2361 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2363 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2364 // Insert it in all hit polyhedrons:
2365 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2366 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2371 // Set back modified connectivity
2372 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2373 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2374 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2377 /************************
2379 ************************/
2380 // Now we have a face-conform mesh.
2382 // Recompute descending
2383 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2384 // Rebuild desc connectivity with orientation this time!!
2385 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2386 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2387 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2388 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2389 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2390 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2391 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2392 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2393 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2394 // std::cout << "writing!\n";
2395 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2396 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2397 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2398 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2399 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2400 const double *bbox2(bboxArr->begin());
2401 int nDesc2Cell=mDesc2->getNumberOfCells();
2402 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2404 // Edges - handle longest first
2405 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2406 DataArrayDouble * lens = lenF->getArray();
2408 // Sort edges by decreasing length:
2409 vector<pair<double,int> > S;
2410 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2412 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2415 sort(S.rbegin(),S.rend()); // reverse sort
2417 vector<bool> hit(nDesc2Cell);
2418 fill(hit.begin(), hit.end(), false);
2420 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2422 int eIdx = (*it).second;
2426 vector<int> candidates, cands2;
2427 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2428 // Keep only candidates colinear with current edge
2430 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2431 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2432 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2433 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2436 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2437 for (int i3=0; i3 < 3; i3++)
2438 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2439 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2440 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2441 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2443 cands2.push_back(*it2);
2445 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2448 // Now rotate edges to bring them on Ox
2449 int startNode = cDesc2[cIDesc2[eIdx]+1];
2450 int endNode = cDesc2[cIDesc2[eIdx]+2];
2451 INTERP_KERNEL::TranslationRotationMatrix rotation;
2452 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2453 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2454 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2455 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2458 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2459 nbElemsNotM1 = tmp->getNbOfElems();
2461 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2462 double * cooPartRef(mPartRef->_coords->getPointer());
2463 double * cooPartCand(mPartCand->_coords->getPointer());
2464 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2465 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2466 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2467 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2470 // Eliminate all edges for which y or z is not null
2471 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2472 vector<int> compo; compo.push_back(1);
2473 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2475 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2476 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2477 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2478 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2479 if (!idsGoodLine->getNumberOfTuples())
2482 // Now the ordering along the Ox axis:
2483 std::vector<int> insidePoints, hitSegs;
2484 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2485 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2486 idsGoodLine->begin(), idsGoodLine->end(),
2487 /*out*/insidePoints, hitSegs);
2488 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2489 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2490 hit[cands2[*its]] = true;
2492 if (!isSplit) // current segment remains in one piece
2495 // Get original node IDs in global coords array
2496 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2497 *iit = nodeMapInv->begin()[*iit];
2499 vector<int> polyIndices, packsIds, facePack;
2500 // For each face implying this edge
2501 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2503 int faceIdx = revDescP2[ii];
2504 // each cell where this face is involved connectivity will be modified:
2505 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2507 // Current face connectivity
2508 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2509 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2511 vector<int> modifiedFace;
2512 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2513 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2514 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2518 // Rebuild 3D connectivity from descending:
2519 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2520 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2521 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2522 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2523 newConn->set3(superIdx, idx, vals);
2524 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2525 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2527 int sz, faceIdx = abs(descP[jj])-1;
2528 bool orient = descP[jj]>0;
2529 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2531 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2534 vector<int> rev(sz-1);
2535 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2536 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2540 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2543 ret = ret->buildUniqueNotSorted();