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).
314 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
316 std::size_t sz(std::distance(connBg,connEnd));
317 if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
318 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
320 INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
321 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
322 unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
323 unsigned nbOfHit(0); // number of fusions operated
324 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
325 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
326 INTERP_KERNEL::NormalizedCellType typeOfSon;
327 std::vector<int> middles;
329 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
331 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
332 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
333 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
334 posEndElt = posBaseElt+1;
336 // Look backward first: are the final edges of the cells colinear with the first ones?
337 // This initializes posBaseElt.
340 for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
342 cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
343 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
344 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
345 bool isColinear=eint->areColinears();
359 const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
360 for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell with one single edge
362 cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #j's connectivity
363 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
364 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
365 bool isColinear(eint->areColinears());
377 //push [posBaseElt,posEndElt) in newConnOfCell using e
378 // 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!
380 // at the beginning of the connectivity (insert type)
381 EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
382 else if((nbOfHit+nbOfTurn) != (nbs-1))
384 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
385 if ((nbOfHit+nbOfTurn) == (nbs-1))
386 // at the end (only quad points to deal with)
387 EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
388 posBaseElt=posEndElt;
392 newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
398 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
400 if(candidates.empty())
402 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
404 const std::vector<int>& pool(intersectEdge1[*it]);
405 int tmp[2]; tmp[0]=start; tmp[1]=stop;
406 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
411 tmp[0]=stop; tmp[1]=start;
412 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
422 * This method performs the 2nd step of Partition of 2D mesh.
423 * This method has 4 inputs :
424 * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
425 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
426 * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
427 * 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'
428 * Nodes end up lying consecutively on a cutted edge.
429 * \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.
430 * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
431 * \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.
432 * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
433 * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
435 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
436 const std::vector<double>& addCoo,
437 const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
439 int offset1=m1->getNumberOfNodes();
440 int ncell=m2->getNumberOfCells();
441 const int *c=m2->getNodalConnectivity()->begin();
442 const int *cI=m2->getNodalConnectivityIndex()->begin();
443 const double *coo=m2->getCoords()->begin();
444 const double *cooBis=m1->getCoords()->begin();
445 int offset2=offset1+m2->getNumberOfNodes();
446 intersectEdge.resize(ncell);
447 for(int i=0;i<ncell;i++,cI++)
449 const std::vector<int>& divs=subDiv[i];
450 int nnode=cI[1]-cI[0]-1;
451 std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
452 std::map<INTERP_KERNEL::Node *, int> mapp22;
453 for(int j=0;j<nnode;j++)
455 INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
456 int nnid=c[(*cI)+j+1];
457 mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
458 mapp22[nn]=nnid+offset1;
460 INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
461 for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
462 ((*it).second.first)->decrRef();
463 std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
464 std::map<INTERP_KERNEL::Node *,int> mapp3;
465 for(std::size_t j=0;j<divs.size();j++)
468 INTERP_KERNEL::Node *tmp=0;
470 tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
472 tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
474 tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
478 e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
479 for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
485 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,
486 MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
488 idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
489 idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
490 int nCells(mesh1D->getNumberOfCells());
491 if(nCells!=(int)intersectEdge2.size())
492 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
493 const DataArrayDouble *coo2(mesh1D->getCoords());
494 const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
495 const double *coo2Ptr(coo2->begin());
496 int offset1(coords1->getNumberOfTuples());
497 int offset2(offset1+coo2->getNumberOfTuples());
498 int offset3(offset2+addCoo.size()/2);
499 std::vector<double> addCooQuad;
500 MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
501 int tmp[4],cicnt(0),kk(0);
502 for(int i=0;i<nCells;i++)
504 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
505 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
506 const std::vector<int>& subEdges(intersectEdge2[i]);
507 int nbSubEdge(subEdges.size()/2);
508 for(int j=0;j<nbSubEdge;j++,kk++)
510 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));
511 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
512 INTERP_KERNEL::Edge *e2Ptr(e2);
513 std::map<int,int>::const_iterator itm;
514 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
516 tmp[0]=INTERP_KERNEL::NORM_SEG3;
517 itm=mergedNodes.find(subEdges[2*j]);
518 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
519 itm=mergedNodes.find(subEdges[2*j+1]);
520 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
521 tmp[3]=offset3+(int)addCooQuad.size()/2;
523 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
525 cOut->insertAtTheEnd(tmp,tmp+4);
526 ciOut->pushBackSilent(cicnt);
530 tmp[0]=INTERP_KERNEL::NORM_SEG2;
531 itm=mergedNodes.find(subEdges[2*j]);
532 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
533 itm=mergedNodes.find(subEdges[2*j+1]);
534 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
536 cOut->insertAtTheEnd(tmp,tmp+3);
537 ciOut->pushBackSilent(cicnt);
540 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
542 idsInRetColinear->pushBackSilent(kk);
543 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
548 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
549 ret->setConnectivity(cOut,ciOut,true);
550 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
551 arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
552 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
553 std::vector<const DataArrayDouble *> coordss(4);
554 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
555 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
560 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
562 std::vector<int> allEdges;
563 for(const int *it2(descBg);it2!=descEnd;it2++)
565 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
567 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
569 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
571 std::size_t nb(allEdges.size());
573 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
574 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
575 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
576 ret->setCoords(coords);
577 ret->allocateCells(1);
578 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
579 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
580 connOut[kk]=allEdges[2*kk];
581 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
585 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
587 const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
588 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
590 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
591 if(sz!=std::distance(descBg,descEnd))
592 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
593 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
594 std::vector<int> allEdges,centers;
595 const double *coordsPtr(coords->begin());
596 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
597 int offset(coords->getNumberOfTuples());
598 for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
600 INTERP_KERNEL::NormalizedCellType typeOfSon;
601 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
602 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
604 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
606 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
608 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
610 {//the current edge has been subsplit -> create corresponding centers.
611 std::size_t nbOfCentersToAppend(edge1.size()/2);
612 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
613 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
614 std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
615 for(std::size_t k=0;k<nbOfCentersToAppend;k++)
618 const double *aa(coordsPtr+2*(*it3++));
619 const double *bb(coordsPtr+2*(*it3++));
620 ee->getMiddleOfPoints(aa,bb,tmpp);
621 addCoo->insertAtTheEnd(tmpp,tmpp+2);
622 centers.push_back(offset+k);
626 std::size_t nb(allEdges.size());
628 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
629 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
630 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
632 ret->setCoords(coords);
635 addCoo->rearrange(2);
636 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
637 ret->setCoords(addCoo);
639 ret->allocateCells(1);
640 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
641 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
642 connOut[kk]=allEdges[2*kk];
643 connOut.insert(connOut.end(),centers.begin(),centers.end());
644 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
649 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
652 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
654 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
656 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
657 if(!cm.isQuadratic())
658 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
660 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
663 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
666 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
668 const INTERP_KERNEL::Edge *ee(*it);
669 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
673 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
676 const double *coo(mesh2D->getCoords()->begin());
677 std::size_t sz(conn.size());
678 std::vector<double> addCoo;
679 std::vector<int> conn2(conn);
680 int offset(mesh2D->getNumberOfNodes());
681 for(std::size_t i=0;i<sz;i++)
684 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
685 addCoo.insert(addCoo.end(),tmp,tmp+2);
686 conn2.push_back(offset+(int)i);
688 mesh2D->getCoords()->rearrange(1);
689 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
690 mesh2D->getCoords()->rearrange(2);
691 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
696 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
698 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
699 * a set of edges defined in \a splitMesh1D.
701 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
702 std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
704 std::size_t nb(edge1Bis.size()/2);
705 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
706 int iEnd(splitMesh1D->getNumberOfCells());
708 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
710 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
711 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
712 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
715 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
716 out0.resize(1); out1.resize(1);
717 std::vector<int>& connOut(out0[0]);
718 connOut.resize(nbOfEdgesOf2DCellSplit);
719 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
720 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
721 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
723 connOut[kk]=edge1Bis[2*kk];
724 edgesPtr[kk]=edge1BisPtr[2*kk];
729 // [i,iEnd[ contains the
730 out0.resize(2); out1.resize(2);
731 std::vector<int>& connOutLeft(out0[0]);
732 std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
733 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
734 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
735 for(std::size_t k=ii;k<jj+1;k++)
736 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
737 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
738 for(int ik=0;ik<iEnd;ik++)
740 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
741 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
744 for(int ik=iEnd-1;ik>=0;ik--)
745 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
746 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
747 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
748 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
749 for(int ik=0;ik<iEnd;ik++)
750 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
751 eright.insert(eright.end(),ees.begin(),ees.end());
759 CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
761 std::vector<int> _edges;
762 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
765 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
767 std::size_t nbe(edges.size());
768 std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
769 for(std::size_t i=0;i<nbe;i++)
771 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
772 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
774 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
775 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
776 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
782 EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
783 EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
784 bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
785 void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
786 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
790 MCAuto<MEDCouplingUMesh> _mesh;
791 MCAuto<INTERP_KERNEL::Edge> _edge;
796 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
798 const MEDCouplingUMesh *mesh(_mesh);
804 { _left++; _right++; return ; }
807 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
808 if((isLeft && isRight) || (!isLeft && !isRight))
809 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
820 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
821 if((isLeft && isRight) || (!isLeft && !isRight))
822 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
837 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
839 const MEDCouplingUMesh *mesh(_mesh);
842 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
845 {// not fully splitting cell case
846 if(mesh2D->getNumberOfCells()==1)
847 {//little optimization. 1 cell no need to find in which cell mesh is !
848 neighbors[0]=offset; neighbors[1]=offset;
853 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
854 int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
856 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
857 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
862 class VectorOfCellInfo
865 VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
866 std::size_t size() const { return _pool.size(); }
867 int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
868 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);
869 const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
870 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
871 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
872 void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
874 int getZePosOfEdgeGivenItsGlobalId(int pos) const;
875 void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
876 const CellInfo& get(int pos) const;
877 CellInfo& get(int pos);
879 std::vector<CellInfo> _pool;
880 MCAuto<MEDCouplingUMesh> _ze_mesh;
881 std::vector<EdgeInfo> _edge_info;
884 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
886 _pool[0]._edges=edges;
887 _pool[0]._edges_ptr=edgesPtr;
890 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
893 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
896 const MEDCouplingUMesh *zeMesh(_ze_mesh);
898 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
899 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
900 return zeMesh->getCellContainingPoint(barys->begin(),eps);
903 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
904 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
905 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
907 get(pos);//to check pos
908 bool isFast(pos==0 && _pool.size()==1);
909 std::size_t sz(edges.size());
910 // dealing with edges
912 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
914 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
916 std::vector<CellInfo> pool(_pool.size()-1+sz);
917 for(std::size_t i=0;i<pos;i++)
919 for(std::size_t j=0;j<sz;j++)
920 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
921 for(int i=pos+1;i<(int)_pool.size();i++)
922 pool[i+sz-1]=_pool[i];
926 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
934 std::vector< MCAuto<MEDCouplingUMesh> > ms;
937 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
941 if(pos<_ze_mesh->getNumberOfCells()-1)
943 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
946 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
947 for(std::size_t j=0;j<ms2.size();j++)
949 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
952 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
954 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
957 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
960 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
962 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
964 if((*it).isInMyRange(pos))
967 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
970 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
972 get(pos);//to perform the sanity check;
973 if(_edge_info.empty())
975 std::size_t sz(_edge_info.size()-1);
976 for(std::size_t i=0;i<sz;i++)
977 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
980 const CellInfo& VectorOfCellInfo::get(int pos) const
982 if(pos<0 || pos>=(int)_pool.size())
983 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
987 CellInfo& VectorOfCellInfo::get(int pos)
989 if(pos<0 || pos>=(int)_pool.size())
990 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
996 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
997 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
999 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1001 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1003 * \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.
1005 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1006 MCAuto<DataArrayInt>& idsLeftRight)
1008 int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1009 if(nbCellsInSplitMesh1D==0)
1010 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1011 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1012 std::size_t nb(allEdges.size()),jj;
1014 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1015 std::vector<int> edge1Bis(nb*2);
1016 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1017 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1018 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1019 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1020 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1022 idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1023 int *idsLeftRightPtr(idsLeftRight->getPointer());
1024 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1026 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1027 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1028 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1029 // of the connectivity.
1030 MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1031 renumb->alloc(nbCellsInSplitMesh1D,1);
1032 const int * renumbP(renumb->begin());
1034 int i, first=cSplitPtr[1];
1035 // Follow 1D line backward as long as it is connected:
1036 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1037 first=cSplitPtr[ciSplitPtr[i]+1];
1038 if (i < nbCellsInSplitMesh1D-1)
1040 // Build circular permutation to shift consecutive edges together
1042 renumb->applyModulus(nbCellsInSplitMesh1D);
1043 splitMesh1D->renumberCells(renumbP, false);
1044 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1045 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1051 for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1052 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1054 for(;iEnd<nbCellsInSplitMesh1D;)
1056 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1062 if(iEnd<nbCellsInSplitMesh1D)
1065 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1066 int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1068 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1069 retTmp->setCoords(splitMesh1D->getCoords());
1070 retTmp->allocateCells();
1072 std::vector< std::vector<int> > out0;
1073 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1075 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1076 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1077 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1078 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1082 for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1083 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1085 return pool.getZeMesh().retn();
1089 * splitMesh1D is an input parameter but might have its cells renumbered.
1091 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1092 const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1093 MCAuto<DataArrayInt>& idsLeftRight)
1095 const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1097 std::vector<int> allEdges;
1098 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1099 for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1101 int edgeId(std::abs(*it)-1);
1102 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1103 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1104 const std::vector<int>& edge1(intersectEdge1[edgeId]);
1106 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1108 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1109 std::size_t sz(edge1.size());
1110 for(std::size_t cnt=0;cnt<sz;cnt++)
1111 allEdgesPtr.push_back(ee);
1114 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1117 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1119 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1120 {//easy case comparison not
1121 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1123 else if(typ1.isQuadratic() && typ2.isQuadratic())
1125 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1128 if(conn1[2]==conn2[2])
1130 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1131 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1135 {//only one is quadratic
1136 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1139 const double *a(0),*bb(0),*be(0);
1140 if(typ1.isQuadratic())
1142 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1146 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1148 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1149 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1155 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1156 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1158 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1160 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1162 if(candidatesIn2DEnd==candidatesIn2DBg)
1163 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1164 const double *coo(mesh2DSplit->getCoords()->begin());
1165 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1166 return *candidatesIn2DBg;
1167 int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1168 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1169 if(cellIdInMesh1DSplitRelative<0)
1170 cur1D->changeOrientationOfCells();
1171 const int *c1D(cur1D->getNodalConnectivity()->begin());
1172 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1173 for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1175 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1176 const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1177 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1178 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1179 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1180 for(unsigned it2=0;it2<sz;it2++)
1182 INTERP_KERNEL::NormalizedCellType typeOfSon;
1183 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1184 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1185 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1189 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1193 * \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.
1194 * 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.
1195 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1196 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1197 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1198 * \param [out] addCoo - nodes to be append at the end
1199 * \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 offseted and value is id in \a m1Desc.
1201 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1202 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)
1204 static const int SPACEDIM=2;
1205 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1206 INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
1207 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1208 // Build BB tree of all edges in the tool mesh (second mesh)
1209 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1210 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1211 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1212 intersectEdge1.resize(nDescCell1);
1213 colinear2.resize(nDescCell2);
1214 subDiv2.resize(nDescCell2);
1215 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1217 std::vector<int> candidates1(1);
1218 int offset1(m1Desc->getNumberOfNodes());
1219 int offset2(offset1+m2Desc->getNumberOfNodes());
1220 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1222 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1223 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1224 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1226 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1227 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1228 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1230 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1231 // 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
1232 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1233 std::set<INTERP_KERNEL::Node *> nodes;
1234 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1235 std::size_t szz(nodes.size());
1236 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1237 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1238 for(std::size_t iii=0;iii<szz;iii++,itt++)
1239 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1240 // end of protection
1241 // Performs edge cutting:
1242 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1247 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1248 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1254 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1255 * It builds the descending connectivity of the two meshes, and then using a binary tree
1256 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1257 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1259 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1260 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1261 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1262 std::vector<double>& addCoo,
1263 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1265 // Build desc connectivity
1266 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1267 desc2=DataArrayInt::New();
1268 descIndx2=DataArrayInt::New();
1269 revDesc2=DataArrayInt::New();
1270 revDescIndx2=DataArrayInt::New();
1271 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1272 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1273 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1274 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1275 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1276 std::map<int,int> notUsedMap;
1277 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1278 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1279 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1283 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1284 * (newly created) nodes corresponding to the edge intersections.
1286 * @param[out] cr, crI connectivity of the resulting mesh
1287 * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1288 * TODO: describe input parameters
1290 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1291 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1292 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1293 const std::vector<double>& addCoords,
1294 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1296 static const int SPACEDIM=2;
1297 const double *coo1(m1->getCoords()->begin());
1298 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1299 int offset1(m1->getNumberOfNodes());
1300 const double *coo2(m2->getCoords()->begin());
1301 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1302 int offset2(offset1+m2->getNumberOfNodes());
1303 int offset3(offset2+((int)addCoords.size())/2);
1304 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1305 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1306 // Here a BBTree on 2D-cells, not on segments:
1307 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1308 int ncell1(m1->getNumberOfCells());
1310 for(int i=0;i<ncell1;i++)
1312 std::vector<int> candidates2;
1313 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1314 std::map<INTERP_KERNEL::Node *,int> mapp;
1315 std::map<int,INTERP_KERNEL::Node *> mappRev;
1316 INTERP_KERNEL::QuadraticPolygon pol1;
1317 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1318 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1319 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1320 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1321 // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
1322 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1323 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1325 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
1326 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1327 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1328 for(it1.first();!it1.finished();it1.next())
1329 edges1.insert(it1.current()->getPtr());
1331 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1332 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1334 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1336 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1337 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1338 // Complete mapping with elements coming from the current cell it2 in mesh2:
1339 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1340 // pol2 is the new QP in the final merged result.
1341 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1342 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1345 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1347 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1348 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1349 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1350 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1352 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1353 // by m2 but that we still want to keep in the final result.
1358 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1360 catch(INTERP_KERNEL::Exception& e)
1362 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();
1363 throw INTERP_KERNEL::Exception(oss.str());
1366 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1367 (*it).second->decrRef();
1371 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1373 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1376 std::size_t sz(conn.size());
1377 std::size_t found(std::numeric_limits<std::size_t>::max());
1378 for(std::size_t i=0;i<sz;i++)
1380 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1381 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]};
1382 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1383 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1384 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1386 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];
1387 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]);
1389 if(dotTest>eps && dotTest<1.-eps)
1395 if(found==std::numeric_limits<std::size_t>::max())
1396 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1397 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1400 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1402 std::size_t sz(conn.size());
1403 std::vector<int> *curPart(&part0);
1404 for(std::size_t i=0;i<sz;i++)
1406 int nextt(conn[(i+1)%sz]);
1407 (*curPart).push_back(nextt);
1408 if(nextt==pt0 || nextt==pt1)
1414 (*curPart).push_back(nextt);
1420 * 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.
1422 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1423 const int *desc, const int *descIndx, const double *coords, double eps,
1424 std::vector<std::vector<int> >& res) const
1426 checkFullyDefined();
1427 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1428 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1429 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1430 int nbOfCells(getNumberOfCells());
1432 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1433 for(int i=0;i<nbOfCells;i++)
1435 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1436 for(int j=0;j<nbOfFaces;j++)
1438 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1439 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1440 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1441 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1442 INTERP_KERNEL::NormalizedCellType cmsId;
1443 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1444 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1445 if(p.first!=-1 && p.second!=-1)
1449 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1450 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1451 std::vector<int> elt1,elt2;
1452 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1453 res.push_back(elt1);
1454 res.push_back(elt2);
1466 * 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).
1468 * \sa MEDCouplingUMesh::split2DCells
1470 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1472 checkConnectivityFullyDefined();
1473 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1474 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1475 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1476 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1477 int prevPosOfCi(ciPtr[0]);
1478 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1480 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1481 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1482 for(int j=0;j<sz;j++)
1484 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1485 for(int k=0;k<sz2;k++)
1486 *cPtr++=subPtr[offset2+k];
1488 *cPtr++=oldConn[prevPosOfCi+j+2];
1491 prevPosOfCi=ciPtr[1];
1492 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1495 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1496 _nodal_connec->decrRef();
1497 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1502 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1504 * \return int - the number of new nodes created.
1505 * \sa MEDCouplingUMesh::split2DCells
1507 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1509 checkConsistencyLight();
1510 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1511 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1512 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1513 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1514 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1515 const double *oldCoordsPtr(getCoords()->begin());
1516 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1517 int prevPosOfCi(ciPtr[0]);
1518 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1520 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1521 for(int j=0;j<sz;j++)
1522 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1523 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1524 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1526 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1530 cPtr[1]=oldConn[prevPosOfCi+2+j];
1531 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1534 std::vector<INTERP_KERNEL::Node *> ns(3);
1535 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1536 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1537 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1538 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1539 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1541 cPtr[1]=subPtr[offset2+k];
1542 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1544 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1547 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1549 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1550 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1553 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1554 _nodal_connec->decrRef();
1555 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1556 addCoo->rearrange(2);
1557 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1559 return addCoo->getNumberOfTuples();
1566 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1567 * returns a result mesh constituted by polygons.
1568 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1569 * all nodes from m2.
1570 * The meshes should be in 2D space. In
1571 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1573 * \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
1574 * 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)
1575 * \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
1576 * 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)
1577 * \param [in] eps - precision used to detect coincident mesh entities.
1578 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1579 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1580 * this array using decrRef() as it is no more needed.
1581 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1582 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1583 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1584 * any cell of \a m2. The caller is to delete this array using decrRef() as
1585 * it is no more needed.
1586 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1587 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1588 * is no more needed.
1589 * \throw If the coordinates array is not set in any of the meshes.
1590 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1591 * \throw If any of the meshes is not a 2D mesh in 2D space.
1593 * \sa conformize2D, mergeNodes
1595 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1596 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1599 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1600 m1->checkFullyDefined();
1601 m2->checkFullyDefined();
1602 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1603 INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
1604 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1605 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1607 // Step 1: compute all edge intersections (new nodes)
1608 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1609 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1610 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1611 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1612 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1613 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1614 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1615 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1616 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1617 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1619 // Step 2: re-order newly created nodes according to the ordering found in m2
1620 std::vector< std::vector<int> > intersectEdge2;
1621 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1622 subDiv2.clear(); dd5=0; dd6=0;
1625 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1626 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1627 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1628 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1630 // Step 4: Prepare final result:
1631 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1632 addCooDa->alloc((int)(addCoo.size())/2,2);
1633 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1634 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1635 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1636 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1637 std::vector<const DataArrayDouble *> coordss(4);
1638 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1639 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1640 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1641 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1642 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1643 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1644 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1645 ret->setConnectivity(conn,connI,true);
1646 ret->setCoords(coo);
1647 cellNb1=c1.retn(); cellNb2=c2.retn();
1652 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1653 * 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
1654 * and finally, in case of quadratic polygon the centers of edges new nodes.
1655 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1657 * \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
1658 * 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)
1659 * \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
1660 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1661 * \param [in] eps - precision used to perform intersections and localization operations.
1662 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1663 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1664 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1665 * 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.
1666 * \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
1667 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1668 * 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.
1670 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1672 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1674 if(!mesh2D || !mesh1D)
1675 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1676 mesh2D->checkFullyDefined();
1677 mesh1D->checkFullyDefined();
1678 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1679 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1680 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1681 // Step 1: compute all edge intersections (new nodes)
1682 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1683 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1684 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1685 INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
1687 // Build desc connectivity
1688 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1689 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1690 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1691 std::map<int,int> mergedNodes;
1692 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1693 // use mergeNodes to fix intersectEdge1
1694 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1696 std::size_t n((*it0).size()/2);
1697 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1698 std::map<int,int>::const_iterator it1;
1699 it1=mergedNodes.find(eltStart);
1700 if(it1!=mergedNodes.end())
1701 (*it0)[0]=(*it1).second;
1702 it1=mergedNodes.find(eltEnd);
1703 if(it1!=mergedNodes.end())
1704 (*it0)[2*n-1]=(*it1).second;
1707 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1708 addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
1709 // Step 2: re-order newly created nodes according to the ordering found in m2
1710 std::vector< std::vector<int> > intersectEdge2;
1711 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1713 // Step 3: compute splitMesh1D
1714 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1715 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1716 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1717 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1718 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1719 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1720 // deal with cells in mesh2D that are not cut but only some of their edges are
1721 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1722 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1723 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1724 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
1725 if(!idsInDesc2DToBeRefined->empty())
1727 DataArrayInt *out0(0),*outi0(0);
1728 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1729 MCAuto<DataArrayInt> outi0s(outi0);
1731 out0s=out0s->buildUnique();
1735 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1736 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1737 MCAuto<DataArrayInt> elts,eltsIndex;
1738 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1739 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1740 if (eltsIndex->getNumberOfTuples() > 1)
1741 eltsIndex2 = eltsIndex->deltaShiftIndex();
1742 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1743 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1744 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1745 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1746 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1747 if((DataArrayInt *)out0s)
1748 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1749 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1750 // OK all is ready to insert in ret2 mesh
1751 if(!untouchedCells->empty())
1752 {// the most easy part, cells in mesh2D not impacted at all
1753 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1754 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1755 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1757 if((DataArrayInt *)out0s)
1758 {// here dealing with cells in out0s but not in cellsToBeModified
1759 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1760 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1761 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1763 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1764 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1766 int offset(ret2->getNumberOfTuples());
1767 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1768 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1769 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1770 int kk(0),*ret3ptr(partOfRet3->getPointer());
1771 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1773 int faceId(std::abs(*it)-1);
1774 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1776 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1779 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1780 ret3ptr[2*kk]=tmp+offset;
1781 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1782 ret3ptr[2*kk+1]=tmp+offset;
1785 {//the current edge is shared by a 2D cell that will be split just after
1786 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1787 ret3ptr[2*kk]=-(*it2+1);
1788 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1789 ret3ptr[2*kk+1]=-(*it2+1);
1793 m1Desc->setCoords(ret1->getCoords());
1794 ret1NonCol->setCoords(ret1->getCoords());
1795 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1796 if(!outMesh2DSplit.empty())
1798 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1799 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1800 (*itt)->setCoords(da);
1803 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1804 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1806 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1807 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1808 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1809 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1810 MCAuto<DataArrayInt> partOfRet3;
1811 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));
1812 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1813 outMesh2DSplit.push_back(splitOfOneCell);
1814 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1815 ret2->pushBackSilent(*it);
1818 std::size_t nbOfMeshes(outMesh2DSplit.size());
1819 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1820 for(std::size_t i=0;i<nbOfMeshes;i++)
1821 tmp[i]=outMesh2DSplit[i];
1823 ret1->getCoords()->setInfoOnComponents(compNames);
1824 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1825 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1827 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1828 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1830 int old2DCellId(-ret3->getIJ(*it,0)-1);
1831 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1832 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
1834 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1837 splitMesh1D=ret1.retn();
1838 splitMesh2D=ret2D.retn();
1839 cellIdInMesh2D=ret2.retn();
1840 cellIdInMesh1D=ret3.retn();
1844 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1845 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1846 * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this this method
1847 * 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).
1848 * 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.
1850 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1851 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1853 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1854 * This method expects that all nodes in \a this are not closer than \a eps.
1855 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1857 * \param [in] eps the relative error to detect merged edges.
1858 * \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
1859 * that the user is expected to deal with.
1861 * \throw If \a this is not coherent.
1862 * \throw If \a this has not spaceDim equal to 2.
1863 * \throw If \a this has not meshDim equal to 2.
1864 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1866 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1868 static const int SPACEDIM=2;
1869 checkConsistencyLight();
1870 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1871 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1872 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1873 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1874 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1875 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1876 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1877 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1878 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1879 std::vector<double> addCoo;
1880 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1881 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1882 INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
1883 for(int i=0;i<nDescCell;i++)
1885 std::vector<int> candidates;
1886 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1887 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1890 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1891 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1892 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1893 INTERP_KERNEL::MergePoints merge;
1894 INTERP_KERNEL::QuadraticPolygon c1,c2;
1895 e1->intersectWith(e2,merge,c1,c2);
1896 e1->decrRef(); e2->decrRef();
1897 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1898 overlapEdge[i].push_back(*it);
1899 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1900 overlapEdge[*it].push_back(i);
1903 // splitting done. sort intersect point in intersectEdge.
1904 std::vector< std::vector<int> > middle(nDescCell);
1905 int nbOf2DCellsToBeSplit(0);
1906 bool middleNeedsToBeUsed(false);
1907 std::vector<bool> cells2DToTreat(nDescCell,false);
1908 for(int i=0;i<nDescCell;i++)
1910 std::vector<int>& isect(intersectEdge[i]);
1911 int sz((int)isect.size());
1914 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1915 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1916 e->sortSubNodesAbs(coords,isect);
1921 int idx0(rdi[i]),idx1(rdi[i+1]);
1923 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1924 if(!cells2DToTreat[rd[idx0]])
1926 cells2DToTreat[rd[idx0]]=true;
1927 nbOf2DCellsToBeSplit++;
1929 // try to reuse at most eventual 'middle' of SEG3
1930 std::vector<int>& mid(middle[i]);
1931 mid.resize(sz+1,-1);
1932 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1934 middleNeedsToBeUsed=true;
1935 const std::vector<int>& candidates(overlapEdge[i]);
1936 std::vector<int> trueCandidates;
1937 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1938 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1939 trueCandidates.push_back(*itc);
1940 int stNode(c[ci[i]+1]),endNode(isect[0]);
1941 for(int j=0;j<sz+1;j++)
1943 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1945 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1946 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1947 { mid[j]=*itc; break; }
1950 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1955 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1956 if(nbOf2DCellsToBeSplit==0)
1959 int *retPtr(ret->getPointer());
1960 for(int i=0;i<nCell;i++)
1961 if(cells2DToTreat[i])
1964 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1965 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1966 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1967 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1968 if(middleNeedsToBeUsed)
1969 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1970 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1971 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1972 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.
1973 setPartOfMySelf(ret->begin(),ret->end(),*modif);
1975 bool areNodesMerged; int newNbOfNodes;
1976 if(nbOfNodesCreated!=0)
1977 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
1983 * 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.
1984 * 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).
1985 * 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
1986 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
1987 * 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
1988 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
1990 * 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
1991 * using new instance, idem for coordinates.
1993 * 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.
1995 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
1997 * \throw If \a this is not coherent.
1998 * \throw If \a this has not spaceDim equal to 2.
1999 * \throw If \a this has not meshDim equal to 2.
2001 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2003 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2005 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2006 checkConsistencyLight();
2007 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2008 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2009 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2010 INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
2011 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2012 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2013 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2014 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2015 const double *coords(_coords->begin());
2016 int *newciptr(newci->getPointer());
2017 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2019 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
2020 ret->pushBackSilent(i);
2021 newciptr[1]=newc->getNumberOfTuples();
2026 if(!appendedCoords->empty())
2028 appendedCoords->rearrange(2);
2029 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2031 setCoords(newCoords);
2034 setConnectivity(newc,newci,true);
2040 * c, cI describe a wire mesh in 3D space, in local numbering
2041 * startNode, endNode in global numbering
2042 *\return true if the segment is indeed split
2044 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2045 const int * c, const int * cI, const int *idsBg, const int *endBg,
2046 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2048 using namespace std;
2050 const int SPACEDIM=3;
2051 typedef pair<double, int> PairDI;
2053 for (const int * it = idsBg; it != endBg; ++it)
2055 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2056 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2057 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2058 x.insert(make_pair(coo[end*SPACEDIM], end));
2061 vector<PairDI> xx(x.begin(), x.end());
2062 sort(xx.begin(),xx.end());
2063 pointIds.reserve(xx.size());
2065 // Keep what is inside [startNode, endNode]:
2067 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2069 const int idx = (*it).second;
2072 if (idx == startNode) go = 1;
2073 if (idx == endNode) go = 2;
2074 if (go) pointIds.push_back(idx);
2077 pointIds.push_back(idx);
2078 if (idx == endNode || idx == startNode)
2082 // vector<int> pointIds2(pointIds.size()+2);
2083 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2084 // pointIds2[0] = startNode;
2085 // pointIds2[pointIds2.size()-1] = endNode;
2088 reverse(pointIds.begin(), pointIds.end());
2090 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2091 for (const int * it = idsBg; it != endBg; ++it)
2093 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2094 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2095 if (itStart == pointIds.end()) continue;
2096 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2097 if (itEnd == pointIds.end()) continue;
2098 if (abs(distance(itEnd, itStart)) != 1) continue;
2099 hitSegs.push_back(*it); // segment is undivided.
2102 return (pointIds.size() > 2); // something else apart start and end node
2105 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2106 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2108 using namespace std;
2109 int dst = distance(sIdxConn, sIdxConnE);
2110 modifiedFace.reserve(dst + insidePoints.size()-2);
2111 modifiedFace.resize(dst);
2112 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2114 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2115 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2116 if (startPos == shortEnd)
2117 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2118 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2119 if (endPos == shortEnd)
2120 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2121 int d = distance(startPos, endPos);
2122 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2123 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2125 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2132 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2133 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2134 * This method performs a conformization of \b this.
2136 * Only polyhedron cells are supported. You can call convertAllToPoly()
2138 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2139 * This method expects that all nodes in \a this are not closer than \a eps.
2140 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2142 * \param [in] eps the relative error to detect merged edges.
2143 * \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
2144 * that the user is expected to deal with.
2146 * \throw If \a this is not coherent.
2147 * \throw If \a this has not spaceDim equal to 3.
2148 * \throw If \a this has not meshDim equal to 3.
2149 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2151 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2153 using namespace std;
2155 static const int SPACEDIM=3;
2156 checkConsistencyLight();
2157 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2158 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2159 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2160 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2162 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2163 const double * coo(_coords->begin());
2164 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2167 /*************************
2169 *************************/
2170 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2171 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2172 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2173 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2174 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2177 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2178 const double *bbox(bboxArr->begin()); getCoords()->begin();
2179 int nDescCell(mDesc->getNumberOfCells());
2180 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2181 // Surfaces - handle biggest first
2182 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2183 DataArrayDouble * surfs = surfF->getArray();
2185 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2186 DataArrayDouble * normals = normalsF->getArray();
2187 const double * normalsP = normals->getConstPointer();
2189 // Sort faces by decreasing surface:
2190 vector< pair<double,int> > S;
2191 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2193 pair<double,int> p = make_pair(surfs->begin()[i], i);
2196 sort(S.rbegin(),S.rend()); // reverse sort
2197 vector<bool> hit(nDescCell);
2198 fill(hit.begin(), hit.end(), false);
2199 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2201 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2203 int faceIdx = (*it).second;
2204 if (hit[faceIdx]) continue;
2206 vector<int> candidates, cands2;
2207 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2208 // Keep only candidates whose normal matches the normal of current face
2209 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2211 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2212 if (*it2 != faceIdx && col)
2213 cands2.push_back(*it2);
2218 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2219 INTERP_KERNEL::TranslationRotationMatrix rotation;
2220 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2221 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2222 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2224 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2225 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2226 double * cooPartRef(mPartRef->_coords->getPointer());
2227 double * cooPartCand(mPartCand->_coords->getPointer());
2228 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2229 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2230 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2231 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2233 // Localize faces in 2D thanks to barycenters
2234 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2235 vector<int> compo; compo.push_back(2);
2236 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2237 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2238 if (!idsGoodPlane->getNumberOfTuples())
2241 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2243 compo[0] = 0; compo.push_back(1);
2244 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2245 mPartRef->changeSpaceDimension(2,0.0);
2246 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2247 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2249 if (!cc->getNumberOfTuples())
2251 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2254 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2255 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2258 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2259 throw INTERP_KERNEL::Exception(oss.str());
2263 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2265 if (!ids->getNumberOfTuples())
2268 double checkSurf=0.0;
2269 const int * idsGoodPlaneP(idsGoodPlane->begin());
2270 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2272 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2273 hit[faceIdx2] = true;
2274 checkSurf += surfs->begin()[faceIdx2];
2276 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2279 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2280 throw INTERP_KERNEL::Exception(oss.str());
2283 // For all polyhedrons using this face, replace connectivity:
2284 vector<int> polyIndices, packsIds, facePack;
2285 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2286 polyIndices.push_back(revDescP[ii]);
2287 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2289 // Current face connectivity
2290 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2291 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2292 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2293 // Deletion of old faces
2295 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2297 if (packsIds[jj] == -1)
2298 // The below should never happen - if a face is used several times, with a different layout of the nodes
2299 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2300 // faces which are actually used only once, by a single cell. This is different for edges below.
2301 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2303 connSla->deletePack(*it2, packsIds[jj]);
2305 // Insertion of new faces:
2306 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2308 // Build pack from the face to insert:
2309 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2311 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2312 // Insert it in all hit polyhedrons:
2313 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2314 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2319 // Set back modified connectivity
2320 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2321 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2322 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2325 /************************
2327 ************************/
2328 // Now we have a face-conform mesh.
2330 // Recompute descending
2331 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2332 // Rebuild desc connectivity with orientation this time!!
2333 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2334 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2335 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2336 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2337 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2338 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2339 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2340 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2341 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2342 // std::cout << "writing!\n";
2343 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2344 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2345 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2346 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2347 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2348 const double *bbox2(bboxArr->begin());
2349 int nDesc2Cell=mDesc2->getNumberOfCells();
2350 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2352 // Edges - handle longest first
2353 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2354 DataArrayDouble * lens = lenF->getArray();
2356 // Sort edges by decreasing length:
2357 vector<pair<double,int> > S;
2358 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2360 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2363 sort(S.rbegin(),S.rend()); // reverse sort
2365 vector<bool> hit(nDesc2Cell);
2366 fill(hit.begin(), hit.end(), false);
2368 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2370 int eIdx = (*it).second;
2374 vector<int> candidates, cands2;
2375 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2376 // Keep only candidates colinear with current edge
2378 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2379 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2380 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2381 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2384 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2385 for (int i3=0; i3 < 3; i3++)
2386 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2387 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2388 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2389 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2391 cands2.push_back(*it2);
2393 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2396 // Now rotate edges to bring them on Ox
2397 int startNode = cDesc2[cIDesc2[eIdx]+1];
2398 int endNode = cDesc2[cIDesc2[eIdx]+2];
2399 INTERP_KERNEL::TranslationRotationMatrix rotation;
2400 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2401 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2402 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2403 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2406 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2407 nbElemsNotM1 = tmp->getNbOfElems();
2409 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2410 double * cooPartRef(mPartRef->_coords->getPointer());
2411 double * cooPartCand(mPartCand->_coords->getPointer());
2412 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2413 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2414 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2415 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2418 // Eliminate all edges for which y or z is not null
2419 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2420 vector<int> compo; compo.push_back(1);
2421 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2423 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2424 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2425 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2426 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2427 if (!idsGoodLine->getNumberOfTuples())
2430 // Now the ordering along the Ox axis:
2431 std::vector<int> insidePoints, hitSegs;
2432 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2433 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2434 idsGoodLine->begin(), idsGoodLine->end(),
2435 /*out*/insidePoints, hitSegs);
2436 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2437 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2438 hit[cands2[*its]] = true;
2440 if (!isSplit) // current segment remains in one piece
2443 // Get original node IDs in global coords array
2444 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2445 *iit = nodeMapInv->begin()[*iit];
2447 vector<int> polyIndices, packsIds, facePack;
2448 // For each face implying this edge
2449 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2451 int faceIdx = revDescP2[ii];
2452 // each cell where this face is involved connectivity will be modified:
2453 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2455 // Current face connectivity
2456 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2457 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2459 vector<int> modifiedFace;
2460 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2461 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2462 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2466 // Rebuild 3D connectivity from descending:
2467 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2468 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2469 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2470 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2471 newConn->set3(superIdx, idx, vals);
2472 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2473 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2475 int sz, faceIdx = abs(descP[jj])-1;
2476 bool orient = descP[jj]>0;
2477 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2479 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2482 vector<int> rev(sz-1);
2483 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2484 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2488 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2491 ret = ret->buildUniqueNotSorted();