1 // Copyright (C) 2007-2019 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, INTERP_KERNEL::NodeWithUsage >& mapp2, const int *bg)
203 INTERP_KERNEL::Edge *ret=0;
205 mapp2[bg[0]].second = INTERP_KERNEL::USAGE_LINEAR;
206 mapp2[bg[1]].second = INTERP_KERNEL::USAGE_LINEAR;
210 case INTERP_KERNEL::NORM_SEG2:
212 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
215 case INTERP_KERNEL::NORM_SEG3:
217 INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
218 INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
219 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
220 // is the SEG3 degenerated, and thus can be reduced to a SEG2?
221 bool colinearity=inters.areColinears();
222 delete e1; delete e2;
224 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
226 ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
227 if (mapp2[bg[2]].second != INTERP_KERNEL::USAGE_LINEAR) // switch the node usage to quadratic only if it is not used as an extreme point for another edge
228 mapp2[bg[2]].second = INTERP_KERNEL::USAGE_QUADRATIC_ONLY;
232 throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
238 * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from
239 * the global mesh 'mDesc'.
240 * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2.
241 * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
243 INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
244 std::map<INTERP_KERNEL::Node *,int>& mapp)
247 std::map<int, INTERP_KERNEL::NodeWithUsage > mapp2; // the last var is a flag specifying if node is an extreme node of the seg (LINEAR) or only a middle for SEG3 (QUADRATIC_ONLY).
248 const double *coo=mDesc->getCoords()->getConstPointer();
249 const int *c=mDesc->getNodalConnectivity()->getConstPointer();
250 const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
252 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
253 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
254 for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
256 INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
257 mapp2[*it2]=INTERP_KERNEL::NodeWithUsage(n,INTERP_KERNEL::USAGE_UNKNOWN);
259 INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
260 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
262 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
263 ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
265 for(std::map<int, INTERP_KERNEL::NodeWithUsage >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
267 if((*it2).second.second == INTERP_KERNEL::USAGE_LINEAR)
268 mapp[(*it2).second.first]=(*it2).first;
269 ((*it2).second.first)->decrRef();
274 INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
278 int locId=nodeId-offset2;
279 return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
283 int locId=nodeId-offset1;
284 return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
286 return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
290 * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
292 void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
293 const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
294 /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
296 for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
298 int eltId1=abs(*desc1)-1;
299 for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
301 std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
302 if(it==mappRev.end())
304 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
316 * 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 ) .
317 * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
318 * \param forbiddenPoints the list of points that should not be removed in the process
320 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset,
321 const std::map<int, bool>& forbiddenPoints,
322 DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
324 std::size_t sz(std::distance(connBg,connEnd));
325 if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
326 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
328 INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
329 INTERP_KERNEL::AutoPtr<int> tmpConn2(new int[sz]);
330 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
331 unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
332 unsigned nbOfHit(0); // number of fusions operated
333 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
334 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
335 INTERP_KERNEL::NormalizedCellType typeOfSon;
336 std::vector<int> middles;
338 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
340 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
341 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
342 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
343 posEndElt = posBaseElt+1;
345 // Look backward first: are the final edges of the cells colinear with the first ones?
346 // This initializes posBaseElt.
349 for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
351 cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn2,typeOfSon);
352 // Identify common point:
353 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
354 auto itE(forbiddenPoints.end());
355 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
357 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
358 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
359 bool isColinear=eint->areColinears();
370 // Update last connectivity
371 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
375 const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
376 for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell with one single edge
378 cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn2,typeOfSon); // get edge #j's connectivity
379 // Identify common point:
380 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
381 auto itE(forbiddenPoints.end());
382 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
384 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
385 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
386 bool isColinear(eint->areColinears());
397 // Update last connectivity
398 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
400 //push [posBaseElt,posEndElt) in newConnOfCell using e
401 // 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!
403 // at the beginning of the connectivity (insert type)
404 EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
405 else if((nbOfHit+nbOfTurn) != (nbs-1))
407 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
408 if ((nbOfHit+nbOfTurn) == (nbs-1))
409 // at the end (only quad points to deal with)
410 EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
411 posBaseElt=posEndElt;
415 newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
421 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
423 if(candidates.empty())
425 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
427 const std::vector<int>& pool(intersectEdge1[*it]);
428 int tmp[2]; tmp[0]=start; tmp[1]=stop;
429 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
434 tmp[0]=stop; tmp[1]=start;
435 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
445 * This method performs the 2nd step of Partition of 2D mesh.
446 * This method has 4 inputs :
447 * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
448 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
449 * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
450 * 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'
451 * Nodes end up lying consecutively on a cutted edge.
452 * \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.
453 * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
454 * \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.
455 * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
456 * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
458 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
459 const std::vector<double>& addCoo,
460 const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
462 int offset1=m1->getNumberOfNodes();
463 int ncell2=m2->getNumberOfCells();
464 const int *c=m2->getNodalConnectivity()->begin();
465 const int *cI=m2->getNodalConnectivityIndex()->begin();
466 const double *coo=m2->getCoords()->begin();
467 const double *cooBis=m1->getCoords()->begin();
468 int offset2=offset1+m2->getNumberOfNodes();
469 intersectEdge.resize(ncell2);
470 for(int i=0;i<ncell2;i++,cI++)
472 const std::vector<int>& divs=subDiv[i];
473 int nnode=cI[1]-cI[0]-1;
474 std::map<int, INTERP_KERNEL::NodeWithUsage > mapp2;
475 std::map<INTERP_KERNEL::Node *, int> mapp22;
476 for(int j=0;j<nnode;j++)
478 INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
479 int nnid=c[(*cI)+j+1];
480 mapp2[nnid]=INTERP_KERNEL::NodeWithUsage(nn,INTERP_KERNEL::USAGE_UNKNOWN);
481 mapp22[nn]=nnid+offset1;
483 INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
484 for(std::map<int, INTERP_KERNEL::NodeWithUsage >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
485 ((*it).second.first)->decrRef();
486 std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
487 std::map<INTERP_KERNEL::Node *,int> mapp3;
488 for(std::size_t j=0;j<divs.size();j++)
491 INTERP_KERNEL::Node *tmp=0;
493 tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
495 tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
497 tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
501 e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
502 for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
508 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,
509 MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
511 idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
512 idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
513 int nCells(mesh1D->getNumberOfCells());
514 if(nCells!=(int)intersectEdge2.size())
515 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
516 const DataArrayDouble *coo2(mesh1D->getCoords());
517 const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
518 const double *coo2Ptr(coo2->begin());
519 int offset1(coords1->getNumberOfTuples());
520 int offset2(offset1+coo2->getNumberOfTuples());
521 int offset3(offset2+addCoo.size()/2);
522 std::vector<double> addCooQuad;
523 MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
524 int tmp[4],cicnt(0),kk(0);
525 for(int i=0;i<nCells;i++)
527 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
528 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
529 const std::vector<int>& subEdges(intersectEdge2[i]);
530 int nbSubEdge(subEdges.size()/2);
531 for(int j=0;j<nbSubEdge;j++,kk++)
533 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));
534 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
535 INTERP_KERNEL::Edge *e2Ptr(e2);
536 std::map<int,int>::const_iterator itm;
537 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
539 tmp[0]=INTERP_KERNEL::NORM_SEG3;
540 itm=mergedNodes.find(subEdges[2*j]);
541 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
542 itm=mergedNodes.find(subEdges[2*j+1]);
543 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
544 tmp[3]=offset3+(int)addCooQuad.size()/2;
546 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
548 cOut->insertAtTheEnd(tmp,tmp+4);
549 ciOut->pushBackSilent(cicnt);
553 tmp[0]=INTERP_KERNEL::NORM_SEG2;
554 itm=mergedNodes.find(subEdges[2*j]);
555 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
556 itm=mergedNodes.find(subEdges[2*j+1]);
557 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
559 cOut->insertAtTheEnd(tmp,tmp+3);
560 ciOut->pushBackSilent(cicnt);
563 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
565 idsInRetColinear->pushBackSilent(kk);
566 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
571 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
572 ret->setConnectivity(cOut,ciOut,true);
573 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
574 arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
575 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,(int)addCooQuad.size()/2,2);
576 std::vector<const DataArrayDouble *> coordss(4);
577 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
578 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
583 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
585 std::vector<int> allEdges;
586 for(const int *it2(descBg);it2!=descEnd;it2++)
588 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
590 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
592 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
594 std::size_t nb(allEdges.size());
596 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
597 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
598 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
599 ret->setCoords(coords);
600 ret->allocateCells(1);
601 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
602 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
603 connOut[kk]=allEdges[2*kk];
604 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
608 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
610 const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
611 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
613 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
614 if(sz!=std::distance(descBg,descEnd))
615 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
616 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
617 std::vector<int> allEdges,centers;
618 const double *coordsPtr(coords->begin());
619 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
620 int offset(coords->getNumberOfTuples());
621 for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
623 INTERP_KERNEL::NormalizedCellType typeOfSon;
624 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
625 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
627 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
629 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
631 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
633 {//the current edge has been subsplit -> create corresponding centers.
634 std::size_t nbOfCentersToAppend(edge1.size()/2);
635 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
636 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
637 std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
638 for(std::size_t k=0;k<nbOfCentersToAppend;k++)
641 const double *aa(coordsPtr+2*(*it3++));
642 const double *bb(coordsPtr+2*(*it3++));
643 ee->getMiddleOfPoints(aa,bb,tmpp);
644 addCoo->insertAtTheEnd(tmpp,tmpp+2);
645 centers.push_back(offset+k);
649 std::size_t nb(allEdges.size());
651 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
652 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
653 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
655 ret->setCoords(coords);
658 addCoo->rearrange(2);
659 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
660 ret->setCoords(addCoo);
662 ret->allocateCells(1);
663 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
664 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
665 connOut[kk]=allEdges[2*kk];
666 connOut.insert(connOut.end(),centers.begin(),centers.end());
667 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
672 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
675 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
677 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
679 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
680 if(!cm.isQuadratic())
681 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
683 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
686 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
689 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
691 const INTERP_KERNEL::Edge *ee(*it);
692 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
696 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
699 const double *coo(mesh2D->getCoords()->begin());
700 std::size_t sz(conn.size());
701 std::vector<double> addCoo;
702 std::vector<int> conn2(conn);
703 int offset(mesh2D->getNumberOfNodes());
704 for(std::size_t i=0;i<sz;i++)
707 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
708 addCoo.insert(addCoo.end(),tmp,tmp+2);
709 conn2.push_back(offset+(int)i);
711 mesh2D->getCoords()->rearrange(1);
712 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
713 mesh2D->getCoords()->rearrange(2);
714 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
719 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
721 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
722 * a set of edges defined in \a splitMesh1D.
724 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
725 std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
727 std::size_t nb(edge1Bis.size()/2);
728 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
729 int iEnd(splitMesh1D->getNumberOfCells());
731 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
733 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
734 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
735 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
738 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
739 out0.resize(1); out1.resize(1);
740 std::vector<int>& connOut(out0[0]);
741 connOut.resize(nbOfEdgesOf2DCellSplit);
742 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
743 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
744 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
746 connOut[kk]=edge1Bis[2*kk];
747 edgesPtr[kk]=edge1BisPtr[2*kk];
752 // [i,iEnd[ contains the
753 out0.resize(2); out1.resize(2);
754 std::vector<int>& connOutLeft(out0[0]);
755 std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
756 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
757 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
758 for(std::size_t k=ii;k<jj+1;k++)
759 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
760 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
761 for(int ik=0;ik<iEnd;ik++)
763 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
764 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
767 for(int ik=iEnd-1;ik>=0;ik--)
768 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
769 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
770 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
771 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
772 for(int ik=0;ik<iEnd;ik++)
773 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
774 eright.insert(eright.end(),ees.begin(),ees.end());
782 CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
784 std::vector<int> _edges;
785 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
788 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
790 std::size_t nbe(edges.size());
791 std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
792 for(std::size_t i=0;i<nbe;i++)
794 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
795 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
797 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
798 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
799 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
805 EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
806 EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
807 bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
808 void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
809 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
813 MCAuto<MEDCouplingUMesh> _mesh;
814 MCAuto<INTERP_KERNEL::Edge> _edge;
819 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
821 const MEDCouplingUMesh *mesh(_mesh);
827 { _left++; _right++; return ; }
830 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
831 if((isLeft && isRight) || (!isLeft && !isRight))
832 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
843 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
844 if((isLeft && isRight) || (!isLeft && !isRight))
845 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
860 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
862 const MEDCouplingUMesh *mesh(_mesh);
865 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
868 {// not fully splitting cell case
869 if(mesh2D->getNumberOfCells()==1)
870 {//little optimization. 1 cell no need to find in which cell mesh is !
871 neighbors[0]=offset; neighbors[1]=offset;
876 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
877 int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
879 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
880 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
885 class VectorOfCellInfo
888 VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
889 std::size_t size() const { return _pool.size(); }
890 int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
891 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);
892 const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
893 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
894 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
895 void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
897 int getZePosOfEdgeGivenItsGlobalId(int pos) const;
898 void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
899 const CellInfo& get(int pos) const;
900 CellInfo& get(int pos);
902 std::vector<CellInfo> _pool;
903 MCAuto<MEDCouplingUMesh> _ze_mesh;
904 std::vector<EdgeInfo> _edge_info;
907 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
909 _pool[0]._edges=edges;
910 _pool[0]._edges_ptr=edgesPtr;
913 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
916 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
919 const MEDCouplingUMesh *zeMesh(_ze_mesh);
921 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
922 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
923 return zeMesh->getCellContainingPoint(barys->begin(),eps);
926 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
927 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
928 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
930 get(pos);//to check pos
931 bool isFast(pos==0 && _pool.size()==1);
932 std::size_t sz(edges.size());
933 // dealing with edges
935 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
937 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
939 std::vector<CellInfo> pool(_pool.size()-1+sz);
940 for(std::size_t i=0;i<pos;i++)
942 for(std::size_t j=0;j<sz;j++)
943 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
944 for(int i=pos+1;i<(int)_pool.size();i++)
945 pool[i+sz-1]=_pool[i];
949 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
957 std::vector< MCAuto<MEDCouplingUMesh> > ms;
960 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
964 if(pos<_ze_mesh->getNumberOfCells()-1)
966 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
969 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
970 for(std::size_t j=0;j<ms2.size();j++)
972 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
975 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
977 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
980 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
983 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
985 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
987 if((*it).isInMyRange(pos))
990 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
993 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
995 get(pos);//to perform the sanity check;
996 if(_edge_info.empty())
998 std::size_t sz(_edge_info.size()-1);
999 for(std::size_t i=0;i<sz;i++)
1000 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
1003 const CellInfo& VectorOfCellInfo::get(int pos) const
1005 if(pos<0 || pos>=(int)_pool.size())
1006 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1010 CellInfo& VectorOfCellInfo::get(int pos)
1012 if(pos<0 || pos>=(int)_pool.size())
1013 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1019 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1020 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1022 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1024 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1026 * \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.
1028 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1029 MCAuto<DataArrayInt>& idsLeftRight)
1031 int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1032 if(nbCellsInSplitMesh1D==0)
1033 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1034 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1035 std::size_t nb(allEdges.size()),jj;
1037 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1038 std::vector<int> edge1Bis(nb*2);
1039 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1040 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1041 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1042 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1043 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1045 idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1046 int *idsLeftRightPtr(idsLeftRight->getPointer());
1047 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1049 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1050 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1051 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1052 // of the connectivity.
1053 MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1054 renumb->alloc(nbCellsInSplitMesh1D,1);
1055 const int * renumbP(renumb->begin());
1057 int i, first=cSplitPtr[1];
1058 // Follow 1D line backward as long as it is connected:
1059 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1060 first=cSplitPtr[ciSplitPtr[i]+1];
1061 if (i < nbCellsInSplitMesh1D-1)
1063 // Build circular permutation to shift consecutive edges together
1065 renumb->applyModulus(nbCellsInSplitMesh1D);
1066 splitMesh1D->renumberCells(renumbP, false);
1067 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1068 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1074 for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1075 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1077 for(;iEnd<nbCellsInSplitMesh1D;)
1079 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1085 if(iEnd<nbCellsInSplitMesh1D)
1088 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1089 int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1091 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1092 retTmp->setCoords(splitMesh1D->getCoords());
1093 retTmp->allocateCells();
1095 std::vector< std::vector<int> > out0;
1096 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1098 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1099 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1100 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1101 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1105 for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1106 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1108 return pool.getZeMesh().retn();
1112 * splitMesh1D is an input parameter but might have its cells renumbered.
1114 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1115 const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1116 MCAuto<DataArrayInt>& idsLeftRight)
1118 const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1120 std::vector<int> allEdges;
1121 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1122 for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1124 int edgeId(std::abs(*it)-1);
1125 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1126 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1127 const std::vector<int>& edge1(intersectEdge1[edgeId]);
1129 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1131 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1132 std::size_t sz(edge1.size());
1133 for(std::size_t cnt=0;cnt<sz;cnt++)
1134 allEdgesPtr.push_back(ee);
1137 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1140 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1142 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1143 {//easy case comparison not
1144 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1146 else if(typ1.isQuadratic() && typ2.isQuadratic())
1148 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1151 if(conn1[2]==conn2[2])
1153 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1154 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1158 {//only one is quadratic
1159 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1162 const double *a(0),*bb(0),*be(0);
1163 if(typ1.isQuadratic())
1165 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1169 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1171 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1172 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1178 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1179 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1181 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1183 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1185 if(candidatesIn2DEnd==candidatesIn2DBg)
1186 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1187 const double *coo(mesh2DSplit->getCoords()->begin());
1188 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1189 return *candidatesIn2DBg;
1190 int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1191 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1192 if(cellIdInMesh1DSplitRelative<0)
1193 cur1D->changeOrientationOfCells();
1194 const int *c1D(cur1D->getNodalConnectivity()->begin());
1195 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1196 for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1198 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1199 const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1200 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1201 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1202 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1203 for(unsigned it2=0;it2<sz;it2++)
1205 INTERP_KERNEL::NormalizedCellType typeOfSon;
1206 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1207 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1208 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1212 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1216 * \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.
1217 * 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.
1218 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1219 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1220 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1221 * \param [out] addCoo - nodes to be append at the end
1222 * \param [out] mergedNodes - gives all pair of nodes of \a m2Desc that have same location than some nodes in \a m1Desc. key is id in \a m2Desc offsetted and value is id in \a m1Desc.
1224 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1225 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)
1227 static const int SPACEDIM=2;
1228 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1229 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1230 // Build BB tree of all edges in the tool mesh (second mesh)
1231 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1232 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1233 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1234 intersectEdge1.resize(nDescCell1);
1235 colinear2.resize(nDescCell2);
1236 subDiv2.resize(nDescCell2);
1237 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1239 std::vector<int> candidates1(1);
1240 int offset1(m1Desc->getNumberOfNodes());
1241 int offset2(offset1+m2Desc->getNumberOfNodes());
1242 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1244 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1245 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1246 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1248 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1249 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1250 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1252 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1253 // 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
1254 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1255 std::set<INTERP_KERNEL::Node *> nodes;
1256 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1257 std::size_t szz(nodes.size());
1258 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1259 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1260 for(std::size_t iii=0;iii<szz;iii++,itt++)
1261 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1262 // end of protection
1263 // Performs edge cutting:
1264 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1269 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1270 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1276 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1277 * It builds the descending connectivity of the two meshes, and then using a binary tree
1278 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1279 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1281 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1282 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1283 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1284 std::vector<double>& addCoo,
1285 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1287 // Build desc connectivity
1288 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1289 desc2=DataArrayInt::New();
1290 descIndx2=DataArrayInt::New();
1291 revDesc2=DataArrayInt::New();
1292 revDescIndx2=DataArrayInt::New();
1293 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1294 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1295 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1296 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1297 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1298 std::map<int,int> notUsedMap;
1299 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1300 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1301 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1305 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1306 * (newly created) nodes corresponding to the edge intersections.
1308 * @param[out] cr, crI connectivity of the resulting mesh
1309 * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1310 * TODO: describe input parameters
1312 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1313 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1314 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1315 const std::vector<double>& addCoords,
1316 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1318 static const int SPACEDIM=2;
1319 const double *coo1(m1->getCoords()->begin());
1320 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1321 int offset1(m1->getNumberOfNodes());
1322 const double *coo2(m2->getCoords()->begin());
1323 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1324 int offset2(offset1+m2->getNumberOfNodes());
1325 int offset3(offset2+((int)addCoords.size())/2);
1326 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1327 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1328 // Here a BBTree on 2D-cells, not on segments:
1329 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1330 int ncell1(m1->getNumberOfCells());
1332 for(int i=0;i<ncell1;i++)
1334 std::vector<int> candidates2;
1335 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1336 std::map<INTERP_KERNEL::Node *,int> mapp;
1337 std::map<int,INTERP_KERNEL::Node *> mappRev;
1338 INTERP_KERNEL::QuadraticPolygon pol1;
1339 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1340 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1341 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1342 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1343 // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1344 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1345 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1347 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
1348 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1349 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1350 for(it1.first();!it1.finished();it1.next())
1351 edges1.insert(it1.current()->getPtr());
1353 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1354 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1356 // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1357 // Again all the additional intersecting nodes are there.
1358 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1360 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1361 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1362 // Complete mapping with elements coming from the current cell it2 in mesh2:
1363 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1364 // pol2 is the new QP in the final merged result.
1365 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1366 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1368 // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1369 for (auto &p: pol2s)
1370 p.cleanDegeneratedConsecutiveEdges();
1371 edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges.
1374 // Now rebuild intersected cells from all this:
1375 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1377 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1378 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1379 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1380 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1382 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1383 // by m2 but that we still want to keep in the final result.
1388 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1390 catch(INTERP_KERNEL::Exception& e)
1392 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();
1393 throw INTERP_KERNEL::Exception(oss.str());
1396 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1397 (*it).second->decrRef();
1401 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1403 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1406 std::size_t sz(conn.size());
1407 std::size_t found(std::numeric_limits<std::size_t>::max());
1408 for(std::size_t i=0;i<sz;i++)
1410 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1411 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]};
1412 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1413 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1414 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1416 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];
1417 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]);
1419 if(dotTest>eps && dotTest<1.-eps)
1425 if(found==std::numeric_limits<std::size_t>::max())
1426 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1427 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1430 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1432 std::size_t sz(conn.size());
1433 std::vector<int> *curPart(&part0);
1434 for(std::size_t i=0;i<sz;i++)
1436 int nextt(conn[(i+1)%sz]);
1437 (*curPart).push_back(nextt);
1438 if(nextt==pt0 || nextt==pt1)
1444 (*curPart).push_back(nextt);
1450 * 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.
1452 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1453 const int *desc, const int *descIndx, const double *coords, double eps,
1454 std::vector<std::vector<int> >& res) const
1456 checkFullyDefined();
1457 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1458 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1459 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1460 int nbOfCells(getNumberOfCells());
1462 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1463 for(int i=0;i<nbOfCells;i++)
1465 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1466 for(int j=0;j<nbOfFaces;j++)
1468 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1469 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1470 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1471 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1472 INTERP_KERNEL::NormalizedCellType cmsId;
1473 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1474 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1475 if(p.first!=-1 && p.second!=-1)
1479 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1480 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1481 std::vector<int> elt1,elt2;
1482 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1483 res.push_back(elt1);
1484 res.push_back(elt2);
1496 * 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).
1498 * \sa MEDCouplingUMesh::split2DCells
1500 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1502 checkConnectivityFullyDefined();
1503 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1504 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1505 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1506 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1507 int prevPosOfCi(ciPtr[0]);
1508 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1510 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1511 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1512 for(int j=0;j<sz;j++)
1514 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1515 for(int k=0;k<sz2;k++)
1516 *cPtr++=subPtr[offset2+k];
1518 *cPtr++=oldConn[prevPosOfCi+j+2];
1521 prevPosOfCi=ciPtr[1];
1522 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1525 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1526 _nodal_connec->decrRef();
1527 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1532 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1534 * \return int - the number of new nodes created.
1535 * \sa MEDCouplingUMesh::split2DCells
1537 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1539 checkConsistencyLight();
1540 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1541 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1542 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1543 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1544 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1545 const double *oldCoordsPtr(getCoords()->begin());
1546 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1547 int prevPosOfCi(ciPtr[0]);
1548 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1550 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1551 for(int j=0;j<sz;j++)
1552 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1553 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1554 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1556 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1560 cPtr[1]=oldConn[prevPosOfCi+2+j];
1561 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1564 std::vector<INTERP_KERNEL::Node *> ns(3);
1565 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1566 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1567 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1568 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1569 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1571 cPtr[1]=subPtr[offset2+k];
1572 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1574 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1577 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1579 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1580 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1583 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1584 _nodal_connec->decrRef();
1585 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1586 addCoo->rearrange(2);
1587 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1589 return addCoo->getNumberOfTuples();
1596 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1597 * returns a result mesh constituted by polygons.
1598 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1599 * all nodes from m2.
1600 * The meshes should be in 2D space. In
1601 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1603 * \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
1604 * 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)
1605 * \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
1606 * 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)
1607 * \param [in] eps - precision used to detect coincident mesh entities.
1608 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1609 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1610 * this array using decrRef() as it is no more needed.
1611 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1612 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1613 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1614 * any cell of \a m2. The caller is to delete this array using decrRef() as
1615 * it is no more needed.
1616 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1617 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1618 * is no more needed.
1619 * \throw If the coordinates array is not set in any of the meshes.
1620 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1621 * \throw If any of the meshes is not a 2D mesh in 2D space.
1623 * \sa conformize2D, mergeNodes
1625 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1626 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1629 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1630 m1->checkFullyDefined();
1631 m2->checkFullyDefined();
1632 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1633 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1634 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1636 // Step 1: compute all edge intersections (new nodes)
1637 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1638 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1639 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1640 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1641 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1642 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1643 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1644 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1645 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1646 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1648 // Step 2: re-order newly created nodes according to the ordering found in m2
1649 std::vector< std::vector<int> > intersectEdge2;
1650 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1651 subDiv2.clear(); dd5=0; dd6=0;
1654 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1655 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1656 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1657 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1659 // Step 4: Prepare final result:
1660 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1661 addCooDa->alloc((int)(addCoo.size())/2,2);
1662 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1663 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1664 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1665 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1666 std::vector<const DataArrayDouble *> coordss(4);
1667 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1668 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1669 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1670 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1671 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1672 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1673 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1674 ret->setConnectivity(conn,connI,true);
1675 ret->setCoords(coo);
1676 cellNb1=c1.retn(); cellNb2=c2.retn();
1681 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1682 * 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
1683 * and finally, in case of quadratic polygon the centers of edges new nodes.
1684 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1686 * \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
1687 * 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)
1688 * \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
1689 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1690 * \param [in] eps - precision used to perform intersections and localization operations.
1691 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1692 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1693 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1694 * 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.
1695 * \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
1696 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1697 * 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.
1699 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1701 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1703 if(!mesh2D || !mesh1D)
1704 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1705 mesh2D->checkFullyDefined();
1706 mesh1D->checkFullyDefined();
1707 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1708 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1709 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1710 // Step 1: compute all edge intersections (new nodes)
1711 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1712 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1713 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1715 // Build desc connectivity
1716 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1717 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1718 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1719 std::map<int,int> mergedNodes;
1720 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1721 // use mergeNodes to fix intersectEdge1
1722 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1724 std::size_t n((*it0).size()/2);
1725 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1726 std::map<int,int>::const_iterator it1;
1727 it1=mergedNodes.find(eltStart);
1728 if(it1!=mergedNodes.end())
1729 (*it0)[0]=(*it1).second;
1730 it1=mergedNodes.find(eltEnd);
1731 if(it1!=mergedNodes.end())
1732 (*it0)[2*n-1]=(*it1).second;
1735 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1736 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
1737 // Step 2: re-order newly created nodes according to the ordering found in m2
1738 std::vector< std::vector<int> > intersectEdge2;
1739 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1741 // Step 3: compute splitMesh1D
1742 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1743 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1744 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1745 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1746 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1747 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1748 // deal with cells in mesh2D that are not cut but only some of their edges are
1749 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1750 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1751 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1752 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
1753 if(!idsInDesc2DToBeRefined->empty())
1755 DataArrayInt *out0(0),*outi0(0);
1756 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1757 MCAuto<DataArrayInt> outi0s(outi0);
1759 out0s=out0s->buildUnique();
1763 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1764 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1765 MCAuto<DataArrayInt> elts,eltsIndex;
1766 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1767 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1768 if (eltsIndex->getNumberOfTuples() > 1)
1769 eltsIndex2 = eltsIndex->deltaShiftIndex();
1770 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1771 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1772 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1773 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1774 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1775 if((DataArrayInt *)out0s)
1776 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1777 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1778 // OK all is ready to insert in ret2 mesh
1779 if(!untouchedCells->empty())
1780 {// the most easy part, cells in mesh2D not impacted at all
1781 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1782 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1783 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1785 if((DataArrayInt *)out0s)
1786 {// here dealing with cells in out0s but not in cellsToBeModified
1787 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1788 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1789 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1791 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1792 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1794 int offset(ret2->getNumberOfTuples());
1795 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1796 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1797 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1798 int kk(0),*ret3ptr(partOfRet3->getPointer());
1799 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1801 int faceId(std::abs(*it)-1);
1802 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1804 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1807 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1808 ret3ptr[2*kk]=tmp+offset;
1809 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1810 ret3ptr[2*kk+1]=tmp+offset;
1813 {//the current edge is shared by a 2D cell that will be split just after
1814 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1815 ret3ptr[2*kk]=-(*it2+1);
1816 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1817 ret3ptr[2*kk+1]=-(*it2+1);
1821 m1Desc->setCoords(ret1->getCoords());
1822 ret1NonCol->setCoords(ret1->getCoords());
1823 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1824 if(!outMesh2DSplit.empty())
1826 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1827 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1828 (*itt)->setCoords(da);
1831 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1832 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1834 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1835 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1836 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1837 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1838 MCAuto<DataArrayInt> partOfRet3;
1839 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));
1840 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1841 outMesh2DSplit.push_back(splitOfOneCell);
1842 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1843 ret2->pushBackSilent(*it);
1846 std::size_t nbOfMeshes(outMesh2DSplit.size());
1847 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1848 for(std::size_t i=0;i<nbOfMeshes;i++)
1849 tmp[i]=outMesh2DSplit[i];
1851 ret1->getCoords()->setInfoOnComponents(compNames);
1852 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1853 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1855 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1856 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1858 int old2DCellId(-ret3->getIJ(*it,0)-1);
1859 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1860 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
1862 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1865 splitMesh1D=ret1.retn();
1866 splitMesh2D=ret2D.retn();
1867 cellIdInMesh2D=ret2.retn();
1868 cellIdInMesh1D=ret3.retn();
1872 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1873 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1874 * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this method
1875 * 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).
1876 * 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.
1878 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1879 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1881 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1882 * This method expects that all nodes in \a this are not closer than \a eps.
1883 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1885 * \param [in] eps the relative error to detect merged edges.
1886 * \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
1887 * that the user is expected to deal with.
1889 * \throw If \a this is not coherent.
1890 * \throw If \a this has not spaceDim equal to 2.
1891 * \throw If \a this has not meshDim equal to 2.
1892 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1894 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1896 static const int SPACEDIM=2;
1897 checkConsistencyLight();
1898 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1899 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1900 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1901 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1902 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1903 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1904 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1905 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1906 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1907 std::vector<double> addCoo;
1908 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1909 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1910 for(int i=0;i<nDescCell;i++)
1912 std::vector<int> candidates;
1913 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1914 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1915 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1917 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1918 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1919 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1920 INTERP_KERNEL::MergePoints merge;
1921 INTERP_KERNEL::QuadraticPolygon c1,c2;
1922 e1->intersectWith(e2,merge,c1,c2);
1923 e1->decrRef(); e2->decrRef();
1924 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1925 overlapEdge[i].push_back(*it);
1926 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1927 overlapEdge[*it].push_back(i);
1930 // splitting done. sort intersect point in intersectEdge.
1931 std::vector< std::vector<int> > middle(nDescCell);
1932 int nbOf2DCellsToBeSplit(0);
1933 bool middleNeedsToBeUsed(false);
1934 std::vector<bool> cells2DToTreat(nDescCell,false);
1935 for(int i=0;i<nDescCell;i++)
1937 std::vector<int>& isect(intersectEdge[i]);
1938 int sz((int)isect.size());
1941 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1942 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1943 e->sortSubNodesAbs(coords,isect);
1948 int idx0(rdi[i]),idx1(rdi[i+1]);
1950 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1951 if(!cells2DToTreat[rd[idx0]])
1953 cells2DToTreat[rd[idx0]]=true;
1954 nbOf2DCellsToBeSplit++;
1956 // try to reuse at most eventual 'middle' of SEG3
1957 std::vector<int>& mid(middle[i]);
1958 mid.resize(sz+1,-1);
1959 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1961 middleNeedsToBeUsed=true;
1962 const std::vector<int>& candidates(overlapEdge[i]);
1963 std::vector<int> trueCandidates;
1964 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1965 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1966 trueCandidates.push_back(*itc);
1967 int stNode(c[ci[i]+1]),endNode(isect[0]);
1968 for(int j=0;j<sz+1;j++)
1970 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1972 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1973 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1974 { mid[j]=*itc; break; }
1977 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1982 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1983 if(nbOf2DCellsToBeSplit==0)
1986 int *retPtr(ret->getPointer());
1987 for(int i=0;i<nCell;i++)
1988 if(cells2DToTreat[i])
1991 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1992 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1993 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1994 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1995 if(middleNeedsToBeUsed)
1996 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1997 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1998 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1999 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.
2000 setPartOfMySelf(ret->begin(),ret->end(),*modif);
2002 bool areNodesMerged; int newNbOfNodes;
2003 if(nbOfNodesCreated!=0)
2004 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2010 * 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.
2011 * 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).
2012 * 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
2013 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2014 * 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
2015 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2017 * 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
2018 * using new instance, idem for coordinates.
2020 * 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.
2022 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
2024 * \throw If \a this is not coherent.
2025 * \throw If \a this has not spaceDim equal to 2.
2026 * \throw If \a this has not meshDim equal to 2.
2028 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2030 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2032 return internalColinearize2D(eps, false);
2036 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2037 * In a given 2D cell, if two edges are colinear and the junction point is used by a third edge, the two edges will not be
2038 * merged, contrary to colinearize2D().
2040 * \sa MEDCouplingUMesh::colinearize2D
2042 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2044 return internalColinearize2D(eps, true);
2049 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2051 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2053 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2054 checkConsistencyLight();
2055 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2056 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2057 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2058 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2059 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2060 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2061 std::map<int, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2064 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2065 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2066 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2067 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2068 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2069 MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2070 MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2071 const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2072 for(auto it = ids->begin(); it != ids->end(); it++)
2073 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2076 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2077 const double *coords(_coords->begin());
2078 int *newciptr(newci->getPointer());
2079 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2081 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2082 ret->pushBackSilent(i);
2083 newciptr[1]=newc->getNumberOfTuples();
2088 if(!appendedCoords->empty())
2090 appendedCoords->rearrange(2);
2091 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2093 setCoords(newCoords);
2096 setConnectivity(newc,newci,true);
2102 * c, cI describe a wire mesh in 3D space, in local numbering
2103 * startNode, endNode in global numbering
2104 *\return true if the segment is indeed split
2106 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2107 const int * c, const int * cI, const int *idsBg, const int *endBg,
2108 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2110 using namespace std;
2112 const int SPACEDIM=3;
2113 typedef pair<double, int> PairDI;
2115 for (const int * it = idsBg; it != endBg; ++it)
2117 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2118 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2119 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2120 x.insert(make_pair(coo[end*SPACEDIM], end));
2123 vector<PairDI> xx(x.begin(), x.end());
2124 sort(xx.begin(),xx.end());
2125 pointIds.reserve(xx.size());
2127 // Keep what is inside [startNode, endNode]:
2129 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2131 const int idx = (*it).second;
2134 if (idx == startNode) go = 1;
2135 if (idx == endNode) go = 2;
2136 if (go) pointIds.push_back(idx);
2139 pointIds.push_back(idx);
2140 if (idx == endNode || idx == startNode)
2144 // vector<int> pointIds2(pointIds.size()+2);
2145 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2146 // pointIds2[0] = startNode;
2147 // pointIds2[pointIds2.size()-1] = endNode;
2150 reverse(pointIds.begin(), pointIds.end());
2152 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2153 for (const int * it = idsBg; it != endBg; ++it)
2155 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2156 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2157 if (itStart == pointIds.end()) continue;
2158 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2159 if (itEnd == pointIds.end()) continue;
2160 if (abs(distance(itEnd, itStart)) != 1) continue;
2161 hitSegs.push_back(*it); // segment is undivided.
2164 return (pointIds.size() > 2); // something else apart start and end node
2167 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2168 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2170 using namespace std;
2171 int dst = distance(sIdxConn, sIdxConnE);
2172 modifiedFace.reserve(dst + insidePoints.size()-2);
2173 modifiedFace.resize(dst);
2174 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2176 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2177 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2178 if (startPos == shortEnd)
2179 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2180 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2181 if (endPos == shortEnd)
2182 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2183 int d = distance(startPos, endPos);
2184 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2185 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2187 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2194 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2195 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2196 * This method performs a conformization of \b this.
2198 * Only polyhedron cells are supported. You can call convertAllToPoly()
2200 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2201 * This method expects that all nodes in \a this are not closer than \a eps.
2202 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2204 * \param [in] eps the relative error to detect merged edges.
2205 * \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
2206 * that the user is expected to deal with.
2208 * \throw If \a this is not coherent.
2209 * \throw If \a this has not spaceDim equal to 3.
2210 * \throw If \a this has not meshDim equal to 3.
2211 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2213 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2215 using namespace std;
2217 static const int SPACEDIM=3;
2218 checkConsistencyLight();
2219 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2220 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2221 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2222 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2224 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2225 const double * coo(_coords->begin());
2226 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2229 /*************************
2231 *************************/
2232 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2233 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2234 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2235 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2236 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2239 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2240 const double *bbox(bboxArr->begin()); getCoords()->begin();
2241 int nDescCell(mDesc->getNumberOfCells());
2242 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2243 // Surfaces - handle biggest first
2244 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2245 DataArrayDouble * surfs = surfF->getArray();
2247 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2248 DataArrayDouble * normals = normalsF->getArray();
2249 const double * normalsP = normals->getConstPointer();
2251 // Sort faces by decreasing surface:
2252 vector< pair<double,int> > S;
2253 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2255 pair<double,int> p = make_pair(surfs->begin()[i], i);
2258 sort(S.rbegin(),S.rend()); // reverse sort
2259 vector<bool> hit(nDescCell);
2260 fill(hit.begin(), hit.end(), false);
2261 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2263 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2265 int faceIdx = (*it).second;
2266 if (hit[faceIdx]) continue;
2268 vector<int> candidates, cands2;
2269 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2270 // Keep only candidates whose normal matches the normal of current face
2271 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2273 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2274 if (*it2 != faceIdx && col)
2275 cands2.push_back(*it2);
2280 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2281 INTERP_KERNEL::TranslationRotationMatrix rotation;
2282 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2283 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2284 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2286 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2287 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2288 double * cooPartRef(mPartRef->_coords->getPointer());
2289 double * cooPartCand(mPartCand->_coords->getPointer());
2290 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2291 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2292 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2293 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2295 // Localize faces in 2D thanks to barycenters
2296 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2297 vector<int> compo; compo.push_back(2);
2298 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2299 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2300 if (!idsGoodPlane->getNumberOfTuples())
2303 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2305 compo[0] = 0; compo.push_back(1);
2306 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2307 mPartRef->changeSpaceDimension(2,0.0);
2308 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2309 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2311 if (!cc->getNumberOfTuples())
2313 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2316 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2317 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2320 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2321 throw INTERP_KERNEL::Exception(oss.str());
2325 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2327 if (!ids->getNumberOfTuples())
2330 double checkSurf=0.0;
2331 const int * idsGoodPlaneP(idsGoodPlane->begin());
2332 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2334 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2335 hit[faceIdx2] = true;
2336 checkSurf += surfs->begin()[faceIdx2];
2338 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2341 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2342 throw INTERP_KERNEL::Exception(oss.str());
2345 // For all polyhedrons using this face, replace connectivity:
2346 vector<int> polyIndices, packsIds, facePack;
2347 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2348 polyIndices.push_back(revDescP[ii]);
2349 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2351 // Current face connectivity
2352 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2353 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2354 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2355 // Deletion of old faces
2357 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2359 if (packsIds[jj] == -1)
2360 // The below should never happen - if a face is used several times, with a different layout of the nodes
2361 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2362 // faces which are actually used only once, by a single cell. This is different for edges below.
2363 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2365 connSla->deletePack(*it2, packsIds[jj]);
2367 // Insertion of new faces:
2368 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2370 // Build pack from the face to insert:
2371 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2373 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2374 // Insert it in all hit polyhedrons:
2375 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2376 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2381 // Set back modified connectivity
2382 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2383 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2384 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2387 /************************
2389 ************************/
2390 // Now we have a face-conform mesh.
2392 // Recompute descending
2393 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2394 // Rebuild desc connectivity with orientation this time!!
2395 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2396 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2397 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2398 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2399 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2400 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2401 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2402 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2403 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2404 // std::cout << "writing!\n";
2405 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2406 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2407 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2408 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2409 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2410 const double *bbox2(bboxArr->begin());
2411 int nDesc2Cell=mDesc2->getNumberOfCells();
2412 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2414 // Edges - handle longest first
2415 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2416 DataArrayDouble * lens = lenF->getArray();
2418 // Sort edges by decreasing length:
2419 vector<pair<double,int> > S;
2420 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2422 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2425 sort(S.rbegin(),S.rend()); // reverse sort
2427 vector<bool> hit(nDesc2Cell);
2428 fill(hit.begin(), hit.end(), false);
2430 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2432 int eIdx = (*it).second;
2436 vector<int> candidates, cands2;
2437 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2438 // Keep only candidates colinear with current edge
2440 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2441 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2442 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2443 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2446 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2447 for (int i3=0; i3 < 3; i3++)
2448 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2449 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2450 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2451 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2453 cands2.push_back(*it2);
2455 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2458 // Now rotate edges to bring them on Ox
2459 int startNode = cDesc2[cIDesc2[eIdx]+1];
2460 int endNode = cDesc2[cIDesc2[eIdx]+2];
2461 INTERP_KERNEL::TranslationRotationMatrix rotation;
2462 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2463 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2464 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2465 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2468 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2469 nbElemsNotM1 = tmp->getNbOfElems();
2471 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2472 double * cooPartRef(mPartRef->_coords->getPointer());
2473 double * cooPartCand(mPartCand->_coords->getPointer());
2474 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2475 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2476 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2477 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2480 // Eliminate all edges for which y or z is not null
2481 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2482 vector<int> compo; compo.push_back(1);
2483 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2485 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2486 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2487 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2488 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2489 if (!idsGoodLine->getNumberOfTuples())
2492 // Now the ordering along the Ox axis:
2493 std::vector<int> insidePoints, hitSegs;
2494 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2495 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2496 idsGoodLine->begin(), idsGoodLine->end(),
2497 /*out*/insidePoints, hitSegs);
2498 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2499 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2500 hit[cands2[*its]] = true;
2502 if (!isSplit) // current segment remains in one piece
2505 // Get original node IDs in global coords array
2506 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2507 *iit = nodeMapInv->begin()[*iit];
2509 vector<int> polyIndices, packsIds, facePack;
2510 // For each face implying this edge
2511 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2513 int faceIdx = revDescP2[ii];
2514 // each cell where this face is involved connectivity will be modified:
2515 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2517 // Current face connectivity
2518 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2519 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2521 vector<int> modifiedFace;
2522 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2523 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2524 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2528 // Rebuild 3D connectivity from descending:
2529 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2530 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2531 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2532 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2533 newConn->set3(superIdx, idx, vals);
2534 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2535 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2537 int sz, faceIdx = abs(descP[jj])-1;
2538 bool orient = descP[jj]>0;
2539 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2541 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2544 vector<int> rev(sz-1);
2545 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2546 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2550 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2553 ret = ret->buildUniqueNotSorted();