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 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1207 // Build BB tree of all edges in the tool mesh (second mesh)
1208 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1209 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1210 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1211 intersectEdge1.resize(nDescCell1);
1212 colinear2.resize(nDescCell2);
1213 subDiv2.resize(nDescCell2);
1214 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1216 std::vector<int> candidates1(1);
1217 int offset1(m1Desc->getNumberOfNodes());
1218 int offset2(offset1+m2Desc->getNumberOfNodes());
1219 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1221 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1222 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1223 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1225 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1226 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1227 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1229 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1230 // 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
1231 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1232 std::set<INTERP_KERNEL::Node *> nodes;
1233 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1234 std::size_t szz(nodes.size());
1235 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1236 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1237 for(std::size_t iii=0;iii<szz;iii++,itt++)
1238 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1239 // end of protection
1240 // Performs edge cutting:
1241 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1246 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1247 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1253 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1254 * It builds the descending connectivity of the two meshes, and then using a binary tree
1255 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1256 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1258 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1259 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1260 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1261 std::vector<double>& addCoo,
1262 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1264 // Build desc connectivity
1265 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1266 desc2=DataArrayInt::New();
1267 descIndx2=DataArrayInt::New();
1268 revDesc2=DataArrayInt::New();
1269 revDescIndx2=DataArrayInt::New();
1270 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1271 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1272 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1273 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1274 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1275 std::map<int,int> notUsedMap;
1276 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1277 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1278 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1282 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1283 * (newly created) nodes corresponding to the edge intersections.
1285 * @param[out] cr, crI connectivity of the resulting mesh
1286 * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1287 * TODO: describe input parameters
1289 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1290 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1291 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1292 const std::vector<double>& addCoords,
1293 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1295 static const int SPACEDIM=2;
1296 const double *coo1(m1->getCoords()->begin());
1297 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1298 int offset1(m1->getNumberOfNodes());
1299 const double *coo2(m2->getCoords()->begin());
1300 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1301 int offset2(offset1+m2->getNumberOfNodes());
1302 int offset3(offset2+((int)addCoords.size())/2);
1303 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1304 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1305 // Here a BBTree on 2D-cells, not on segments:
1306 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1307 int ncell1(m1->getNumberOfCells());
1309 for(int i=0;i<ncell1;i++)
1311 std::vector<int> candidates2;
1312 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1313 std::map<INTERP_KERNEL::Node *,int> mapp;
1314 std::map<int,INTERP_KERNEL::Node *> mappRev;
1315 INTERP_KERNEL::QuadraticPolygon pol1;
1316 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1317 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1318 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1319 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1320 // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
1321 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1322 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1324 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
1325 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1326 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1327 for(it1.first();!it1.finished();it1.next())
1328 edges1.insert(it1.current()->getPtr());
1330 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1331 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1333 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1335 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1336 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1337 // Complete mapping with elements coming from the current cell it2 in mesh2:
1338 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1339 // pol2 is the new QP in the final merged result.
1340 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1341 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1344 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1346 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1347 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1348 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1349 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1351 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1352 // by m2 but that we still want to keep in the final result.
1357 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1359 catch(INTERP_KERNEL::Exception& e)
1361 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();
1362 throw INTERP_KERNEL::Exception(oss.str());
1365 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1366 (*it).second->decrRef();
1370 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1372 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1375 std::size_t sz(conn.size());
1376 std::size_t found(std::numeric_limits<std::size_t>::max());
1377 for(std::size_t i=0;i<sz;i++)
1379 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1380 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]};
1381 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1382 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1383 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1385 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];
1386 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]);
1388 if(dotTest>eps && dotTest<1.-eps)
1394 if(found==std::numeric_limits<std::size_t>::max())
1395 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1396 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1399 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1401 std::size_t sz(conn.size());
1402 std::vector<int> *curPart(&part0);
1403 for(std::size_t i=0;i<sz;i++)
1405 int nextt(conn[(i+1)%sz]);
1406 (*curPart).push_back(nextt);
1407 if(nextt==pt0 || nextt==pt1)
1413 (*curPart).push_back(nextt);
1419 * 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.
1421 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1422 const int *desc, const int *descIndx, const double *coords, double eps,
1423 std::vector<std::vector<int> >& res) const
1425 checkFullyDefined();
1426 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1427 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1428 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1429 int nbOfCells(getNumberOfCells());
1431 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1432 for(int i=0;i<nbOfCells;i++)
1434 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1435 for(int j=0;j<nbOfFaces;j++)
1437 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1438 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1439 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1440 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1441 INTERP_KERNEL::NormalizedCellType cmsId;
1442 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1443 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1444 if(p.first!=-1 && p.second!=-1)
1448 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1449 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1450 std::vector<int> elt1,elt2;
1451 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1452 res.push_back(elt1);
1453 res.push_back(elt2);
1465 * 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).
1467 * \sa MEDCouplingUMesh::split2DCells
1469 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1471 checkConnectivityFullyDefined();
1472 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1473 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1474 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1475 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1476 int prevPosOfCi(ciPtr[0]);
1477 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1479 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1480 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1481 for(int j=0;j<sz;j++)
1483 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1484 for(int k=0;k<sz2;k++)
1485 *cPtr++=subPtr[offset2+k];
1487 *cPtr++=oldConn[prevPosOfCi+j+2];
1490 prevPosOfCi=ciPtr[1];
1491 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1494 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1495 _nodal_connec->decrRef();
1496 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1501 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1503 * \return int - the number of new nodes created.
1504 * \sa MEDCouplingUMesh::split2DCells
1506 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1508 checkConsistencyLight();
1509 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1510 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1511 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1512 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1513 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1514 const double *oldCoordsPtr(getCoords()->begin());
1515 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1516 int prevPosOfCi(ciPtr[0]);
1517 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1519 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1520 for(int j=0;j<sz;j++)
1521 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1522 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1523 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1525 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1529 cPtr[1]=oldConn[prevPosOfCi+2+j];
1530 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1533 std::vector<INTERP_KERNEL::Node *> ns(3);
1534 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1535 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1536 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1537 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1538 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1540 cPtr[1]=subPtr[offset2+k];
1541 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1543 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1546 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1548 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1549 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1552 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1553 _nodal_connec->decrRef();
1554 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1555 addCoo->rearrange(2);
1556 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1558 return addCoo->getNumberOfTuples();
1565 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1566 * returns a result mesh constituted by polygons.
1567 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1568 * all nodes from m2.
1569 * The meshes should be in 2D space. In
1570 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1572 * \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
1573 * 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)
1574 * \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
1575 * 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)
1576 * \param [in] eps - precision used to detect coincident mesh entities.
1577 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1578 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1579 * this array using decrRef() as it is no more needed.
1580 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1581 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1582 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1583 * any cell of \a m2. The caller is to delete this array using decrRef() as
1584 * it is no more needed.
1585 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1586 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1587 * is no more needed.
1588 * \throw If the coordinates array is not set in any of the meshes.
1589 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1590 * \throw If any of the meshes is not a 2D mesh in 2D space.
1592 * \sa conformize2D, mergeNodes
1594 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1595 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1598 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1599 m1->checkFullyDefined();
1600 m2->checkFullyDefined();
1601 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1602 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1603 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1605 // Step 1: compute all edge intersections (new nodes)
1606 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1607 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1608 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1609 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1610 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1611 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1612 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1613 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1614 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1615 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1617 // Step 2: re-order newly created nodes according to the ordering found in m2
1618 std::vector< std::vector<int> > intersectEdge2;
1619 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1620 subDiv2.clear(); dd5=0; dd6=0;
1623 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1624 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1625 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1626 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1628 // Step 4: Prepare final result:
1629 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1630 addCooDa->alloc((int)(addCoo.size())/2,2);
1631 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1632 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1633 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1634 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1635 std::vector<const DataArrayDouble *> coordss(4);
1636 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1637 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1638 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1639 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1640 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1641 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1642 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1643 ret->setConnectivity(conn,connI,true);
1644 ret->setCoords(coo);
1645 cellNb1=c1.retn(); cellNb2=c2.retn();
1650 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1651 * 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
1652 * and finally, in case of quadratic polygon the centers of edges new nodes.
1653 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1655 * \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
1656 * 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)
1657 * \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
1658 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1659 * \param [in] eps - precision used to perform intersections and localization operations.
1660 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1661 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1662 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1663 * 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.
1664 * \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
1665 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1666 * 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.
1668 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1670 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1672 if(!mesh2D || !mesh1D)
1673 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1674 mesh2D->checkFullyDefined();
1675 mesh1D->checkFullyDefined();
1676 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1677 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1678 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1679 // Step 1: compute all edge intersections (new nodes)
1680 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1681 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1682 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1684 // Build desc connectivity
1685 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1686 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1687 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1688 std::map<int,int> mergedNodes;
1689 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1690 // use mergeNodes to fix intersectEdge1
1691 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1693 std::size_t n((*it0).size()/2);
1694 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1695 std::map<int,int>::const_iterator it1;
1696 it1=mergedNodes.find(eltStart);
1697 if(it1!=mergedNodes.end())
1698 (*it0)[0]=(*it1).second;
1699 it1=mergedNodes.find(eltEnd);
1700 if(it1!=mergedNodes.end())
1701 (*it0)[2*n-1]=(*it1).second;
1704 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1705 addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
1706 // Step 2: re-order newly created nodes according to the ordering found in m2
1707 std::vector< std::vector<int> > intersectEdge2;
1708 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1710 // Step 3: compute splitMesh1D
1711 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1712 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1713 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1714 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1715 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1716 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1717 // deal with cells in mesh2D that are not cut but only some of their edges are
1718 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1719 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1720 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1721 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
1722 if(!idsInDesc2DToBeRefined->empty())
1724 DataArrayInt *out0(0),*outi0(0);
1725 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1726 MCAuto<DataArrayInt> outi0s(outi0);
1728 out0s=out0s->buildUnique();
1732 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1733 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1734 MCAuto<DataArrayInt> elts,eltsIndex;
1735 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1736 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1737 if (eltsIndex->getNumberOfTuples() > 1)
1738 eltsIndex2 = eltsIndex->deltaShiftIndex();
1739 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1740 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1741 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1742 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1743 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1744 if((DataArrayInt *)out0s)
1745 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1746 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1747 // OK all is ready to insert in ret2 mesh
1748 if(!untouchedCells->empty())
1749 {// the most easy part, cells in mesh2D not impacted at all
1750 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1751 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1752 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1754 if((DataArrayInt *)out0s)
1755 {// here dealing with cells in out0s but not in cellsToBeModified
1756 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1757 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1758 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1760 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1761 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1763 int offset(ret2->getNumberOfTuples());
1764 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1765 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1766 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1767 int kk(0),*ret3ptr(partOfRet3->getPointer());
1768 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1770 int faceId(std::abs(*it)-1);
1771 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1773 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1776 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1777 ret3ptr[2*kk]=tmp+offset;
1778 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1779 ret3ptr[2*kk+1]=tmp+offset;
1782 {//the current edge is shared by a 2D cell that will be split just after
1783 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1784 ret3ptr[2*kk]=-(*it2+1);
1785 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1786 ret3ptr[2*kk+1]=-(*it2+1);
1790 m1Desc->setCoords(ret1->getCoords());
1791 ret1NonCol->setCoords(ret1->getCoords());
1792 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1793 if(!outMesh2DSplit.empty())
1795 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1796 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1797 (*itt)->setCoords(da);
1800 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1801 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1803 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1804 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1805 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1806 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1807 MCAuto<DataArrayInt> partOfRet3;
1808 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));
1809 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1810 outMesh2DSplit.push_back(splitOfOneCell);
1811 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1812 ret2->pushBackSilent(*it);
1815 std::size_t nbOfMeshes(outMesh2DSplit.size());
1816 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1817 for(std::size_t i=0;i<nbOfMeshes;i++)
1818 tmp[i]=outMesh2DSplit[i];
1820 ret1->getCoords()->setInfoOnComponents(compNames);
1821 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1822 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1824 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1825 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1827 int old2DCellId(-ret3->getIJ(*it,0)-1);
1828 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1829 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
1831 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1834 splitMesh1D=ret1.retn();
1835 splitMesh2D=ret2D.retn();
1836 cellIdInMesh2D=ret2.retn();
1837 cellIdInMesh1D=ret3.retn();
1841 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1842 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1843 * 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
1844 * 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).
1845 * 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.
1847 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1848 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1850 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1851 * This method expects that all nodes in \a this are not closer than \a eps.
1852 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1854 * \param [in] eps the relative error to detect merged edges.
1855 * \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
1856 * that the user is expected to deal with.
1858 * \throw If \a this is not coherent.
1859 * \throw If \a this has not spaceDim equal to 2.
1860 * \throw If \a this has not meshDim equal to 2.
1861 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1863 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1865 static const int SPACEDIM=2;
1866 checkConsistencyLight();
1867 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1868 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1869 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1870 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1871 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1872 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1873 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1874 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1875 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1876 std::vector<double> addCoo;
1877 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1878 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1879 for(int i=0;i<nDescCell;i++)
1881 std::vector<int> candidates;
1882 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1883 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1884 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1886 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1887 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1888 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1889 INTERP_KERNEL::MergePoints merge;
1890 INTERP_KERNEL::QuadraticPolygon c1,c2;
1891 e1->intersectWith(e2,merge,c1,c2);
1892 e1->decrRef(); e2->decrRef();
1893 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1894 overlapEdge[i].push_back(*it);
1895 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1896 overlapEdge[*it].push_back(i);
1899 // splitting done. sort intersect point in intersectEdge.
1900 std::vector< std::vector<int> > middle(nDescCell);
1901 int nbOf2DCellsToBeSplit(0);
1902 bool middleNeedsToBeUsed(false);
1903 std::vector<bool> cells2DToTreat(nDescCell,false);
1904 for(int i=0;i<nDescCell;i++)
1906 std::vector<int>& isect(intersectEdge[i]);
1907 int sz((int)isect.size());
1910 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1911 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1912 e->sortSubNodesAbs(coords,isect);
1917 int idx0(rdi[i]),idx1(rdi[i+1]);
1919 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1920 if(!cells2DToTreat[rd[idx0]])
1922 cells2DToTreat[rd[idx0]]=true;
1923 nbOf2DCellsToBeSplit++;
1925 // try to reuse at most eventual 'middle' of SEG3
1926 std::vector<int>& mid(middle[i]);
1927 mid.resize(sz+1,-1);
1928 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1930 middleNeedsToBeUsed=true;
1931 const std::vector<int>& candidates(overlapEdge[i]);
1932 std::vector<int> trueCandidates;
1933 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1934 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1935 trueCandidates.push_back(*itc);
1936 int stNode(c[ci[i]+1]),endNode(isect[0]);
1937 for(int j=0;j<sz+1;j++)
1939 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1941 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1942 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1943 { mid[j]=*itc; break; }
1946 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1951 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1952 if(nbOf2DCellsToBeSplit==0)
1955 int *retPtr(ret->getPointer());
1956 for(int i=0;i<nCell;i++)
1957 if(cells2DToTreat[i])
1960 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1961 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1962 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1963 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1964 if(middleNeedsToBeUsed)
1965 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1966 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1967 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1968 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.
1969 setPartOfMySelf(ret->begin(),ret->end(),*modif);
1971 bool areNodesMerged; int newNbOfNodes;
1972 if(nbOfNodesCreated!=0)
1973 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
1979 * 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.
1980 * 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).
1981 * 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
1982 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
1983 * 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
1984 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
1986 * 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
1987 * using new instance, idem for coordinates.
1989 * 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.
1991 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
1993 * \throw If \a this is not coherent.
1994 * \throw If \a this has not spaceDim equal to 2.
1995 * \throw If \a this has not meshDim equal to 2.
1997 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
1999 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2001 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2002 checkConsistencyLight();
2003 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2004 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2005 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2006 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2007 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2008 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2009 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2010 const double *coords(_coords->begin());
2011 int *newciptr(newci->getPointer());
2012 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2014 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
2015 ret->pushBackSilent(i);
2016 newciptr[1]=newc->getNumberOfTuples();
2021 if(!appendedCoords->empty())
2023 appendedCoords->rearrange(2);
2024 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2026 setCoords(newCoords);
2029 setConnectivity(newc,newci,true);
2035 * c, cI describe a wire mesh in 3D space, in local numbering
2036 * startNode, endNode in global numbering
2037 *\return true if the segment is indeed split
2039 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2040 const int * c, const int * cI, const int *idsBg, const int *endBg,
2041 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2043 using namespace std;
2045 const int SPACEDIM=3;
2046 typedef pair<double, int> PairDI;
2048 for (const int * it = idsBg; it != endBg; ++it)
2050 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2051 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2052 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2053 x.insert(make_pair(coo[end*SPACEDIM], end));
2056 vector<PairDI> xx(x.begin(), x.end());
2057 sort(xx.begin(),xx.end());
2058 pointIds.reserve(xx.size());
2060 // Keep what is inside [startNode, endNode]:
2062 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2064 const int idx = (*it).second;
2067 if (idx == startNode) go = 1;
2068 if (idx == endNode) go = 2;
2069 if (go) pointIds.push_back(idx);
2072 pointIds.push_back(idx);
2073 if (idx == endNode || idx == startNode)
2077 // vector<int> pointIds2(pointIds.size()+2);
2078 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2079 // pointIds2[0] = startNode;
2080 // pointIds2[pointIds2.size()-1] = endNode;
2083 reverse(pointIds.begin(), pointIds.end());
2085 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2086 for (const int * it = idsBg; it != endBg; ++it)
2088 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2089 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2090 if (itStart == pointIds.end()) continue;
2091 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2092 if (itEnd == pointIds.end()) continue;
2093 if (abs(distance(itEnd, itStart)) != 1) continue;
2094 hitSegs.push_back(*it); // segment is undivided.
2097 return (pointIds.size() > 2); // something else apart start and end node
2100 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2101 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2103 using namespace std;
2104 int dst = distance(sIdxConn, sIdxConnE);
2105 modifiedFace.reserve(dst + insidePoints.size()-2);
2106 modifiedFace.resize(dst);
2107 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2109 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2110 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2111 if (startPos == shortEnd)
2112 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2113 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2114 if (endPos == shortEnd)
2115 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2116 int d = distance(startPos, endPos);
2117 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2118 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2120 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2127 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2128 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2129 * This method performs a conformization of \b this.
2131 * Only polyhedron cells are supported. You can call convertAllToPoly()
2133 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2134 * This method expects that all nodes in \a this are not closer than \a eps.
2135 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2137 * \param [in] eps the relative error to detect merged edges.
2138 * \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
2139 * that the user is expected to deal with.
2141 * \throw If \a this is not coherent.
2142 * \throw If \a this has not spaceDim equal to 3.
2143 * \throw If \a this has not meshDim equal to 3.
2144 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2146 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2148 using namespace std;
2150 static const int SPACEDIM=3;
2151 checkConsistencyLight();
2152 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2153 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2154 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2155 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2157 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2158 const double * coo(_coords->begin());
2159 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2162 /*************************
2164 *************************/
2165 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2166 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2167 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2168 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2169 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2172 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2173 const double *bbox(bboxArr->begin()); getCoords()->begin();
2174 int nDescCell(mDesc->getNumberOfCells());
2175 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2176 // Surfaces - handle biggest first
2177 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2178 DataArrayDouble * surfs = surfF->getArray();
2180 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2181 DataArrayDouble * normals = normalsF->getArray();
2182 const double * normalsP = normals->getConstPointer();
2184 // Sort faces by decreasing surface:
2185 vector< pair<double,int> > S;
2186 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2188 pair<double,int> p = make_pair(surfs->begin()[i], i);
2191 sort(S.rbegin(),S.rend()); // reverse sort
2192 vector<bool> hit(nDescCell);
2193 fill(hit.begin(), hit.end(), false);
2194 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2196 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2198 int faceIdx = (*it).second;
2199 if (hit[faceIdx]) continue;
2201 vector<int> candidates, cands2;
2202 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2203 // Keep only candidates whose normal matches the normal of current face
2204 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2206 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2207 if (*it2 != faceIdx && col)
2208 cands2.push_back(*it2);
2213 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2214 INTERP_KERNEL::TranslationRotationMatrix rotation;
2215 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2216 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2217 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2219 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2220 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2221 double * cooPartRef(mPartRef->_coords->getPointer());
2222 double * cooPartCand(mPartCand->_coords->getPointer());
2223 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2224 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2225 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2226 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2228 // Localize faces in 2D thanks to barycenters
2229 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2230 vector<int> compo; compo.push_back(2);
2231 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2232 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2233 if (!idsGoodPlane->getNumberOfTuples())
2236 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2238 compo[0] = 0; compo.push_back(1);
2239 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2240 mPartRef->changeSpaceDimension(2,0.0);
2241 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2242 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2244 if (!cc->getNumberOfTuples())
2246 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2249 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2250 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2253 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2254 throw INTERP_KERNEL::Exception(oss.str());
2258 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2260 if (!ids->getNumberOfTuples())
2263 double checkSurf=0.0;
2264 const int * idsGoodPlaneP(idsGoodPlane->begin());
2265 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2267 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2268 hit[faceIdx2] = true;
2269 checkSurf += surfs->begin()[faceIdx2];
2271 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2274 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2275 throw INTERP_KERNEL::Exception(oss.str());
2278 // For all polyhedrons using this face, replace connectivity:
2279 vector<int> polyIndices, packsIds, facePack;
2280 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2281 polyIndices.push_back(revDescP[ii]);
2282 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2284 // Current face connectivity
2285 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2286 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2287 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2288 // Deletion of old faces
2290 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2292 if (packsIds[jj] == -1)
2293 // The below should never happen - if a face is used several times, with a different layout of the nodes
2294 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2295 // faces which are actually used only once, by a single cell. This is different for edges below.
2296 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2298 connSla->deletePack(*it2, packsIds[jj]);
2300 // Insertion of new faces:
2301 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2303 // Build pack from the face to insert:
2304 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2306 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2307 // Insert it in all hit polyhedrons:
2308 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2309 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2314 // Set back modified connectivity
2315 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2316 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2317 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2320 /************************
2322 ************************/
2323 // Now we have a face-conform mesh.
2325 // Recompute descending
2326 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2327 // Rebuild desc connectivity with orientation this time!!
2328 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2329 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2330 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2331 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2332 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2333 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2334 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2335 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2336 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2337 // std::cout << "writing!\n";
2338 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2339 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2340 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2341 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2342 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2343 const double *bbox2(bboxArr->begin());
2344 int nDesc2Cell=mDesc2->getNumberOfCells();
2345 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2347 // Edges - handle longest first
2348 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2349 DataArrayDouble * lens = lenF->getArray();
2351 // Sort edges by decreasing length:
2352 vector<pair<double,int> > S;
2353 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2355 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2358 sort(S.rbegin(),S.rend()); // reverse sort
2360 vector<bool> hit(nDesc2Cell);
2361 fill(hit.begin(), hit.end(), false);
2363 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2365 int eIdx = (*it).second;
2369 vector<int> candidates, cands2;
2370 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2371 // Keep only candidates colinear with current edge
2373 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2374 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2375 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2376 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2379 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2380 for (int i3=0; i3 < 3; i3++)
2381 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2382 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2383 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2384 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2386 cands2.push_back(*it2);
2388 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2391 // Now rotate edges to bring them on Ox
2392 int startNode = cDesc2[cIDesc2[eIdx]+1];
2393 int endNode = cDesc2[cIDesc2[eIdx]+2];
2394 INTERP_KERNEL::TranslationRotationMatrix rotation;
2395 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2396 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2397 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2398 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2401 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2402 nbElemsNotM1 = tmp->getNbOfElems();
2404 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2405 double * cooPartRef(mPartRef->_coords->getPointer());
2406 double * cooPartCand(mPartCand->_coords->getPointer());
2407 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2408 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2409 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2410 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2413 // Eliminate all edges for which y or z is not null
2414 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2415 vector<int> compo; compo.push_back(1);
2416 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2418 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2419 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2420 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2421 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2422 if (!idsGoodLine->getNumberOfTuples())
2425 // Now the ordering along the Ox axis:
2426 std::vector<int> insidePoints, hitSegs;
2427 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2428 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2429 idsGoodLine->begin(), idsGoodLine->end(),
2430 /*out*/insidePoints, hitSegs);
2431 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2432 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2433 hit[cands2[*its]] = true;
2435 if (!isSplit) // current segment remains in one piece
2438 // Get original node IDs in global coords array
2439 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2440 *iit = nodeMapInv->begin()[*iit];
2442 vector<int> polyIndices, packsIds, facePack;
2443 // For each face implying this edge
2444 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2446 int faceIdx = revDescP2[ii];
2447 // each cell where this face is involved connectivity will be modified:
2448 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2450 // Current face connectivity
2451 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2452 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2454 vector<int> modifiedFace;
2455 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2456 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2457 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2461 // Rebuild 3D connectivity from descending:
2462 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2463 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2464 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2465 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2466 newConn->set3(superIdx, idx, vals);
2467 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2468 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2470 int sz, faceIdx = abs(descP[jj])-1;
2471 bool orient = descP[jj]>0;
2472 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2474 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2477 vector<int> rev(sz-1);
2478 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2479 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2483 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2486 ret = ret->buildUniqueNotSorted();