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);
1363 // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1364 for (auto &p: pol2s)
1365 p.cleanDegeneratedConsecutiveEdges();
1366 edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges.
1369 // Now rebuild intersected cells from all this:
1370 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1372 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1373 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1374 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1375 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1377 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1378 // by m2 but that we still want to keep in the final result.
1383 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1385 catch(INTERP_KERNEL::Exception& e)
1387 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();
1388 throw INTERP_KERNEL::Exception(oss.str());
1391 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1392 (*it).second->decrRef();
1396 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1398 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1401 std::size_t sz(conn.size());
1402 std::size_t found(std::numeric_limits<std::size_t>::max());
1403 for(std::size_t i=0;i<sz;i++)
1405 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1406 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]};
1407 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1408 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1409 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1411 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];
1412 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]);
1414 if(dotTest>eps && dotTest<1.-eps)
1420 if(found==std::numeric_limits<std::size_t>::max())
1421 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1422 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1425 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1427 std::size_t sz(conn.size());
1428 std::vector<int> *curPart(&part0);
1429 for(std::size_t i=0;i<sz;i++)
1431 int nextt(conn[(i+1)%sz]);
1432 (*curPart).push_back(nextt);
1433 if(nextt==pt0 || nextt==pt1)
1439 (*curPart).push_back(nextt);
1445 * 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.
1447 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1448 const int *desc, const int *descIndx, const double *coords, double eps,
1449 std::vector<std::vector<int> >& res) const
1451 checkFullyDefined();
1452 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1453 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1454 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1455 int nbOfCells(getNumberOfCells());
1457 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1458 for(int i=0;i<nbOfCells;i++)
1460 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1461 for(int j=0;j<nbOfFaces;j++)
1463 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1464 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1465 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1466 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1467 INTERP_KERNEL::NormalizedCellType cmsId;
1468 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1469 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1470 if(p.first!=-1 && p.second!=-1)
1474 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1475 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1476 std::vector<int> elt1,elt2;
1477 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1478 res.push_back(elt1);
1479 res.push_back(elt2);
1491 * 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).
1493 * \sa MEDCouplingUMesh::split2DCells
1495 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1497 checkConnectivityFullyDefined();
1498 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1499 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1500 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1501 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1502 int prevPosOfCi(ciPtr[0]);
1503 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1505 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1506 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1507 for(int j=0;j<sz;j++)
1509 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1510 for(int k=0;k<sz2;k++)
1511 *cPtr++=subPtr[offset2+k];
1513 *cPtr++=oldConn[prevPosOfCi+j+2];
1516 prevPosOfCi=ciPtr[1];
1517 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1520 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1521 _nodal_connec->decrRef();
1522 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1527 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1529 * \return int - the number of new nodes created.
1530 * \sa MEDCouplingUMesh::split2DCells
1532 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1534 checkConsistencyLight();
1535 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1536 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1537 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1538 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1539 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1540 const double *oldCoordsPtr(getCoords()->begin());
1541 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1542 int prevPosOfCi(ciPtr[0]);
1543 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1545 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1546 for(int j=0;j<sz;j++)
1547 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1548 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1549 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1551 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1555 cPtr[1]=oldConn[prevPosOfCi+2+j];
1556 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1559 std::vector<INTERP_KERNEL::Node *> ns(3);
1560 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1561 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1562 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1563 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1564 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1566 cPtr[1]=subPtr[offset2+k];
1567 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1569 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1572 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1574 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1575 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1578 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1579 _nodal_connec->decrRef();
1580 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1581 addCoo->rearrange(2);
1582 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1584 return addCoo->getNumberOfTuples();
1591 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1592 * returns a result mesh constituted by polygons.
1593 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1594 * all nodes from m2.
1595 * The meshes should be in 2D space. In
1596 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1598 * \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
1599 * 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)
1600 * \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
1601 * 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)
1602 * \param [in] eps - precision used to detect coincident mesh entities.
1603 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1604 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1605 * this array using decrRef() as it is no more needed.
1606 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1607 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1608 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1609 * any cell of \a m2. The caller is to delete this array using decrRef() as
1610 * it is no more needed.
1611 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1612 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1613 * is no more needed.
1614 * \throw If the coordinates array is not set in any of the meshes.
1615 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1616 * \throw If any of the meshes is not a 2D mesh in 2D space.
1618 * \sa conformize2D, mergeNodes
1620 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1621 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1624 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1625 m1->checkFullyDefined();
1626 m2->checkFullyDefined();
1627 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1628 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1629 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1631 // Step 1: compute all edge intersections (new nodes)
1632 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1633 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1634 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1635 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1636 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1637 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1638 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1639 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1640 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1641 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1643 // Step 2: re-order newly created nodes according to the ordering found in m2
1644 std::vector< std::vector<int> > intersectEdge2;
1645 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1646 subDiv2.clear(); dd5=0; dd6=0;
1649 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1650 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1651 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1652 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1654 // Step 4: Prepare final result:
1655 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1656 addCooDa->alloc((int)(addCoo.size())/2,2);
1657 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1658 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1659 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1660 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1661 std::vector<const DataArrayDouble *> coordss(4);
1662 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1663 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1664 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1665 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1666 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1667 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1668 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1669 ret->setConnectivity(conn,connI,true);
1670 ret->setCoords(coo);
1671 cellNb1=c1.retn(); cellNb2=c2.retn();
1676 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1677 * 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
1678 * and finally, in case of quadratic polygon the centers of edges new nodes.
1679 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1681 * \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
1682 * 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)
1683 * \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
1684 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1685 * \param [in] eps - precision used to perform intersections and localization operations.
1686 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1687 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1688 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1689 * 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.
1690 * \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
1691 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1692 * 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.
1694 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1696 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1698 if(!mesh2D || !mesh1D)
1699 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1700 mesh2D->checkFullyDefined();
1701 mesh1D->checkFullyDefined();
1702 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1703 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1704 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1705 // Step 1: compute all edge intersections (new nodes)
1706 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1707 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1708 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1710 // Build desc connectivity
1711 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1712 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1713 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1714 std::map<int,int> mergedNodes;
1715 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1716 // use mergeNodes to fix intersectEdge1
1717 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1719 std::size_t n((*it0).size()/2);
1720 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1721 std::map<int,int>::const_iterator it1;
1722 it1=mergedNodes.find(eltStart);
1723 if(it1!=mergedNodes.end())
1724 (*it0)[0]=(*it1).second;
1725 it1=mergedNodes.find(eltEnd);
1726 if(it1!=mergedNodes.end())
1727 (*it0)[2*n-1]=(*it1).second;
1730 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1731 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
1732 // Step 2: re-order newly created nodes according to the ordering found in m2
1733 std::vector< std::vector<int> > intersectEdge2;
1734 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1736 // Step 3: compute splitMesh1D
1737 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1738 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1739 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1740 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1741 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1742 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1743 // deal with cells in mesh2D that are not cut but only some of their edges are
1744 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1745 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1746 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1747 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
1748 if(!idsInDesc2DToBeRefined->empty())
1750 DataArrayInt *out0(0),*outi0(0);
1751 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1752 MCAuto<DataArrayInt> outi0s(outi0);
1754 out0s=out0s->buildUnique();
1758 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1759 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1760 MCAuto<DataArrayInt> elts,eltsIndex;
1761 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1762 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1763 if (eltsIndex->getNumberOfTuples() > 1)
1764 eltsIndex2 = eltsIndex->deltaShiftIndex();
1765 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1766 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1767 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1768 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1769 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1770 if((DataArrayInt *)out0s)
1771 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1772 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1773 // OK all is ready to insert in ret2 mesh
1774 if(!untouchedCells->empty())
1775 {// the most easy part, cells in mesh2D not impacted at all
1776 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1777 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1778 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1780 if((DataArrayInt *)out0s)
1781 {// here dealing with cells in out0s but not in cellsToBeModified
1782 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1783 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1784 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1786 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1787 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1789 int offset(ret2->getNumberOfTuples());
1790 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1791 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1792 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1793 int kk(0),*ret3ptr(partOfRet3->getPointer());
1794 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1796 int faceId(std::abs(*it)-1);
1797 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1799 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1802 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1803 ret3ptr[2*kk]=tmp+offset;
1804 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1805 ret3ptr[2*kk+1]=tmp+offset;
1808 {//the current edge is shared by a 2D cell that will be split just after
1809 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1810 ret3ptr[2*kk]=-(*it2+1);
1811 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1812 ret3ptr[2*kk+1]=-(*it2+1);
1816 m1Desc->setCoords(ret1->getCoords());
1817 ret1NonCol->setCoords(ret1->getCoords());
1818 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1819 if(!outMesh2DSplit.empty())
1821 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1822 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1823 (*itt)->setCoords(da);
1826 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1827 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1829 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1830 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1831 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1832 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1833 MCAuto<DataArrayInt> partOfRet3;
1834 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));
1835 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1836 outMesh2DSplit.push_back(splitOfOneCell);
1837 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1838 ret2->pushBackSilent(*it);
1841 std::size_t nbOfMeshes(outMesh2DSplit.size());
1842 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1843 for(std::size_t i=0;i<nbOfMeshes;i++)
1844 tmp[i]=outMesh2DSplit[i];
1846 ret1->getCoords()->setInfoOnComponents(compNames);
1847 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1848 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1850 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1851 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1853 int old2DCellId(-ret3->getIJ(*it,0)-1);
1854 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1855 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
1857 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1860 splitMesh1D=ret1.retn();
1861 splitMesh2D=ret2D.retn();
1862 cellIdInMesh2D=ret2.retn();
1863 cellIdInMesh1D=ret3.retn();
1867 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1868 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1869 * 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
1870 * 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).
1871 * 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.
1873 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1874 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1876 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1877 * This method expects that all nodes in \a this are not closer than \a eps.
1878 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1880 * \param [in] eps the relative error to detect merged edges.
1881 * \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
1882 * that the user is expected to deal with.
1884 * \throw If \a this is not coherent.
1885 * \throw If \a this has not spaceDim equal to 2.
1886 * \throw If \a this has not meshDim equal to 2.
1887 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1889 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1891 static const int SPACEDIM=2;
1892 checkConsistencyLight();
1893 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1894 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1895 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1896 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1897 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1898 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1899 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1900 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1901 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1902 std::vector<double> addCoo;
1903 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1904 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1905 for(int i=0;i<nDescCell;i++)
1907 std::vector<int> candidates;
1908 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1909 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1910 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1912 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1913 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1914 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1915 INTERP_KERNEL::MergePoints merge;
1916 INTERP_KERNEL::QuadraticPolygon c1,c2;
1917 e1->intersectWith(e2,merge,c1,c2);
1918 e1->decrRef(); e2->decrRef();
1919 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1920 overlapEdge[i].push_back(*it);
1921 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1922 overlapEdge[*it].push_back(i);
1925 // splitting done. sort intersect point in intersectEdge.
1926 std::vector< std::vector<int> > middle(nDescCell);
1927 int nbOf2DCellsToBeSplit(0);
1928 bool middleNeedsToBeUsed(false);
1929 std::vector<bool> cells2DToTreat(nDescCell,false);
1930 for(int i=0;i<nDescCell;i++)
1932 std::vector<int>& isect(intersectEdge[i]);
1933 int sz((int)isect.size());
1936 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1937 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1938 e->sortSubNodesAbs(coords,isect);
1943 int idx0(rdi[i]),idx1(rdi[i+1]);
1945 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1946 if(!cells2DToTreat[rd[idx0]])
1948 cells2DToTreat[rd[idx0]]=true;
1949 nbOf2DCellsToBeSplit++;
1951 // try to reuse at most eventual 'middle' of SEG3
1952 std::vector<int>& mid(middle[i]);
1953 mid.resize(sz+1,-1);
1954 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1956 middleNeedsToBeUsed=true;
1957 const std::vector<int>& candidates(overlapEdge[i]);
1958 std::vector<int> trueCandidates;
1959 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1960 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1961 trueCandidates.push_back(*itc);
1962 int stNode(c[ci[i]+1]),endNode(isect[0]);
1963 for(int j=0;j<sz+1;j++)
1965 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1967 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1968 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1969 { mid[j]=*itc; break; }
1972 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1977 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1978 if(nbOf2DCellsToBeSplit==0)
1981 int *retPtr(ret->getPointer());
1982 for(int i=0;i<nCell;i++)
1983 if(cells2DToTreat[i])
1986 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1987 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1988 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1989 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1990 if(middleNeedsToBeUsed)
1991 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1992 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1993 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1994 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.
1995 setPartOfMySelf(ret->begin(),ret->end(),*modif);
1997 bool areNodesMerged; int newNbOfNodes;
1998 if(nbOfNodesCreated!=0)
1999 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2005 * 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.
2006 * 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).
2007 * 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
2008 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2009 * 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
2010 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2012 * 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
2013 * using new instance, idem for coordinates.
2015 * 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.
2017 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
2019 * \throw If \a this is not coherent.
2020 * \throw If \a this has not spaceDim equal to 2.
2021 * \throw If \a this has not meshDim equal to 2.
2023 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2025 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2027 return internalColinearize2D(eps, false);
2031 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2032 * 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
2033 * merged, contrary to colinearize2D().
2035 * \sa MEDCouplingUMesh::colinearize2D
2037 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2039 return internalColinearize2D(eps, true);
2044 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2046 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2048 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2049 checkConsistencyLight();
2050 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2051 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2052 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2053 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2054 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2055 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2056 std::map<int, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2059 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2060 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2061 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2062 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2063 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2064 MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2065 MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2066 const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2067 for(auto it = ids->begin(); it != ids->end(); it++)
2068 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2071 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2072 const double *coords(_coords->begin());
2073 int *newciptr(newci->getPointer());
2074 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2076 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2077 ret->pushBackSilent(i);
2078 newciptr[1]=newc->getNumberOfTuples();
2083 if(!appendedCoords->empty())
2085 appendedCoords->rearrange(2);
2086 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2088 setCoords(newCoords);
2091 setConnectivity(newc,newci,true);
2097 * c, cI describe a wire mesh in 3D space, in local numbering
2098 * startNode, endNode in global numbering
2099 *\return true if the segment is indeed split
2101 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2102 const int * c, const int * cI, const int *idsBg, const int *endBg,
2103 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2105 using namespace std;
2107 const int SPACEDIM=3;
2108 typedef pair<double, int> PairDI;
2110 for (const int * it = idsBg; it != endBg; ++it)
2112 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2113 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2114 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2115 x.insert(make_pair(coo[end*SPACEDIM], end));
2118 vector<PairDI> xx(x.begin(), x.end());
2119 sort(xx.begin(),xx.end());
2120 pointIds.reserve(xx.size());
2122 // Keep what is inside [startNode, endNode]:
2124 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2126 const int idx = (*it).second;
2129 if (idx == startNode) go = 1;
2130 if (idx == endNode) go = 2;
2131 if (go) pointIds.push_back(idx);
2134 pointIds.push_back(idx);
2135 if (idx == endNode || idx == startNode)
2139 // vector<int> pointIds2(pointIds.size()+2);
2140 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2141 // pointIds2[0] = startNode;
2142 // pointIds2[pointIds2.size()-1] = endNode;
2145 reverse(pointIds.begin(), pointIds.end());
2147 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2148 for (const int * it = idsBg; it != endBg; ++it)
2150 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2151 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2152 if (itStart == pointIds.end()) continue;
2153 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2154 if (itEnd == pointIds.end()) continue;
2155 if (abs(distance(itEnd, itStart)) != 1) continue;
2156 hitSegs.push_back(*it); // segment is undivided.
2159 return (pointIds.size() > 2); // something else apart start and end node
2162 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2163 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2165 using namespace std;
2166 int dst = distance(sIdxConn, sIdxConnE);
2167 modifiedFace.reserve(dst + insidePoints.size()-2);
2168 modifiedFace.resize(dst);
2169 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2171 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2172 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2173 if (startPos == shortEnd)
2174 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2175 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2176 if (endPos == shortEnd)
2177 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2178 int d = distance(startPos, endPos);
2179 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2180 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2182 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2189 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2190 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2191 * This method performs a conformization of \b this.
2193 * Only polyhedron cells are supported. You can call convertAllToPoly()
2195 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2196 * This method expects that all nodes in \a this are not closer than \a eps.
2197 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2199 * \param [in] eps the relative error to detect merged edges.
2200 * \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
2201 * that the user is expected to deal with.
2203 * \throw If \a this is not coherent.
2204 * \throw If \a this has not spaceDim equal to 3.
2205 * \throw If \a this has not meshDim equal to 3.
2206 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2208 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2210 using namespace std;
2212 static const int SPACEDIM=3;
2213 checkConsistencyLight();
2214 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2215 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2216 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2217 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2219 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2220 const double * coo(_coords->begin());
2221 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2224 /*************************
2226 *************************/
2227 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2228 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2229 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2230 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2231 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2234 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2235 const double *bbox(bboxArr->begin()); getCoords()->begin();
2236 int nDescCell(mDesc->getNumberOfCells());
2237 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2238 // Surfaces - handle biggest first
2239 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2240 DataArrayDouble * surfs = surfF->getArray();
2242 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2243 DataArrayDouble * normals = normalsF->getArray();
2244 const double * normalsP = normals->getConstPointer();
2246 // Sort faces by decreasing surface:
2247 vector< pair<double,int> > S;
2248 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2250 pair<double,int> p = make_pair(surfs->begin()[i], i);
2253 sort(S.rbegin(),S.rend()); // reverse sort
2254 vector<bool> hit(nDescCell);
2255 fill(hit.begin(), hit.end(), false);
2256 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2258 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2260 int faceIdx = (*it).second;
2261 if (hit[faceIdx]) continue;
2263 vector<int> candidates, cands2;
2264 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2265 // Keep only candidates whose normal matches the normal of current face
2266 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2268 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2269 if (*it2 != faceIdx && col)
2270 cands2.push_back(*it2);
2275 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2276 INTERP_KERNEL::TranslationRotationMatrix rotation;
2277 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2278 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2279 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2281 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2282 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2283 double * cooPartRef(mPartRef->_coords->getPointer());
2284 double * cooPartCand(mPartCand->_coords->getPointer());
2285 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2286 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2287 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2288 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2290 // Localize faces in 2D thanks to barycenters
2291 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2292 vector<int> compo; compo.push_back(2);
2293 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2294 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2295 if (!idsGoodPlane->getNumberOfTuples())
2298 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2300 compo[0] = 0; compo.push_back(1);
2301 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2302 mPartRef->changeSpaceDimension(2,0.0);
2303 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2304 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2306 if (!cc->getNumberOfTuples())
2308 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2311 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2312 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2315 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2316 throw INTERP_KERNEL::Exception(oss.str());
2320 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2322 if (!ids->getNumberOfTuples())
2325 double checkSurf=0.0;
2326 const int * idsGoodPlaneP(idsGoodPlane->begin());
2327 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2329 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2330 hit[faceIdx2] = true;
2331 checkSurf += surfs->begin()[faceIdx2];
2333 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2336 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2337 throw INTERP_KERNEL::Exception(oss.str());
2340 // For all polyhedrons using this face, replace connectivity:
2341 vector<int> polyIndices, packsIds, facePack;
2342 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2343 polyIndices.push_back(revDescP[ii]);
2344 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2346 // Current face connectivity
2347 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2348 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2349 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2350 // Deletion of old faces
2352 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2354 if (packsIds[jj] == -1)
2355 // The below should never happen - if a face is used several times, with a different layout of the nodes
2356 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2357 // faces which are actually used only once, by a single cell. This is different for edges below.
2358 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2360 connSla->deletePack(*it2, packsIds[jj]);
2362 // Insertion of new faces:
2363 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2365 // Build pack from the face to insert:
2366 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2368 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2369 // Insert it in all hit polyhedrons:
2370 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2371 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2376 // Set back modified connectivity
2377 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2378 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2379 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2382 /************************
2384 ************************/
2385 // Now we have a face-conform mesh.
2387 // Recompute descending
2388 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2389 // Rebuild desc connectivity with orientation this time!!
2390 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2391 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2392 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2393 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2394 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2395 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2396 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2397 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2398 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2399 // std::cout << "writing!\n";
2400 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2401 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2402 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2403 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2404 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2405 const double *bbox2(bboxArr->begin());
2406 int nDesc2Cell=mDesc2->getNumberOfCells();
2407 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2409 // Edges - handle longest first
2410 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2411 DataArrayDouble * lens = lenF->getArray();
2413 // Sort edges by decreasing length:
2414 vector<pair<double,int> > S;
2415 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2417 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2420 sort(S.rbegin(),S.rend()); // reverse sort
2422 vector<bool> hit(nDesc2Cell);
2423 fill(hit.begin(), hit.end(), false);
2425 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2427 int eIdx = (*it).second;
2431 vector<int> candidates, cands2;
2432 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2433 // Keep only candidates colinear with current edge
2435 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2436 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2437 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2438 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2441 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2442 for (int i3=0; i3 < 3; i3++)
2443 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2444 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2445 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2446 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2448 cands2.push_back(*it2);
2450 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2453 // Now rotate edges to bring them on Ox
2454 int startNode = cDesc2[cIDesc2[eIdx]+1];
2455 int endNode = cDesc2[cIDesc2[eIdx]+2];
2456 INTERP_KERNEL::TranslationRotationMatrix rotation;
2457 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2458 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2459 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2460 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2463 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2464 nbElemsNotM1 = tmp->getNbOfElems();
2466 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2467 double * cooPartRef(mPartRef->_coords->getPointer());
2468 double * cooPartCand(mPartCand->_coords->getPointer());
2469 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2470 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2471 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2472 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2475 // Eliminate all edges for which y or z is not null
2476 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2477 vector<int> compo; compo.push_back(1);
2478 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2480 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2481 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2482 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2483 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2484 if (!idsGoodLine->getNumberOfTuples())
2487 // Now the ordering along the Ox axis:
2488 std::vector<int> insidePoints, hitSegs;
2489 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2490 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2491 idsGoodLine->begin(), idsGoodLine->end(),
2492 /*out*/insidePoints, hitSegs);
2493 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2494 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2495 hit[cands2[*its]] = true;
2497 if (!isSplit) // current segment remains in one piece
2500 // Get original node IDs in global coords array
2501 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2502 *iit = nodeMapInv->begin()[*iit];
2504 vector<int> polyIndices, packsIds, facePack;
2505 // For each face implying this edge
2506 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2508 int faceIdx = revDescP2[ii];
2509 // each cell where this face is involved connectivity will be modified:
2510 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2512 // Current face connectivity
2513 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2514 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2516 vector<int> modifiedFace;
2517 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2518 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2519 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2523 // Rebuild 3D connectivity from descending:
2524 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2525 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2526 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2527 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2528 newConn->set3(superIdx, idx, vals);
2529 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2530 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2532 int sz, faceIdx = abs(descP[jj])-1;
2533 bool orient = descP[jj]>0;
2534 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2536 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2539 vector<int> rev(sz-1);
2540 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2541 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2545 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2548 ret = ret->buildUniqueNotSorted();