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::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMeshWithTree(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
275 std::map<INTERP_KERNEL::Node *,int>& mapp,
276 const BBTreePts<2,int> & nodeTree,
277 const std::map<int, INTERP_KERNEL::Node *>& mapRev)
280 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).
281 const double *coo=mDesc->getCoords()->getConstPointer();
282 const int *c=mDesc->getNodalConnectivity()->getConstPointer();
283 const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
285 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
286 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
287 for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
289 INTERP_KERNEL::Node *n;
290 // Look for a potential node to merge
291 std::vector<int> candNode;
292 nodeTree.getElementsAroundPoint(coo+2*(*it2), candNode);
293 if (candNode.size() > 2)
294 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MEDCouplingUMeshBuildQPFromMeshWithTree(): some nodes are not properly merged (within eps) in input mesh!");
295 bool node_created=false;
298 auto itt=mapRev.find(candNode[0]);
299 if (itt != mapRev.end()) // we might hit a node which is in the coords array but not used in the connectivity in which case it won't be in the revMap
307 n = new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
308 mapp2[*it2]=INTERP_KERNEL::NodeWithUsage(n,INTERP_KERNEL::USAGE_UNKNOWN);
310 INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
311 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
313 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
314 ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1)); // this call will set quad points to false in the map
316 for(std::map<int, INTERP_KERNEL::NodeWithUsage >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
318 if((*it2).second.second == INTERP_KERNEL::USAGE_LINEAR)
319 mapp[(*it2).second.first]=(*it2).first;
320 ((*it2).second.first)->decrRef();
325 INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
329 int locId=nodeId-offset2;
330 return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
334 int locId=nodeId-offset1;
335 return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
337 return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
341 * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
343 void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
344 const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
345 /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
347 for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
349 int eltId1=abs(*desc1)-1;
350 for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
352 std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
353 if(it==mappRev.end())
355 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
367 * 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 ) .
368 * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
369 * \param forbiddenPoints the list of points that should not be removed in the process
371 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset,
372 const std::map<int, bool>& forbiddenPoints,
373 DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
375 std::size_t sz(std::distance(connBg,connEnd));
376 if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
377 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
379 INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
380 INTERP_KERNEL::AutoPtr<int> tmpConn2(new int[sz]);
381 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
382 unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
383 unsigned nbOfHit(0); // number of fusions operated
384 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
385 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
386 INTERP_KERNEL::NormalizedCellType typeOfSon;
387 std::vector<int> middles;
389 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
391 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
392 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
393 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
394 posEndElt = posBaseElt+1;
396 // Look backward first: are the final edges of the cells colinear with the first ones?
397 // This initializes posBaseElt.
400 for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
402 cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn2,typeOfSon);
403 // Identify common point:
404 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
405 auto itE(forbiddenPoints.end());
406 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
408 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
409 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
410 bool isColinear=eint->areColinears();
421 // Update last connectivity
422 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
426 const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
427 for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell with one single edge
429 cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn2,typeOfSon); // get edge #j's connectivity
430 // Identify common point:
431 int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
432 auto itE(forbiddenPoints.end());
433 if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
435 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
436 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
437 bool isColinear(eint->areColinears());
448 // Update last connectivity
449 std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
451 //push [posBaseElt,posEndElt) in newConnOfCell using e
452 // 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!
454 // at the beginning of the connectivity (insert type)
455 EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
456 else if((nbOfHit+nbOfTurn) != (nbs-1))
458 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
459 if ((nbOfHit+nbOfTurn) == (nbs-1))
460 // at the end (only quad points to deal with)
461 EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
462 posBaseElt=posEndElt;
466 newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
472 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
474 if(candidates.empty())
476 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
478 const std::vector<int>& pool(intersectEdge1[*it]);
479 int tmp[2]; tmp[0]=start; tmp[1]=stop;
480 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
485 tmp[0]=stop; tmp[1]=start;
486 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
496 * This method performs the 2nd step of Partition of 2D mesh.
497 * This method has 4 inputs :
498 * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
499 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
500 * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
501 * 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'
502 * Nodes end up lying consecutively on a cutted edge.
503 * \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.
504 * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
505 * \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.
506 * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
507 * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
509 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
510 const std::vector<double>& addCoo,
511 const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
513 int offset1=m1->getNumberOfNodes();
514 int ncell2=m2->getNumberOfCells();
515 const int *c=m2->getNodalConnectivity()->begin();
516 const int *cI=m2->getNodalConnectivityIndex()->begin();
517 const double *coo=m2->getCoords()->begin();
518 const double *cooBis=m1->getCoords()->begin();
519 int offset2=offset1+m2->getNumberOfNodes();
520 intersectEdge.resize(ncell2);
521 for(int i=0;i<ncell2;i++,cI++)
523 const std::vector<int>& divs=subDiv[i];
524 int nnode=cI[1]-cI[0]-1;
525 std::map<int, INTERP_KERNEL::NodeWithUsage > mapp2;
526 std::map<INTERP_KERNEL::Node *, int> mapp22;
527 for(int j=0;j<nnode;j++)
529 INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
530 int nnid=c[(*cI)+j+1];
531 mapp2[nnid]=INTERP_KERNEL::NodeWithUsage(nn,INTERP_KERNEL::USAGE_UNKNOWN);
532 mapp22[nn]=nnid+offset1;
534 INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
535 for(std::map<int, INTERP_KERNEL::NodeWithUsage >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
536 ((*it).second.first)->decrRef();
537 std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
538 std::map<INTERP_KERNEL::Node *,int> mapp3;
539 for(std::size_t j=0;j<divs.size();j++)
542 INTERP_KERNEL::Node *tmp=0;
544 tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
546 tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
548 tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
552 e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
553 for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
559 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,
560 MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
562 idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
563 idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
564 int nCells(mesh1D->getNumberOfCells());
565 if(nCells!=(int)intersectEdge2.size())
566 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
567 const DataArrayDouble *coo2(mesh1D->getCoords());
568 const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
569 const double *coo2Ptr(coo2->begin());
570 int offset1(coords1->getNumberOfTuples());
571 int offset2(offset1+coo2->getNumberOfTuples());
572 int offset3(offset2+addCoo.size()/2);
573 std::vector<double> addCooQuad;
574 MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
575 int tmp[4],cicnt(0),kk(0);
576 for(int i=0;i<nCells;i++)
578 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
579 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
580 const std::vector<int>& subEdges(intersectEdge2[i]);
581 int nbSubEdge(subEdges.size()/2);
582 for(int j=0;j<nbSubEdge;j++,kk++)
584 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));
585 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
586 INTERP_KERNEL::Edge *e2Ptr(e2);
587 std::map<int,int>::const_iterator itm;
588 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
590 tmp[0]=INTERP_KERNEL::NORM_SEG3;
591 itm=mergedNodes.find(subEdges[2*j]);
592 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
593 itm=mergedNodes.find(subEdges[2*j+1]);
594 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
595 tmp[3]=offset3+(int)addCooQuad.size()/2;
597 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
599 cOut->insertAtTheEnd(tmp,tmp+4);
600 ciOut->pushBackSilent(cicnt);
604 tmp[0]=INTERP_KERNEL::NORM_SEG2;
605 itm=mergedNodes.find(subEdges[2*j]);
606 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
607 itm=mergedNodes.find(subEdges[2*j+1]);
608 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
610 cOut->insertAtTheEnd(tmp,tmp+3);
611 ciOut->pushBackSilent(cicnt);
614 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
616 idsInRetColinear->pushBackSilent(kk);
617 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
622 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
623 ret->setConnectivity(cOut,ciOut,true);
624 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
625 arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
626 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,(int)addCooQuad.size()/2,2);
627 std::vector<const DataArrayDouble *> coordss(4);
628 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
629 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
634 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
636 std::vector<int> allEdges;
637 for(const int *it2(descBg);it2!=descEnd;it2++)
639 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
641 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
643 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
645 std::size_t nb(allEdges.size());
647 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
648 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
649 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
650 ret->setCoords(coords);
651 ret->allocateCells(1);
652 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
653 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
654 connOut[kk]=allEdges[2*kk];
655 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
659 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
661 const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
662 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
664 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
665 if(sz!=std::distance(descBg,descEnd))
666 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
667 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
668 std::vector<int> allEdges,centers;
669 const double *coordsPtr(coords->begin());
670 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
671 int offset(coords->getNumberOfTuples());
672 for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
674 INTERP_KERNEL::NormalizedCellType typeOfSon;
675 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
676 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
678 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
680 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
682 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
684 {//the current edge has been subsplit -> create corresponding centers.
685 std::size_t nbOfCentersToAppend(edge1.size()/2);
686 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
687 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
688 std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
689 for(std::size_t k=0;k<nbOfCentersToAppend;k++)
692 const double *aa(coordsPtr+2*(*it3++));
693 const double *bb(coordsPtr+2*(*it3++));
694 ee->getMiddleOfPoints(aa,bb,tmpp);
695 addCoo->insertAtTheEnd(tmpp,tmpp+2);
696 centers.push_back(offset+k);
700 std::size_t nb(allEdges.size());
702 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
703 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
704 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
706 ret->setCoords(coords);
709 addCoo->rearrange(2);
710 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
711 ret->setCoords(addCoo);
713 ret->allocateCells(1);
714 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
715 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
716 connOut[kk]=allEdges[2*kk];
717 connOut.insert(connOut.end(),centers.begin(),centers.end());
718 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
723 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
726 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
728 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
730 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
731 if(!cm.isQuadratic())
732 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
734 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
737 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
740 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
742 const INTERP_KERNEL::Edge *ee(*it);
743 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
747 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
750 const double *coo(mesh2D->getCoords()->begin());
751 std::size_t sz(conn.size());
752 std::vector<double> addCoo;
753 std::vector<int> conn2(conn);
754 int offset(mesh2D->getNumberOfNodes());
755 for(std::size_t i=0;i<sz;i++)
758 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
759 addCoo.insert(addCoo.end(),tmp,tmp+2);
760 conn2.push_back(offset+(int)i);
762 mesh2D->getCoords()->rearrange(1);
763 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
764 mesh2D->getCoords()->rearrange(2);
765 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
770 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
772 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
773 * a set of edges defined in \a splitMesh1D.
775 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
776 std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
778 std::size_t nb(edge1Bis.size()/2);
779 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
780 int iEnd(splitMesh1D->getNumberOfCells());
782 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
784 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
785 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
786 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
789 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
790 out0.resize(1); out1.resize(1);
791 std::vector<int>& connOut(out0[0]);
792 connOut.resize(nbOfEdgesOf2DCellSplit);
793 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
794 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
795 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
797 connOut[kk]=edge1Bis[2*kk];
798 edgesPtr[kk]=edge1BisPtr[2*kk];
803 // [i,iEnd[ contains the
804 out0.resize(2); out1.resize(2);
805 std::vector<int>& connOutLeft(out0[0]);
806 std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
807 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
808 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
809 for(std::size_t k=ii;k<jj+1;k++)
810 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
811 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
812 for(int ik=0;ik<iEnd;ik++)
814 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
815 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
818 for(int ik=iEnd-1;ik>=0;ik--)
819 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
820 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
821 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
822 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
823 for(int ik=0;ik<iEnd;ik++)
824 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
825 eright.insert(eright.end(),ees.begin(),ees.end());
833 CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
835 std::vector<int> _edges;
836 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
839 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
841 std::size_t nbe(edges.size());
842 std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
843 for(std::size_t i=0;i<nbe;i++)
845 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
846 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
848 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
849 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
850 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
856 EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
857 EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
858 bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
859 void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
860 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
864 MCAuto<MEDCouplingUMesh> _mesh;
865 MCAuto<INTERP_KERNEL::Edge> _edge;
870 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
872 const MEDCouplingUMesh *mesh(_mesh);
878 { _left++; _right++; return ; }
879 if (_right > pos && _left != pos)
880 { _right++; return ; }
883 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
884 if((isLeft && isRight) || (!isLeft && !isRight))
885 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
896 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
897 if((isLeft && isRight) || (!isLeft && !isRight))
898 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
913 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
915 const MEDCouplingUMesh *mesh(_mesh);
918 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
921 {// not fully splitting cell case
922 if(mesh2D->getNumberOfCells()==1)
923 {//little optimization. 1 cell no need to find in which cell mesh is !
924 neighbors[0]=offset; neighbors[1]=offset;
929 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
930 int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
932 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
933 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
938 class VectorOfCellInfo
941 VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
942 std::size_t size() const { return _pool.size(); }
943 int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
944 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);
945 const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
946 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
947 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
948 void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
950 int getZePosOfEdgeGivenItsGlobalId(int pos) const;
951 void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
952 const CellInfo& get(int pos) const;
953 CellInfo& get(int pos);
955 std::vector<CellInfo> _pool;
956 MCAuto<MEDCouplingUMesh> _ze_mesh;
957 std::vector<EdgeInfo> _edge_info;
960 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
962 _pool[0]._edges=edges;
963 _pool[0]._edges_ptr=edgesPtr;
966 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
969 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
972 const MEDCouplingUMesh *zeMesh(_ze_mesh);
974 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
975 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
976 return zeMesh->getCellContainingPoint(barys->begin(),eps);
979 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
980 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
981 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
983 get(pos);//to check pos
984 bool isFast(pos==0 && _pool.size()==1);
985 std::size_t sz(edges.size());
986 // dealing with edges
988 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
990 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
992 std::vector<CellInfo> pool(_pool.size()-1+sz);
993 for(std::size_t i=0;i<pos;i++)
995 for(std::size_t j=0;j<sz;j++)
996 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
997 for(int i=pos+1;i<(int)_pool.size();i++)
998 pool[i+sz-1]=_pool[i];
1002 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
1010 std::vector< MCAuto<MEDCouplingUMesh> > ms;
1013 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
1017 if(pos<_ze_mesh->getNumberOfCells()-1)
1019 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
1022 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
1023 for(std::size_t j=0;j<ms2.size();j++)
1025 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
1028 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
1030 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
1033 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
1036 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
1038 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
1040 if((*it).isInMyRange(pos))
1043 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
1046 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
1048 get(pos);//to perform the sanity check;
1049 if(_edge_info.empty())
1051 std::size_t sz(_edge_info.size()-1);
1052 for(std::size_t i=0;i<sz;i++)
1053 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
1056 const CellInfo& VectorOfCellInfo::get(int pos) const
1058 if(pos<0 || pos>=(int)_pool.size())
1059 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1063 CellInfo& VectorOfCellInfo::get(int pos)
1065 if(pos<0 || pos>=(int)_pool.size())
1066 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1072 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1073 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1075 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1077 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1079 * \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.
1081 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1082 MCAuto<DataArrayInt>& idsLeftRight)
1084 int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1085 if(nbCellsInSplitMesh1D==0)
1086 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1087 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1088 std::size_t nb(allEdges.size()),jj;
1090 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1091 std::vector<int> edge1Bis(nb*2);
1092 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1093 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1094 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1095 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1096 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1098 idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1099 int *idsLeftRightPtr(idsLeftRight->getPointer());
1100 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1102 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1103 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1104 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1105 // of the connectivity.
1106 MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1107 renumb->alloc(nbCellsInSplitMesh1D,1);
1108 const int * renumbP(renumb->begin());
1110 int i, first=cSplitPtr[1];
1111 // Follow 1D line backward as long as it is connected:
1112 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1113 first=cSplitPtr[ciSplitPtr[i]+1];
1114 if (i < nbCellsInSplitMesh1D-1)
1116 // Build circular permutation to shift consecutive edges together
1118 renumb->applyModulus(nbCellsInSplitMesh1D);
1119 splitMesh1D->renumberCells(renumbP, false);
1120 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1121 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1127 for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1128 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1130 for(;iEnd<nbCellsInSplitMesh1D;)
1132 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1138 if(iEnd<nbCellsInSplitMesh1D)
1141 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1142 int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1144 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1145 retTmp->setCoords(splitMesh1D->getCoords());
1146 retTmp->allocateCells();
1148 std::vector< std::vector<int> > out0;
1149 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1151 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1152 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1153 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1154 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1158 for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1159 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1161 return pool.getZeMesh().retn();
1165 * splitMesh1D is an input parameter but might have its cells renumbered.
1167 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1168 const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1169 MCAuto<DataArrayInt>& idsLeftRight)
1171 const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1173 std::vector<int> allEdges;
1174 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1175 for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1177 int edgeId(std::abs(*it)-1);
1178 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1179 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1180 const std::vector<int>& edge1(intersectEdge1[edgeId]);
1182 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1184 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1185 std::size_t sz(edge1.size());
1186 for(std::size_t cnt=0;cnt<sz;cnt++)
1187 allEdgesPtr.push_back(ee);
1190 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1193 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1195 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1196 {//easy case comparison not
1197 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1199 else if(typ1.isQuadratic() && typ2.isQuadratic())
1201 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1204 if(conn1[2]==conn2[2])
1206 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1207 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1211 {//only one is quadratic
1212 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1215 const double *a(0),*bb(0),*be(0);
1216 if(typ1.isQuadratic())
1218 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1222 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1224 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1225 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1231 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1232 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1234 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1236 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1238 if(candidatesIn2DEnd==candidatesIn2DBg)
1239 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1240 const double *coo(mesh2DSplit->getCoords()->begin());
1241 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1242 return *candidatesIn2DBg;
1243 int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1244 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1245 if(cellIdInMesh1DSplitRelative<0)
1246 cur1D->changeOrientationOfCells();
1247 const int *c1D(cur1D->getNodalConnectivity()->begin());
1248 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1249 for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1251 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1252 const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1253 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1254 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1255 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1256 for(unsigned it2=0;it2<sz;it2++)
1258 INTERP_KERNEL::NormalizedCellType typeOfSon;
1259 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1260 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1261 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1265 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1269 * \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.
1270 * 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.
1271 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1272 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1273 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1274 * \param [out] addCoo - nodes to be append at the end
1275 * \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.
1277 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1278 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)
1280 static const int SPACEDIM=2;
1281 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1282 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1283 // Build BB tree of all edges in the tool mesh (second mesh)
1284 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1285 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1286 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1287 intersectEdge1.resize(nDescCell1);
1288 colinear2.resize(nDescCell2);
1289 subDiv2.resize(nDescCell2);
1290 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1291 BBTreePts<SPACEDIM,int> treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps);
1293 std::vector<int> candidates1(1);
1294 int offset1(m1Desc->getNumberOfNodes());
1295 int offset2(offset1+m2Desc->getNumberOfNodes());
1296 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1298 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1299 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1300 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1302 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1303 std::map<int, INTERP_KERNEL::Node *> revMap2;
1304 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1305 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1307 for (auto& kv : map2)
1308 revMap2[kv.second] = kv.first;
1310 // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged.
1311 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2);
1312 // 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
1313 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1314 std::set<INTERP_KERNEL::Node *> nodes;
1315 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1316 std::size_t szz(nodes.size());
1317 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1318 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1319 for(std::size_t iii=0;iii<szz;iii++,itt++)
1320 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1321 // end of protection
1322 // Performs edge cutting:
1323 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1328 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1329 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1335 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1336 * It builds the descending connectivity of the two meshes, and then using a binary tree
1337 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1338 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1340 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1341 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1342 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1343 std::vector<double>& addCoo,
1344 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1346 // Build desc connectivity
1347 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1348 desc2=DataArrayInt::New();
1349 descIndx2=DataArrayInt::New();
1350 revDesc2=DataArrayInt::New();
1351 revDescIndx2=DataArrayInt::New();
1352 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1353 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1354 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1355 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1356 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1357 std::map<int,int> notUsedMap;
1358 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1359 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1360 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1364 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1365 * (newly created) nodes corresponding to the edge intersections.
1367 * @param[out] cr, crI connectivity of the resulting mesh
1368 * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1369 * TODO: describe input parameters
1371 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1372 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1373 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1374 const std::vector<double>& addCoords,
1375 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1377 static const int SPACEDIM=2;
1378 const double *coo1(m1->getCoords()->begin());
1379 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1380 int offset1(m1->getNumberOfNodes());
1381 const double *coo2(m2->getCoords()->begin());
1382 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1383 int offset2(offset1+m2->getNumberOfNodes());
1384 int offset3(offset2+((int)addCoords.size())/2);
1385 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1386 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1387 // Here a BBTree on 2D-cells, not on segments:
1388 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1389 int ncell1(m1->getNumberOfCells());
1391 for(int i=0;i<ncell1;i++)
1393 std::vector<int> candidates2;
1394 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1395 std::map<INTERP_KERNEL::Node *,int> mapp;
1396 std::map<int,INTERP_KERNEL::Node *> mappRev;
1397 INTERP_KERNEL::QuadraticPolygon pol1;
1398 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1399 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1400 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1401 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1402 // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1403 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1404 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1406 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
1407 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1408 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1409 for(it1.first();!it1.finished();it1.next())
1410 edges1.insert(it1.current()->getPtr());
1412 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1413 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1415 // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1416 // Again all the additional intersecting nodes are there.
1417 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1419 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1420 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1421 // Complete mapping with elements coming from the current cell it2 in mesh2:
1422 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1423 // pol2 is the new QP in the final merged result.
1424 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1425 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1427 // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1428 for (auto &p: pol2s)
1429 p.cleanDegeneratedConsecutiveEdges();
1430 edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges.
1433 // Now rebuild intersected cells from all this:
1434 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1436 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1437 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1438 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1439 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1441 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1442 // by m2 but that we still want to keep in the final result.
1447 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1449 catch(INTERP_KERNEL::Exception& e)
1451 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();
1452 throw INTERP_KERNEL::Exception(oss.str());
1455 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1456 (*it).second->decrRef();
1460 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1462 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1465 std::size_t sz(conn.size());
1466 std::size_t found(std::numeric_limits<std::size_t>::max());
1467 for(std::size_t i=0;i<sz;i++)
1469 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1470 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]};
1471 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1472 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1473 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1475 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];
1476 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]);
1478 if(dotTest>eps && dotTest<1.-eps)
1484 if(found==std::numeric_limits<std::size_t>::max())
1485 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1486 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1489 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1491 std::size_t sz(conn.size());
1492 std::vector<int> *curPart(&part0);
1493 for(std::size_t i=0;i<sz;i++)
1495 int nextt(conn[(i+1)%sz]);
1496 (*curPart).push_back(nextt);
1497 if(nextt==pt0 || nextt==pt1)
1503 (*curPart).push_back(nextt);
1509 * 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.
1511 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1512 const int *desc, const int *descIndx, const double *coords, double eps,
1513 std::vector<std::vector<int> >& res) const
1515 checkFullyDefined();
1516 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1517 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1518 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1519 int nbOfCells(getNumberOfCells());
1521 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1522 for(int i=0;i<nbOfCells;i++)
1524 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1525 for(int j=0;j<nbOfFaces;j++)
1527 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1528 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1529 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1530 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1531 INTERP_KERNEL::NormalizedCellType cmsId;
1532 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1533 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1534 if(p.first!=-1 && p.second!=-1)
1538 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1539 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1540 std::vector<int> elt1,elt2;
1541 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1542 res.push_back(elt1);
1543 res.push_back(elt2);
1555 * 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).
1557 * \sa MEDCouplingUMesh::split2DCells
1559 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1561 checkConnectivityFullyDefined();
1562 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1563 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1564 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1565 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1566 int prevPosOfCi(ciPtr[0]);
1567 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1569 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1570 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1571 for(int j=0;j<sz;j++)
1573 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1574 for(int k=0;k<sz2;k++)
1575 *cPtr++=subPtr[offset2+k];
1577 *cPtr++=oldConn[prevPosOfCi+j+2];
1580 prevPosOfCi=ciPtr[1];
1581 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1584 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1585 _nodal_connec->decrRef();
1586 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1591 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1593 * \return int - the number of new nodes created.
1594 * \sa MEDCouplingUMesh::split2DCells
1596 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1598 checkConsistencyLight();
1599 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1600 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1601 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1602 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1603 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1604 const double *oldCoordsPtr(getCoords()->begin());
1605 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1606 int prevPosOfCi(ciPtr[0]);
1607 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1609 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1610 for(int j=0;j<sz;j++)
1611 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1612 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1613 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1615 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1619 cPtr[1]=oldConn[prevPosOfCi+2+j];
1620 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1623 std::vector<INTERP_KERNEL::Node *> ns(3);
1624 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1625 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1626 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1627 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1628 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1630 cPtr[1]=subPtr[offset2+k];
1631 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1633 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1636 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1638 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1639 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1642 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1643 _nodal_connec->decrRef();
1644 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1645 addCoo->rearrange(2);
1646 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1648 return addCoo->getNumberOfTuples();
1655 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1656 * returns a result mesh constituted by polygons.
1657 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1658 * all nodes from m2.
1659 * The meshes should be in 2D space. In
1660 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1662 * \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
1663 * 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)
1664 * \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
1665 * 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)
1666 * \param [in] eps - precision used to detect coincident mesh entities.
1667 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1668 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1669 * this array using decrRef() as it is no more needed.
1670 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1671 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1672 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1673 * any cell of \a m2. The caller is to delete this array using decrRef() as
1674 * it is no more needed.
1675 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1676 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1677 * is no more needed.
1678 * \throw If the coordinates array is not set in any of the meshes.
1679 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1680 * \throw If any of the meshes is not a 2D mesh in 2D space.
1682 * \sa conformize2D, mergeNodes
1684 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1685 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1688 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1689 m1->checkFullyDefined();
1690 m2->checkFullyDefined();
1691 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1692 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1693 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1695 // Step 1: compute all edge intersections (new nodes)
1696 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1697 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1698 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1699 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1700 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1701 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1702 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1703 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1704 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1705 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1707 // Step 2: re-order newly created nodes according to the ordering found in m2
1708 std::vector< std::vector<int> > intersectEdge2;
1709 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1710 subDiv2.clear(); dd5=0; dd6=0;
1713 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1714 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1715 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1716 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1718 // Step 4: Prepare final result:
1719 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1720 addCooDa->alloc((int)(addCoo.size())/2,2);
1721 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1722 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1723 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1724 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1725 std::vector<const DataArrayDouble *> coordss(4);
1726 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1727 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1728 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1729 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1730 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1731 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1732 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1733 ret->setConnectivity(conn,connI,true);
1734 ret->setCoords(coo);
1735 cellNb1=c1.retn(); cellNb2=c2.retn();
1740 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1741 * 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
1742 * and finally, in case of quadratic polygon the centers of edges new nodes.
1743 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1745 * \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
1746 * 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)
1747 * \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
1748 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1749 * \param [in] eps - precision used to perform intersections and localization operations.
1750 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1751 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1752 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1753 * 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.
1754 * \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
1755 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1756 * 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.
1758 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1760 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1762 if(!mesh2D || !mesh1D)
1763 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1764 mesh2D->checkFullyDefined();
1765 mesh1D->checkFullyDefined();
1766 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1767 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1768 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1769 // Step 1: compute all edge intersections (new nodes)
1770 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1771 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1772 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1774 // Build desc connectivity
1775 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1776 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1777 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1778 std::map<int,int> mergedNodes;
1779 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1780 // use mergeNodes to fix intersectEdge1
1781 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1783 std::size_t n((*it0).size()/2);
1784 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1785 std::map<int,int>::const_iterator it1;
1786 it1=mergedNodes.find(eltStart);
1787 if(it1!=mergedNodes.end())
1788 (*it0)[0]=(*it1).second;
1789 it1=mergedNodes.find(eltEnd);
1790 if(it1!=mergedNodes.end())
1791 (*it0)[2*n-1]=(*it1).second;
1794 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1795 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
1796 // Step 2: re-order newly created nodes according to the ordering found in m2
1797 std::vector< std::vector<int> > intersectEdge2;
1798 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1800 // Step 3: compute splitMesh1D
1801 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1802 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1803 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1804 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1805 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1806 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1807 // deal with cells in mesh2D that are not cut but only some of their edges are
1808 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1809 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1810 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1811 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
1812 if(!idsInDesc2DToBeRefined->empty())
1814 DataArrayInt *out0(0),*outi0(0);
1815 DataArrayInt::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1816 MCAuto<DataArrayInt> outi0s(outi0);
1818 out0s=out0s->buildUnique();
1822 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1823 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1824 MCAuto<DataArrayInt> elts,eltsIndex;
1825 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1826 MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1827 if (eltsIndex->getNumberOfTuples() > 1)
1828 eltsIndex2 = eltsIndex->deltaShiftIndex();
1829 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1830 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1831 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1832 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1833 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1834 if((DataArrayInt *)out0s)
1835 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1836 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1837 // OK all is ready to insert in ret2 mesh
1838 if(!untouchedCells->empty())
1839 {// the most easy part, cells in mesh2D not impacted at all
1840 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1841 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1842 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1844 if((DataArrayInt *)out0s)
1845 {// here dealing with cells in out0s but not in cellsToBeModified
1846 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1847 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1848 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1850 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1851 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1853 int offset(ret2->getNumberOfTuples());
1854 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1855 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1856 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1857 int kk(0),*ret3ptr(partOfRet3->getPointer());
1858 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1860 int faceId(std::abs(*it)-1);
1861 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1863 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1866 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1867 ret3ptr[2*kk]=tmp+offset;
1868 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1869 ret3ptr[2*kk+1]=tmp+offset;
1872 {//the current edge is shared by a 2D cell that will be split just after
1873 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1874 ret3ptr[2*kk]=-(*it2+1);
1875 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1876 ret3ptr[2*kk+1]=-(*it2+1);
1880 m1Desc->setCoords(ret1->getCoords());
1881 ret1NonCol->setCoords(ret1->getCoords());
1882 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1883 if(!outMesh2DSplit.empty())
1885 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1886 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1887 (*itt)->setCoords(da);
1890 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1891 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1893 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1894 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1895 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1896 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1897 MCAuto<DataArrayInt> partOfRet3;
1898 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));
1899 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1900 outMesh2DSplit.push_back(splitOfOneCell);
1901 for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1902 ret2->pushBackSilent(*it);
1905 std::size_t nbOfMeshes(outMesh2DSplit.size());
1906 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1907 for(std::size_t i=0;i<nbOfMeshes;i++)
1908 tmp[i]=outMesh2DSplit[i];
1910 ret1->getCoords()->setInfoOnComponents(compNames);
1911 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1912 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1914 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1915 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1917 int old2DCellId(-ret3->getIJ(*it,0)-1);
1918 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1919 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
1921 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1924 splitMesh1D=ret1.retn();
1925 splitMesh2D=ret2D.retn();
1926 cellIdInMesh2D=ret2.retn();
1927 cellIdInMesh1D=ret3.retn();
1931 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1932 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1933 * 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
1934 * 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).
1935 * 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.
1937 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1938 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1940 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1941 * This method expects that all nodes in \a this are not closer than \a eps.
1942 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1944 * \param [in] eps the relative error to detect merged edges.
1945 * \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
1946 * that the user is expected to deal with.
1948 * \throw If \a this is not coherent.
1949 * \throw If \a this has not spaceDim equal to 2.
1950 * \throw If \a this has not meshDim equal to 2.
1951 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1953 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1955 static const int SPACEDIM=2;
1956 checkConsistencyLight();
1957 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1958 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1959 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1960 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1961 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1962 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1963 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1964 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1965 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1966 std::vector<double> addCoo;
1967 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1968 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1969 for(int i=0;i<nDescCell;i++)
1971 std::vector<int> candidates;
1972 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1973 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1974 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1976 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1977 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1978 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1979 INTERP_KERNEL::MergePoints merge;
1980 INTERP_KERNEL::QuadraticPolygon c1,c2;
1981 e1->intersectWith(e2,merge,c1,c2);
1982 e1->decrRef(); e2->decrRef();
1983 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1984 overlapEdge[i].push_back(*it);
1985 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1986 overlapEdge[*it].push_back(i);
1989 // splitting done. sort intersect point in intersectEdge.
1990 std::vector< std::vector<int> > middle(nDescCell);
1991 int nbOf2DCellsToBeSplit(0);
1992 bool middleNeedsToBeUsed(false);
1993 std::vector<bool> cells2DToTreat(nDescCell,false);
1994 for(int i=0;i<nDescCell;i++)
1996 std::vector<int>& isect(intersectEdge[i]);
1997 int sz((int)isect.size());
2000 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
2001 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
2002 e->sortSubNodesAbs(coords,isect);
2007 int idx0(rdi[i]),idx1(rdi[i+1]);
2009 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
2010 if(!cells2DToTreat[rd[idx0]])
2012 cells2DToTreat[rd[idx0]]=true;
2013 nbOf2DCellsToBeSplit++;
2015 // try to reuse at most eventual 'middle' of SEG3
2016 std::vector<int>& mid(middle[i]);
2017 mid.resize(sz+1,-1);
2018 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
2020 middleNeedsToBeUsed=true;
2021 const std::vector<int>& candidates(overlapEdge[i]);
2022 std::vector<int> trueCandidates;
2023 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
2024 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
2025 trueCandidates.push_back(*itc);
2026 int stNode(c[ci[i]+1]),endNode(isect[0]);
2027 for(int j=0;j<sz+1;j++)
2029 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
2031 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
2032 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
2033 { mid[j]=*itc; break; }
2036 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
2041 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
2042 if(nbOf2DCellsToBeSplit==0)
2045 int *retPtr(ret->getPointer());
2046 for(int i=0;i<nCell;i++)
2047 if(cells2DToTreat[i])
2050 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
2051 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
2052 DataArrayInt::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
2053 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
2054 if(middleNeedsToBeUsed)
2055 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
2056 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
2057 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
2058 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.
2059 setPartOfMySelf(ret->begin(),ret->end(),*modif);
2061 bool areNodesMerged; int newNbOfNodes;
2062 if(nbOfNodesCreated!=0)
2063 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2069 * 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.
2070 * 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).
2071 * 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
2072 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2073 * 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
2074 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2076 * 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
2077 * using new instance, idem for coordinates.
2079 * 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.
2081 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
2083 * \throw If \a this is not coherent.
2084 * \throw If \a this has not spaceDim equal to 2.
2085 * \throw If \a this has not meshDim equal to 2.
2087 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2089 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2091 return internalColinearize2D(eps, false);
2095 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2096 * 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
2097 * merged, contrary to colinearize2D().
2099 * \sa MEDCouplingUMesh::colinearize2D
2101 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2103 return internalColinearize2D(eps, true);
2108 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2110 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2112 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2113 checkConsistencyLight();
2114 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2115 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2116 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2117 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2118 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2119 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2120 std::map<int, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2123 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2124 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2125 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2126 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2127 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2128 MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2129 MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2130 const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2131 for(auto it = ids->begin(); it != ids->end(); it++)
2132 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2135 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2136 const double *coords(_coords->begin());
2137 int *newciptr(newci->getPointer());
2138 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2140 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2141 ret->pushBackSilent(i);
2142 newciptr[1]=newc->getNumberOfTuples();
2147 if(!appendedCoords->empty())
2149 appendedCoords->rearrange(2);
2150 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2152 setCoords(newCoords);
2155 setConnectivity(newc,newci,true);
2161 * c, cI describe a wire mesh in 3D space, in local numbering
2162 * startNode, endNode in global numbering
2163 *\return true if the segment is indeed split
2165 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2166 const int * c, const int * cI, const int *idsBg, const int *endBg,
2167 std::vector<int> & pointIds, std::vector<int> & hitSegs)
2169 using namespace std;
2171 const int SPACEDIM=3;
2172 typedef pair<double, int> PairDI;
2174 for (const int * it = idsBg; it != endBg; ++it)
2176 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2177 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2178 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2179 x.insert(make_pair(coo[end*SPACEDIM], end));
2182 vector<PairDI> xx(x.begin(), x.end());
2183 sort(xx.begin(),xx.end());
2184 pointIds.reserve(xx.size());
2186 // Keep what is inside [startNode, endNode]:
2188 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2190 const int idx = (*it).second;
2193 if (idx == startNode) go = 1;
2194 if (idx == endNode) go = 2;
2195 if (go) pointIds.push_back(idx);
2198 pointIds.push_back(idx);
2199 if (idx == endNode || idx == startNode)
2203 // vector<int> pointIds2(pointIds.size()+2);
2204 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2205 // pointIds2[0] = startNode;
2206 // pointIds2[pointIds2.size()-1] = endNode;
2209 reverse(pointIds.begin(), pointIds.end());
2211 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2212 for (const int * it = idsBg; it != endBg; ++it)
2214 int start = c[cI[*it]+1], end = c[cI[*it]+2];
2215 vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2216 if (itStart == pointIds.end()) continue;
2217 vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2218 if (itEnd == pointIds.end()) continue;
2219 if (abs(distance(itEnd, itStart)) != 1) continue;
2220 hitSegs.push_back(*it); // segment is undivided.
2223 return (pointIds.size() > 2); // something else apart start and end node
2226 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2227 const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2229 using namespace std;
2230 int dst = distance(sIdxConn, sIdxConnE);
2231 modifiedFace.reserve(dst + insidePoints.size()-2);
2232 modifiedFace.resize(dst);
2233 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2235 vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2236 vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2237 if (startPos == shortEnd)
2238 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2239 vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2240 if (endPos == shortEnd)
2241 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2242 int d = distance(startPos, endPos);
2243 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2244 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2246 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2253 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2254 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2255 * This method performs a conformization of \b this.
2257 * Only polyhedron cells are supported. You can call convertAllToPoly()
2259 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2260 * This method expects that all nodes in \a this are not closer than \a eps.
2261 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2263 * \param [in] eps the relative error to detect merged edges.
2264 * \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
2265 * that the user is expected to deal with.
2267 * \throw If \a this is not coherent.
2268 * \throw If \a this has not spaceDim equal to 3.
2269 * \throw If \a this has not meshDim equal to 3.
2270 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2272 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2274 using namespace std;
2276 static const int SPACEDIM=3;
2277 checkConsistencyLight();
2278 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2279 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2280 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2281 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2283 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2284 const double * coo(_coords->begin());
2285 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2288 /*************************
2290 *************************/
2291 MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2292 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2293 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2294 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2295 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2298 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2299 const double *bbox(bboxArr->begin()); getCoords()->begin();
2300 int nDescCell(mDesc->getNumberOfCells());
2301 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2302 // Surfaces - handle biggest first
2303 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2304 DataArrayDouble * surfs = surfF->getArray();
2306 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2307 DataArrayDouble * normals = normalsF->getArray();
2308 const double * normalsP = normals->getConstPointer();
2310 // Sort faces by decreasing surface:
2311 vector< pair<double,int> > S;
2312 for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2314 pair<double,int> p = make_pair(surfs->begin()[i], i);
2317 sort(S.rbegin(),S.rend()); // reverse sort
2318 vector<bool> hit(nDescCell);
2319 fill(hit.begin(), hit.end(), false);
2320 vector<int> hitPoly; // the final result: which 3D cells have been modified.
2322 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2324 int faceIdx = (*it).second;
2325 if (hit[faceIdx]) continue;
2327 vector<int> candidates, cands2;
2328 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2329 // Keep only candidates whose normal matches the normal of current face
2330 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2332 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2333 if (*it2 != faceIdx && col)
2334 cands2.push_back(*it2);
2339 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2340 INTERP_KERNEL::TranslationRotationMatrix rotation;
2341 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2342 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2343 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2345 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2346 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2347 double * cooPartRef(mPartRef->_coords->getPointer());
2348 double * cooPartCand(mPartCand->_coords->getPointer());
2349 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2350 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2351 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2352 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2354 // Localize faces in 2D thanks to barycenters
2355 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2356 vector<int> compo; compo.push_back(2);
2357 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2358 MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2359 if (!idsGoodPlane->getNumberOfTuples())
2362 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2364 compo[0] = 0; compo.push_back(1);
2365 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2366 mPartRef->changeSpaceDimension(2,0.0);
2367 MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2368 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2370 if (!cc->getNumberOfTuples())
2372 MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2375 MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2376 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2379 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2380 throw INTERP_KERNEL::Exception(oss.str());
2384 MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2386 if (!ids->getNumberOfTuples())
2389 double checkSurf=0.0;
2390 const int * idsGoodPlaneP(idsGoodPlane->begin());
2391 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2393 int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2394 hit[faceIdx2] = true;
2395 checkSurf += surfs->begin()[faceIdx2];
2397 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2400 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2401 throw INTERP_KERNEL::Exception(oss.str());
2404 // For all polyhedrons using this face, replace connectivity:
2405 vector<int> polyIndices, packsIds, facePack;
2406 for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2407 polyIndices.push_back(revDescP[ii]);
2408 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2410 // Current face connectivity
2411 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2412 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2413 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2414 // Deletion of old faces
2416 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2418 if (packsIds[jj] == -1)
2419 // The below should never happen - if a face is used several times, with a different layout of the nodes
2420 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2421 // faces which are actually used only once, by a single cell. This is different for edges below.
2422 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2424 connSla->deletePack(*it2, packsIds[jj]);
2426 // Insertion of new faces:
2427 for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2429 // Build pack from the face to insert:
2430 int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2432 const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2433 // Insert it in all hit polyhedrons:
2434 for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2435 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2440 // Set back modified connectivity
2441 MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2442 MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2443 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2446 /************************
2448 ************************/
2449 // Now we have a face-conform mesh.
2451 // Recompute descending
2452 MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2453 // Rebuild desc connectivity with orientation this time!!
2454 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2455 const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2456 const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2457 const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2458 MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2459 MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2460 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2461 MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2462 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2463 // std::cout << "writing!\n";
2464 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2465 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2466 const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2467 const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2468 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2469 const double *bbox2(bboxArr->begin());
2470 int nDesc2Cell=mDesc2->getNumberOfCells();
2471 BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2473 // Edges - handle longest first
2474 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2475 DataArrayDouble * lens = lenF->getArray();
2477 // Sort edges by decreasing length:
2478 vector<pair<double,int> > S;
2479 for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2481 pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2484 sort(S.rbegin(),S.rend()); // reverse sort
2486 vector<bool> hit(nDesc2Cell);
2487 fill(hit.begin(), hit.end(), false);
2489 for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2491 int eIdx = (*it).second;
2495 vector<int> candidates, cands2;
2496 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2497 // Keep only candidates colinear with current edge
2499 unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2500 for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2501 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2502 for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2505 unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2506 for (int i3=0; i3 < 3; i3++)
2507 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2508 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2509 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2510 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2512 cands2.push_back(*it2);
2514 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2517 // Now rotate edges to bring them on Ox
2518 int startNode = cDesc2[cIDesc2[eIdx]+1];
2519 int endNode = cDesc2[cIDesc2[eIdx]+2];
2520 INTERP_KERNEL::TranslationRotationMatrix rotation;
2521 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2522 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2523 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2524 MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2527 MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2528 nbElemsNotM1 = tmp->getNbOfElems();
2530 MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2531 double * cooPartRef(mPartRef->_coords->getPointer());
2532 double * cooPartCand(mPartCand->_coords->getPointer());
2533 for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2534 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2535 for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2536 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2539 // Eliminate all edges for which y or z is not null
2540 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2541 vector<int> compo; compo.push_back(1);
2542 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2544 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2545 MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2546 MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2547 MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2548 if (!idsGoodLine->getNumberOfTuples())
2551 // Now the ordering along the Ox axis:
2552 std::vector<int> insidePoints, hitSegs;
2553 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2554 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2555 idsGoodLine->begin(), idsGoodLine->end(),
2556 /*out*/insidePoints, hitSegs);
2557 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2558 for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2559 hit[cands2[*its]] = true;
2561 if (!isSplit) // current segment remains in one piece
2564 // Get original node IDs in global coords array
2565 for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2566 *iit = nodeMapInv->begin()[*iit];
2568 vector<int> polyIndices, packsIds, facePack;
2569 // For each face implying this edge
2570 for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2572 int faceIdx = revDescP2[ii];
2573 // each cell where this face is involved connectivity will be modified:
2574 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2576 // Current face connectivity
2577 const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2578 const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2580 vector<int> modifiedFace;
2581 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2582 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2583 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2587 // Rebuild 3D connectivity from descending:
2588 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2589 MCAuto<DataArrayInt> superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2590 MCAuto<DataArrayInt> idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0);
2591 MCAuto<DataArrayInt> vals(DataArrayInt::New()); vals->alloc(0);
2592 newConn->set3(superIdx, idx, vals);
2593 for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2594 for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2596 int sz, faceIdx = abs(descP[jj])-1;
2597 bool orient = descP[jj]>0;
2598 const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2600 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2603 vector<int> rev(sz-1);
2604 for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2605 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2609 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2612 ret = ret->buildUniqueNotSorted();