1 // Copyright (C) 2007-2021 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 mcIdType InternalAddPoint(const INTERP_KERNEL::Edge *e, mcIdType id, const double *coo, mcIdType startId, mcIdType endId, DataArrayDouble& addCoo, mcIdType& nodesCnter)
59 mcIdType ret(nodesCnter++);
61 e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
62 addCoo.insertAtTheEnd(newPt,newPt+2);
67 mcIdType InternalAddPointOriented(const INTERP_KERNEL::Edge *e, mcIdType id, const double *coo, mcIdType startId, mcIdType endId, DataArrayDouble& addCoo, mcIdType& nodesCnter)
73 mcIdType 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 mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<mcIdType>& middles)
85 int trueStart(start>=0?start:nbOfEdges+start);
86 tmp[0]=ToIdType(linOrArc?INTERP_KERNEL::NORM_QPOLYG:INTERP_KERNEL::NORM_POLYGON); tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
87 newConnOfCell->insertAtTheEnd(tmp,tmp+3);
92 mcIdType 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 mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<mcIdType>& middles)
103 mcIdType tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
104 newConnOfCell->pushBackSilent(tmpEnd);
109 mcIdType 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 mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<mcIdType>& 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 mcIdType tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]); // % to handle last seg.
126 mcIdType 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>,mcIdType>& m, mcIdType forbVal0, mcIdType forbVal1, std::vector<mcIdType>& isect)
137 MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
138 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType>::const_iterator it(m.find(nTmp));
140 throw INTERP_KERNEL::Exception("Internal error in remapping !");
141 mcIdType v((*it).second);
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>,mcIdType>& m, mcIdType forbVal0, mcIdType forbVal1, std::vector<mcIdType>& 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 mcIdType *bg, const double *coords2D, std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType>& 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<mcIdType, INTERP_KERNEL::NodeWithUsage >& mapp2, const mcIdType *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<mcIdType>& candidates,
244 std::map<INTERP_KERNEL::Node *,mcIdType>& mapp)
247 std::map<mcIdType, 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 mcIdType *c=mDesc->getNodalConnectivity()->getConstPointer();
250 const mcIdType *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
251 std::set<mcIdType> s;
252 for(std::vector<mcIdType>::const_iterator it=candidates.begin();it!=candidates.end();it++)
253 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
254 for(std::set<mcIdType>::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<mcIdType>::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<mcIdType, 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<mcIdType>& candidates,
275 std::map<INTERP_KERNEL::Node *,mcIdType>& mapp,
276 const BBTreePts<2,mcIdType> & nodeTree,
277 const std::map<mcIdType, INTERP_KERNEL::Node *>& mapRev)
280 std::map<mcIdType, 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 mcIdType *c=mDesc->getNodalConnectivity()->getConstPointer();
283 const mcIdType *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
284 std::set<mcIdType> s;
285 for(std::vector<mcIdType>::const_iterator it=candidates.begin();it!=candidates.end();it++)
286 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
287 for(std::set<mcIdType>::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<mcIdType> 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<mcIdType>::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<mcIdType, 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(mcIdType nodeId, const double *coo1, mcIdType offset1, const double *coo2, mcIdType offset2, const std::vector<double>& addCoo)
329 mcIdType locId=nodeId-offset2;
330 return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
334 mcIdType 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, mcIdType offset1, const double *coo2, mcIdType offset2, const std::vector<double>& addCoo,
344 const mcIdType *desc1Bg, const mcIdType *desc1End, const std::vector<std::vector<mcIdType> >& intesctEdges1,
345 /*output*/std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, std::map<mcIdType,INTERP_KERNEL::Node *>& mappRev)
347 for(const mcIdType *desc1=desc1Bg;desc1!=desc1End;desc1++)
349 mcIdType eltId1=std::abs(*desc1)-1;
350 for(std::vector<mcIdType>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
352 std::map<mcIdType,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 mcIdType *connBg, const mcIdType *connEnd, mcIdType offset,
372 const std::map<mcIdType, bool>& forbiddenPoints,
373 DataArrayIdType *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<mcIdType> tmpConn(new mcIdType[sz]);
380 INTERP_KERNEL::AutoPtr<mcIdType> tmpConn2(new mcIdType[sz]);
381 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
382 unsigned nbs(cm.getNumberOfSons2(connBg+1,ToIdType(sz)));
383 unsigned nbOfHit(0); // number of fusions operated
384 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
385 const std::size_t 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<mcIdType> middles;
389 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
391 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,ToIdType(sz),tmpConn,typeOfSon);
392 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> 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,ToIdType(sz),tmpConn2,typeOfSon);
403 // Identify common point:
404 mcIdType commPoint = std::find((mcIdType *)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((mcIdType *)tmpConn2, tmpConn2+sz, (mcIdType *)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(j,connBg+1,ToIdType(sz),tmpConn2,typeOfSon); // get edge #j's connectivity
430 // Identify common point:
431 mcIdType commPoint = std::find((mcIdType *)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((mcIdType *)tmpConn2, tmpConn2+sz, (mcIdType *)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,nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
456 else if((nbOfHit+nbOfTurn) != (nbs-1))
458 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,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,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<mcIdType> >& intersectEdge1, const std::vector<mcIdType>& candidates, mcIdType start, mcIdType stop, mcIdType& retVal)
474 if(candidates.empty())
476 for(std::vector<mcIdType>::const_iterator it=candidates.begin();it!=candidates.end();it++)
478 const std::vector<mcIdType>& pool(intersectEdge1[*it]);
479 mcIdType 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<mcIdType> >& subDiv, std::vector< std::vector<mcIdType> >& intersectEdge)
513 mcIdType offset1=m1->getNumberOfNodes();
514 mcIdType ncell2=m2->getNumberOfCells();
515 const mcIdType *c=m2->getNodalConnectivity()->begin();
516 const mcIdType *cI=m2->getNodalConnectivityIndex()->begin();
517 const double *coo=m2->getCoords()->begin();
518 const double *cooBis=m1->getCoords()->begin();
519 mcIdType offset2=offset1+m2->getNumberOfNodes();
520 intersectEdge.resize(ncell2);
521 for(mcIdType i=0;i<ncell2;i++,cI++)
523 const std::vector<mcIdType>& divs=subDiv[i];
524 mcIdType nnode=cI[1]-cI[0]-1;
525 std::map<mcIdType, INTERP_KERNEL::NodeWithUsage > mapp2;
526 std::map<INTERP_KERNEL::Node *, mcIdType> mapp22;
527 for(mcIdType 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 mcIdType 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<mcIdType, 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 *,mcIdType> 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++)
560 * Build the final 1D mesh resulting from the newly created points after intersection with the segments of the descending 2D mesh.
561 * @param[out] idsInRetColinear IDs of edges in the result (ret) that are colinears to one of the segment of the descending 2D mesh. Indexing scheme
562 * is the one of the ret 1D mesh.
563 * @param[out] idsInMesh1DForIdsInRetColinear same IDs as above in the descending 2D mesh
565 MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<mcIdType> >& intersectEdge2,
566 const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<mcIdType,mcIdType>& mergedNodes,
567 const std::vector< std::vector<mcIdType> >& colinear2, const std::vector< std::vector<mcIdType> >& intersectEdge1,
568 MCAuto<DataArrayIdType>& idsInRetColinear, MCAuto<DataArrayIdType>& idsInMesh1DForIdsInRetColinear)
570 idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1);
571 idsInMesh1DForIdsInRetColinear=DataArrayIdType::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
572 mcIdType nCells=mesh1D->getNumberOfCells();
573 if(nCells!=ToIdType(intersectEdge2.size()))
574 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
575 const DataArrayDouble *coo2(mesh1D->getCoords());
576 const mcIdType *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
577 const double *coo2Ptr(coo2->begin());
578 mcIdType offset1(coords1->getNumberOfTuples());
579 mcIdType offset2(offset1+coo2->getNumberOfTuples());
580 mcIdType offset3(offset2+ToIdType(addCoo.size())/2);
581 std::vector<double> addCooQuad;
582 MCAuto<DataArrayIdType> cOut(DataArrayIdType::New()),ciOut(DataArrayIdType::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
583 mcIdType tmp[4],cicnt(0),kk(0);
584 for(mcIdType i=0;i<nCells;i++)
586 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
587 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
588 const std::vector<mcIdType>& subEdges(intersectEdge2[i]);
589 mcIdType nbSubEdge=ToIdType(subEdges.size()/2);
590 for(mcIdType j=0;j<nbSubEdge;j++,kk++)
592 MCAuto<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),
593 n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
594 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
595 INTERP_KERNEL::Edge *e2Ptr(e2);
596 std::map<mcIdType,mcIdType>::const_iterator itm;
597 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
599 tmp[0]=INTERP_KERNEL::NORM_SEG3;
600 itm=mergedNodes.find(subEdges[2*j]);
601 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
602 itm=mergedNodes.find(subEdges[2*j+1]);
603 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
604 tmp[3]=offset3+ToIdType(addCooQuad.size()/2);
606 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
608 cOut->insertAtTheEnd(tmp,tmp+4);
609 ciOut->pushBackSilent(cicnt);
613 tmp[0]=INTERP_KERNEL::NORM_SEG2;
614 itm=mergedNodes.find(subEdges[2*j]);
615 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
616 itm=mergedNodes.find(subEdges[2*j+1]);
617 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
619 cOut->insertAtTheEnd(tmp,tmp+3);
620 ciOut->pushBackSilent(cicnt);
623 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
625 idsInRetColinear->pushBackSilent(kk);
626 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
632 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
633 ret->setConnectivity(cOut,ciOut,true);
634 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
635 arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2);
636 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,addCooQuad.size()/2,2);
637 std::vector<const DataArrayDouble *> coordss(4);
638 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
639 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
644 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
646 std::vector<mcIdType> allEdges;
647 for(const mcIdType *it2(descBg);it2!=descEnd;it2++)
649 const std::vector<mcIdType>& edge1(intersectEdge1[std::abs(*it2)-1]);
651 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
653 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
655 std::size_t nb(allEdges.size());
657 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
658 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
659 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
660 ret->setCoords(coords);
661 ret->allocateCells(1);
662 std::vector<mcIdType> connOut(nbOfEdgesOf2DCellSplit);
663 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
664 connOut[kk]=allEdges[2*kk];
665 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(connOut.size()),&connOut[0]);
669 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
671 const mcIdType *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
672 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
674 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
675 if(sz!=std::distance(descBg,descEnd))
676 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
677 INTERP_KERNEL::AutoPtr<mcIdType> tmpPtr(new mcIdType[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
678 std::vector<mcIdType> allEdges,centers;
679 const double *coordsPtr(coords->begin());
680 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
681 mcIdType offset(coords->getNumberOfTuples());
682 for(const mcIdType *it2(descBg);it2!=descEnd;it2++,ii++)
684 INTERP_KERNEL::NormalizedCellType typeOfSon;
685 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
686 const std::vector<mcIdType>& edge1(intersectEdge1[std::abs(*it2)-1]);
688 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
690 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
692 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
694 {//the current edge has been subsplit -> create corresponding centers.
695 mcIdType nbOfCentersToAppend=ToIdType(edge1.size()/2);
696 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
697 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
698 std::vector<mcIdType>::const_iterator it3(allEdges.end()-edge1.size());
699 for(mcIdType k=0;k<nbOfCentersToAppend;k++)
702 const double *aa(coordsPtr+2*(*it3++));
703 const double *bb(coordsPtr+2*(*it3++));
704 ee->getMiddleOfPoints(aa,bb,tmpp);
705 addCoo->insertAtTheEnd(tmpp,tmpp+2);
706 centers.push_back(offset+k);
710 std::size_t nb(allEdges.size());
712 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
713 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
714 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
716 ret->setCoords(coords);
719 addCoo->rearrange(2);
720 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
721 ret->setCoords(addCoo);
723 ret->allocateCells(1);
724 std::vector<mcIdType> connOut(nbOfEdgesOf2DCellSplit);
725 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
726 connOut[kk]=allEdges[2*kk];
727 connOut.insert(connOut.end(),centers.begin(),centers.end());
728 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(connOut.size()),&connOut[0]);
733 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
736 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
738 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
740 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
741 if(!cm.isQuadratic())
742 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
744 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
747 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<mcIdType>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
750 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
752 const INTERP_KERNEL::Edge *ee(*it);
753 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
757 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(conn.size()),&conn[0]);
760 const double *coo(mesh2D->getCoords()->begin());
761 std::size_t sz(conn.size());
762 std::vector<double> addCoo;
763 std::vector<mcIdType> conn2(conn);
764 mcIdType offset(mesh2D->getNumberOfNodes());
765 for(std::size_t i=0;i<sz;i++)
768 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
769 addCoo.insert(addCoo.end(),tmp,tmp+2);
770 conn2.push_back(offset+ToIdType(i));
772 mesh2D->getCoords()->rearrange(1);
773 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
774 mesh2D->getCoords()->rearrange(2);
775 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(conn2.size()),&conn2[0]);
780 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
782 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
783 * a set of edges defined in \a splitMesh1D.
785 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<mcIdType>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
786 std::vector< std::vector<mcIdType> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
788 std::size_t nb(edge1Bis.size()/2);
789 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
790 mcIdType iEnd=splitMesh1D->getNumberOfCells();
792 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
794 const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
795 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
796 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
799 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
800 out0.resize(1); out1.resize(1);
801 std::vector<mcIdType>& connOut(out0[0]);
802 connOut.resize(nbOfEdgesOf2DCellSplit);
803 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
804 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
805 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
807 connOut[kk]=edge1Bis[2*kk];
808 edgesPtr[kk]=edge1BisPtr[2*kk];
813 // [i,iEnd[ contains the
814 out0.resize(2); out1.resize(2);
815 std::vector<mcIdType>& connOutLeft(out0[0]);
816 std::vector<mcIdType>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
817 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
818 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
819 for(std::size_t k=ii;k<jj+1;k++)
820 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
821 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
822 for(mcIdType ik=0;ik<iEnd;ik++)
824 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
825 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
828 for(mcIdType ik=iEnd-1;ik>=0;ik--)
829 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
830 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
831 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
832 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
833 for(mcIdType ik=0;ik<iEnd;ik++)
834 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
835 eright.insert(eright.end(),ees.begin(),ees.end());
843 CellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
845 std::vector<mcIdType> _edges;
846 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
849 CellInfo::CellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
851 std::size_t nbe(edges.size());
852 std::vector<mcIdType> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
853 for(std::size_t i=0;i<nbe;i++)
855 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
856 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
858 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
859 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
860 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
866 EdgeInfo(mcIdType istart, mcIdType iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
867 EdgeInfo(mcIdType istart, mcIdType iend, mcIdType pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
868 bool isInMyRange(mcIdType pos) const { return pos>=_istart && pos<_iend; }
869 void somethingHappendAt(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
870 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, mcIdType offset, mcIdType neighbors[2]) const;
874 MCAuto<MEDCouplingUMesh> _mesh;
875 MCAuto<INTERP_KERNEL::Edge> _edge;
876 mcIdType _left; // index (local numbering) of the left 2D cell bordering the edge '_edge'
877 mcIdType _right; // same as above, right side.
881 * Update indices of left and right 2D cell bordering the current edge.
883 void EdgeInfo::somethingHappendAt(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
885 const MEDCouplingUMesh *mesh(_mesh);
891 { _left++; _right++; return ; }
892 if (_right > pos && _left != pos)
893 { _right++; return ; }
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 # 1 !");
909 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
910 if((isLeft && isRight) || (!isLeft && !isRight))
911 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
926 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, mcIdType offset, mcIdType neighbors[2]) const
928 const MEDCouplingUMesh *mesh(_mesh);
931 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
934 {// not fully splitting cell case
935 if(mesh2D->getNumberOfCells()==1)
936 {//little optimization. 1 cell no need to find in which cell mesh is !
937 neighbors[0]=offset; neighbors[1]=offset;
942 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
943 mcIdType cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
945 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
946 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
951 class VectorOfCellInfo
954 VectorOfCellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
955 std::size_t size() const { return _pool.size(); }
956 mcIdType getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
957 void setMeshAt(mcIdType pos, const MCAuto<MEDCouplingUMesh>& mesh, mcIdType istart, mcIdType iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<mcIdType> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs);
958 const std::vector<mcIdType>& getConnOf(mcIdType pos) const { return get(pos)._edges; }
959 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(mcIdType pos) const { return get(pos)._edges_ptr; }
960 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
961 void feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const;
963 mcIdType getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const;
964 void updateEdgeInfo(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
965 const CellInfo& get(mcIdType pos) const;
966 CellInfo& get(mcIdType pos);
968 std::vector<CellInfo> _pool; // for a newly created 2D cell, the list of edges ToIdType( and edges ptr constiuing it
969 MCAuto<MEDCouplingUMesh> _ze_mesh; // the aggregated mesh
970 std::vector<EdgeInfo> _edge_info; // for each new edge added when cuting the 2D cell, the information on left and right bordering 2D cell
973 VectorOfCellInfo::VectorOfCellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
975 _pool[0]._edges=edges;
976 _pool[0]._edges_ptr=edgesPtr;
979 mcIdType VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
982 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
985 const MEDCouplingUMesh *zeMesh(_ze_mesh);
987 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
988 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
989 return zeMesh->getCellContainingPoint(barys->begin(),eps);
992 void VectorOfCellInfo::setMeshAt(mcIdType pos, const MCAuto<MEDCouplingUMesh>& mesh, mcIdType istart, mcIdType iend,
993 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<mcIdType> >& edges,
994 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
996 get(pos);//to check pos
997 bool isFast(pos==0 && _pool.size()==1);
998 std::size_t sz(edges.size());
999 // dealing with edges
1001 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
1003 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
1005 std::vector<CellInfo> pool(_pool.size()-1+sz);
1006 for(mcIdType i=0;i<pos;i++)
1008 for(std::size_t j=0;j<sz;j++)
1009 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
1010 for(std::size_t i=pos+1;i<_pool.size();i++)
1011 pool[i+sz-1]=_pool[i];
1015 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
1023 std::vector< MCAuto<MEDCouplingUMesh> > ms;
1026 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
1030 if(pos<_ze_mesh->getNumberOfCells()-1)
1032 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
1035 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
1036 for(std::size_t j=0;j<ms2.size();j++)
1038 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
1041 void VectorOfCellInfo::feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const
1043 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
1046 mcIdType VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const
1049 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
1051 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
1053 if((*it).isInMyRange(pos))
1056 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
1059 void VectorOfCellInfo::updateEdgeInfo(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
1061 get(pos);//to perform the sanity check;
1062 if(_edge_info.empty())
1064 std::size_t sz(_edge_info.size()-1);
1065 for(std::size_t i=0;i<sz;i++)
1066 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
1069 const CellInfo& VectorOfCellInfo::get(mcIdType pos) const
1071 if(pos<0 || pos>=ToIdType(_pool.size()))
1072 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1076 CellInfo& VectorOfCellInfo::get(mcIdType pos)
1078 if(pos<0 || pos>=ToIdType(_pool.size()))
1079 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1085 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1086 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1088 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1090 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1092 * \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.
1094 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<mcIdType>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, mcIdType offset,
1095 MCAuto<DataArrayIdType>& idsLeftRight)
1097 mcIdType nbCellsInSplitMesh1D=splitMesh1D->getNumberOfCells();
1098 if(nbCellsInSplitMesh1D==0)
1099 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1100 const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1101 std::size_t nb(allEdges.size()),jj;
1103 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1104 std::vector<mcIdType> edge1Bis(nb*2);
1105 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1106 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1107 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1108 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1109 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1111 idsLeftRight=DataArrayIdType::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1112 mcIdType *idsLeftRightPtr(idsLeftRight->getPointer());
1113 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1115 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1116 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1117 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1118 // of the connectivity.
1119 MCAuto <DataArrayIdType> renumb(DataArrayIdType::New());
1120 renumb->alloc(nbCellsInSplitMesh1D,1);
1121 const mcIdType * renumbP(renumb->begin());
1123 mcIdType i, first=cSplitPtr[1];
1124 // Follow 1D line backward as long as it is connected:
1125 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1126 first=cSplitPtr[ciSplitPtr[i]+1];
1127 if (i < nbCellsInSplitMesh1D-1)
1129 // Build circular permutation to shift consecutive edges together
1131 renumb->applyModulus(nbCellsInSplitMesh1D);
1132 splitMesh1D->renumberCells(renumbP, false);
1133 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1134 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1139 // The 1D first piece is used to intersect the 2D cell resulting in max two 2D cells.
1140 // The next 1D piece is localised (getPositionOf()) into this previous cut.
1141 // The result of the next intersection replaces the former single 2D cell that has been cut in the
1142 // pool. The neighbourhood information detained by pool._edge_info is also updated so that left and right
1143 // adjacent 2D cell of a 1D piece is kept up to date.
1144 // And so on and so forth.
1145 for(mcIdType iStart=0;iStart<nbCellsInSplitMesh1D;)
1146 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1147 mcIdType iEnd(iStart);
1148 for(;iEnd<nbCellsInSplitMesh1D;)
1150 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1156 if(iEnd<nbCellsInSplitMesh1D)
1159 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1160 mcIdType pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1162 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1163 retTmp->setCoords(splitMesh1D->getCoords());
1164 retTmp->allocateCells();
1166 std::vector< std::vector<mcIdType> > out0;
1167 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1169 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1170 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1171 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1172 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1176 for(mcIdType mm=0;mm<nbCellsInSplitMesh1D;mm++)
1177 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1179 return pool.getZeMesh().retn();
1183 * splitMesh1D is an input parameter but might have its cells renumbered.
1185 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, mcIdType cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1186 const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1, mcIdType offset,
1187 MCAuto<DataArrayIdType>& idsLeftRight)
1189 const mcIdType *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1191 std::vector<mcIdType> allEdges;
1192 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1193 for(const mcIdType *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1195 mcIdType edgeId(std::abs(*it)-1);
1196 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
1197 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1198 const std::vector<mcIdType>& edge1(intersectEdge1[edgeId]);
1200 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1202 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1203 std::size_t sz(edge1.size());
1204 for(std::size_t cnt=0;cnt<sz;cnt++)
1205 allEdgesPtr.push_back(ee);
1208 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1211 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const mcIdType *conn1, const INTERP_KERNEL::CellModel& typ2, const mcIdType *conn2, double eps)
1213 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1214 {//easy case comparison not
1215 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1217 else if(typ1.isQuadratic() && typ2.isQuadratic())
1219 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1222 if(conn1[2]==conn2[2])
1224 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1225 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1229 {//only one is quadratic
1230 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1233 const double *a(0),*bb(0),*be(0);
1234 if(typ1.isQuadratic())
1236 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1240 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1242 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1243 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1249 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1250 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1252 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1254 mcIdType FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const mcIdType *candidatesIn2DBg, const mcIdType *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, mcIdType cellIdInMesh1DSplitRelative, double eps)
1256 if(candidatesIn2DEnd==candidatesIn2DBg)
1257 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1258 const double *coo(mesh2DSplit->getCoords()->begin());
1259 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1260 return *candidatesIn2DBg;
1261 mcIdType edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1262 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1263 if(cellIdInMesh1DSplitRelative<0)
1264 cur1D->changeOrientationOfCells();
1265 const mcIdType *c1D(cur1D->getNodalConnectivity()->begin());
1266 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1267 for(const mcIdType *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1269 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1270 const mcIdType *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1271 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1272 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1273 INTERP_KERNEL::AutoPtr<mcIdType> tmpPtr(new mcIdType[ci[1]-ci[0]]);
1274 for(unsigned it2=0;it2<sz;it2++)
1276 INTERP_KERNEL::NormalizedCellType typeOfSon;
1277 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1278 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1279 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1283 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1287 * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of mcIdType given resp start and stop.
1288 * 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.
1289 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1290 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1291 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1292 * \param [out] addCoo - nodes to be append at the end
1293 * \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.
1295 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1296 std::vector< std::vector<mcIdType> >& intersectEdge1, std::vector< std::vector<mcIdType> >& colinear2, std::vector< std::vector<mcIdType> >& subDiv2, std::vector<double>& addCoo, std::map<mcIdType,mcIdType>& mergedNodes)
1298 static const int SPACEDIM=2;
1299 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1300 const mcIdType *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1301 // Build BB tree of all edges in the tool mesh (second mesh)
1302 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1303 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1304 mcIdType nDescCell1=m1Desc->getNumberOfCells(),nDescCell2=m2Desc->getNumberOfCells();
1305 intersectEdge1.resize(nDescCell1);
1306 colinear2.resize(nDescCell2);
1307 subDiv2.resize(nDescCell2);
1308 BBTree<SPACEDIM,mcIdType> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1309 BBTreePts<SPACEDIM,mcIdType> treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps);
1311 std::vector<mcIdType> candidates1(1);
1312 mcIdType offset1(m1Desc->getNumberOfNodes());
1313 mcIdType offset2(offset1+m2Desc->getNumberOfNodes());
1314 for(mcIdType i=0;i<nDescCell1;i++) // for all edges in the first mesh
1316 std::vector<mcIdType> candidates2; // edges of mesh2 candidate for intersection
1317 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1318 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1320 std::map<INTERP_KERNEL::Node *,mcIdType> map1,map2;
1321 std::map<mcIdType, INTERP_KERNEL::Node *> revMap2;
1322 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1323 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1325 for (auto& kv : map2)
1326 revMap2[kv.second] = kv.first;
1328 // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged.
1329 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2);
1330 // 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
1331 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1332 std::set<INTERP_KERNEL::Node *> nodes;
1333 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1334 std::size_t szz(nodes.size());
1335 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1336 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1337 for(std::size_t iii=0;iii<szz;iii++,itt++)
1338 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1339 // end of protection
1340 // Performs edge cutting:
1341 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1346 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1347 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1353 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1354 * It builds the descending connectivity of the two meshes, and then using a binary tree
1355 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1356 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1358 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1359 std::vector< std::vector<mcIdType> >& intersectEdge1, std::vector< std::vector<mcIdType> >& colinear2, std::vector< std::vector<mcIdType> >& subDiv2,
1360 MEDCouplingUMesh *& m1Desc, DataArrayIdType *&desc1, DataArrayIdType *&descIndx1, DataArrayIdType *&revDesc1, DataArrayIdType *&revDescIndx1,
1361 std::vector<double>& addCoo,
1362 MEDCouplingUMesh *& m2Desc, DataArrayIdType *&desc2, DataArrayIdType *&descIndx2, DataArrayIdType *&revDesc2, DataArrayIdType *&revDescIndx2)
1364 // Build desc connectivity
1365 desc1=DataArrayIdType::New(); descIndx1=DataArrayIdType::New(); revDesc1=DataArrayIdType::New(); revDescIndx1=DataArrayIdType::New();
1366 desc2=DataArrayIdType::New();
1367 descIndx2=DataArrayIdType::New();
1368 revDesc2=DataArrayIdType::New();
1369 revDescIndx2=DataArrayIdType::New();
1370 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1371 MCAuto<DataArrayIdType> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1372 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1373 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1374 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1375 std::map<mcIdType,mcIdType> notUsedMap;
1376 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1377 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1378 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1382 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1383 * (newly created) nodes corresponding to the edge intersections.
1385 * @param[out] cr connectivity of the resulting mesh
1386 * @param[out] crI connectivity of the resulting mesh
1387 * @param[out] cNb1 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1388 * @param[out] cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1389 * TODO: describe input parameters
1391 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const mcIdType *desc1, const mcIdType *descIndx1,
1392 const std::vector<std::vector<mcIdType> >& intesctEdges1, const std::vector< std::vector<mcIdType> >& colinear2,
1393 const MEDCouplingUMesh *m2, const mcIdType *desc2, const mcIdType *descIndx2, const std::vector<std::vector<mcIdType> >& intesctEdges2,
1394 const std::vector<double>& addCoords,
1395 std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& cr, std::vector<mcIdType>& crI, std::vector<mcIdType>& cNb1, std::vector<mcIdType>& cNb2)
1397 static const int SPACEDIM=2;
1398 const double *coo1(m1->getCoords()->begin());
1399 const mcIdType *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1400 mcIdType offset1(m1->getNumberOfNodes());
1401 const double *coo2(m2->getCoords()->begin());
1402 const mcIdType *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1403 mcIdType offset2(offset1+m2->getNumberOfNodes());
1404 mcIdType offset3(offset2+ToIdType(addCoords.size())/2);
1405 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1406 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1407 // Here a BBTree on 2D-cells, not on segments:
1408 BBTree<SPACEDIM,mcIdType> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1409 mcIdType ncell1=m1->getNumberOfCells();
1411 for(mcIdType i=0;i<ncell1;i++)
1413 std::vector<mcIdType> candidates2;
1414 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1415 std::map<INTERP_KERNEL::Node *,mcIdType> mapp;
1416 std::map<mcIdType,INTERP_KERNEL::Node *> mappRev;
1417 INTERP_KERNEL::QuadraticPolygon pol1;
1418 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1419 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1420 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1421 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1422 // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1423 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1424 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1426 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
1427 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1428 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1429 for(it1.first();!it1.finished();it1.next())
1430 edges1.insert(it1.current()->getPtr());
1432 std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1433 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1435 // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1436 // Again all the additional intersecting nodes are there.
1437 for(std::vector<mcIdType>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1439 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1440 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1441 // Complete mapping with elements coming from the current cell it2 in mesh2:
1442 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1443 // pol2 is the new QP in the final merged result.
1444 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1445 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1447 // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1448 for (auto &p: pol2s)
1449 p.cleanDegeneratedConsecutiveEdges();
1450 edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges.
1453 // Now rebuild intersected cells from all this:
1454 for(std::vector<mcIdType>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1456 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1457 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1458 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1459 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1461 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1462 // by m2 but that we still want to keep in the final result.
1467 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1469 catch(INTERP_KERNEL::Exception& e)
1471 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();
1472 throw INTERP_KERNEL::Exception(oss.str());
1475 for(std::map<mcIdType,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1476 (*it).second->decrRef();
1480 void InsertNodeInConnIfNecessary(mcIdType nodeIdToInsert, std::vector<mcIdType>& conn, const double *coords, double eps)
1482 std::vector<mcIdType>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1485 std::size_t sz(conn.size());
1486 std::size_t found(std::numeric_limits<std::size_t>::max());
1487 for(std::size_t i=0;i<sz;i++)
1489 mcIdType pt0(conn[i]),pt1(conn[(i+1)%sz]);
1490 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]};
1491 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1492 std::transform(v1,v1+3,v1,std::bind(std::multiplies<double>(),std::placeholders::_1,1./normm));
1493 std::transform(v2,v2+3,v2,std::bind(std::multiplies<double>(),std::placeholders::_1,1./normm));
1495 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];
1496 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]);
1498 if(dotTest>eps && dotTest<1.-eps)
1504 if(found==std::numeric_limits<std::size_t>::max())
1505 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1506 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1509 void SplitIntoToPart(const std::vector<mcIdType>& conn, mcIdType pt0, mcIdType pt1, std::vector<mcIdType>& part0, std::vector<mcIdType>& part1)
1511 std::size_t sz(conn.size());
1512 std::vector<mcIdType> *curPart(&part0);
1513 for(std::size_t i=0;i<sz;i++)
1515 mcIdType nextt(conn[(i+1)%sz]);
1516 (*curPart).push_back(nextt);
1517 if(nextt==pt0 || nextt==pt1)
1523 (*curPart).push_back(nextt);
1529 * 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.
1531 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1532 const mcIdType *desc, const mcIdType *descIndx, const double *coords, double eps,
1533 std::vector<std::vector<mcIdType> >& res) const
1535 checkFullyDefined();
1536 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1537 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1538 const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1539 mcIdType nbOfCells=getNumberOfCells();
1541 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1542 for(mcIdType i=0;i<nbOfCells;i++)
1544 mcIdType offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1545 for(int j=0;j<nbOfFaces;j++)
1547 const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1548 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1549 mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1550 INTERP_KERNEL::AutoPtr<mcIdType> tmp(new mcIdType[sz]);
1551 INTERP_KERNEL::NormalizedCellType cmsId;
1552 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1553 std::vector<mcIdType> elt((mcIdType *)tmp,(mcIdType *)tmp+nbOfNodesSon);
1554 if(p.first!=-1 && p.second!=-1)
1558 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1559 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1560 std::vector<mcIdType> elt1,elt2;
1561 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1562 res.push_back(elt1);
1563 res.push_back(elt2);
1575 * 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).
1577 * \sa MEDCouplingUMesh::split2DCells
1579 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI)
1581 checkConnectivityFullyDefined();
1582 mcIdType ncells=getNumberOfCells();
1583 mcIdType lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1584 MCAuto<DataArrayIdType> c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach);
1585 const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1586 mcIdType *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1587 mcIdType prevPosOfCi(ciPtr[0]);
1588 for(mcIdType i=0;i<ncells;i++,ciPtr++,descIPtr++)
1590 mcIdType offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1591 *cPtr++=ToIdType(INTERP_KERNEL::NORM_POLYGON); *cPtr++=oldConn[prevPosOfCi+1];
1592 for(mcIdType j=0;j<sz;j++)
1594 mcIdType offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1595 for(mcIdType k=0;k<sz2;k++)
1596 *cPtr++=subPtr[offset2+k];
1598 *cPtr++=oldConn[prevPosOfCi+j+2];
1601 prevPosOfCi=ciPtr[1];
1602 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1605 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1606 _nodal_connec->decrRef();
1607 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1612 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1614 * \return mcIdType - the number of new nodes created.
1615 * \sa MEDCouplingUMesh::split2DCells
1617 mcIdType MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI, const DataArrayIdType *mid, const DataArrayIdType *midI)
1619 checkConsistencyLight();
1620 mcIdType ncells=getNumberOfCells();
1621 mcIdType lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples());
1622 mcIdType nodesCnt(getNumberOfNodes());
1623 MCAuto<DataArrayIdType> c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach);
1624 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1625 const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1626 const mcIdType *midPtr(mid->begin()),*midIPtr(midI->begin());
1627 const double *oldCoordsPtr(getCoords()->begin());
1628 mcIdType *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1629 mcIdType prevPosOfCi(ciPtr[0]);
1630 for(mcIdType i=0;i<ncells;i++,ciPtr++,descIPtr++)
1632 mcIdType offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1633 for(mcIdType j=0;j<sz;j++)
1634 { mcIdType sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1635 *cPtr++=ToIdType(INTERP_KERNEL::NORM_QPOLYG); cPtr[0]=oldConn[prevPosOfCi+1];
1636 for(mcIdType j=0;j<sz;j++)//loop over subedges of oldConn
1638 mcIdType offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1642 cPtr[1]=oldConn[prevPosOfCi+2+j];
1643 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1646 std::vector<INTERP_KERNEL::Node *> ns(3);
1647 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1648 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1649 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1650 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1651 for(mcIdType k=0;k<sz2;k++)//loop over subsplit of current subedge
1653 cPtr[1]=subPtr[offset2+k];
1654 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1656 mcIdType tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1659 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1661 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1662 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1665 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1666 _nodal_connec->decrRef();
1667 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1668 addCoo->rearrange(2);
1669 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1671 return addCoo->getNumberOfTuples();
1678 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1679 * returns a result mesh constituted by polygons.
1680 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1681 * all nodes from m2.
1682 * The meshes should be in 2D space. In
1683 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1685 * \b WARNING: the two meshes should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and
1686 * orientCorrectly2DCells() can be used for this.
1687 * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells)
1688 * \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
1689 * 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)
1690 * \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
1691 * 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)
1692 * \param [in] eps - precision used to detect coincident mesh entities.
1693 * \param [out] cellNb1 - a new instance of DataArrayIdType holding for each result
1694 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1695 * this array using decrRef() as it is no more needed.
1696 * \param [out] cellNb2 - a new instance of DataArrayIdType holding for each result
1697 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1698 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1699 * any cell of \a m2. The caller is to delete this array using decrRef() as
1700 * it is no more needed.
1701 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1702 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1703 * is no more needed.
1704 * \throw If the coordinates array is not set in any of the meshes.
1705 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1706 * \throw If any of the meshes is not a 2D mesh in 2D space.
1708 * \sa conformize2D, mergeNodes
1710 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1711 double eps, DataArrayIdType *&cellNb1, DataArrayIdType *&cellNb2)
1714 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1715 m1->checkFullyDefined();
1716 m2->checkFullyDefined();
1717 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1718 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1719 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1721 // Step 1: compute all edge intersections (new nodes)
1722 std::vector< std::vector<mcIdType> > intersectEdge1, colinear2, subDiv2;
1723 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1724 DataArrayIdType *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1725 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1726 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1727 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1728 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1729 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1730 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1731 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1733 // Step 2: re-order newly created nodes according to the ordering found in m2
1734 std::vector< std::vector<mcIdType> > intersectEdge2;
1735 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1736 subDiv2.clear(); dd5=0; dd6=0;
1739 std::vector<mcIdType> cr,crI; //no DataArrayIdType because interface with Geometric2D
1740 std::vector<mcIdType> cNb1,cNb2; //no DataArrayIdType because interface with Geometric2D
1741 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1742 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1744 // Step 4: Prepare final result:
1745 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1746 addCooDa->alloc(addCoo.size()/2,2);
1747 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1748 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1749 addCoordsQuadraticDa->alloc(addCoordsQuadratic.size()/2,2);
1750 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1751 std::vector<const DataArrayDouble *> coordss(4);
1752 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1753 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1754 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1755 MCAuto<DataArrayIdType> conn(DataArrayIdType::New()); conn->alloc(cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1756 MCAuto<DataArrayIdType> connI(DataArrayIdType::New()); connI->alloc(crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1757 MCAuto<DataArrayIdType> c1(DataArrayIdType::New()); c1->alloc(cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1758 MCAuto<DataArrayIdType> c2(DataArrayIdType::New()); c2->alloc(cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1759 ret->setConnectivity(conn,connI,true);
1760 ret->setCoords(coo);
1761 cellNb1=c1.retn(); cellNb2=c2.retn();
1766 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1767 * 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
1768 * and finally, in case of quadratic polygon the centers of edges new nodes.
1769 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1771 * \b WARNING: the 2D mesh should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and
1772 * orientCorrectly2DCells() can be used for this.
1773 * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells)
1775 * \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
1776 * 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)
1777 * \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
1778 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1779 * \param [in] eps - precision used to perform intersections and localization operations.
1780 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1781 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1782 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1783 * 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.
1784 * \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
1785 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1786 * 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.
1788 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1790 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayIdType *&cellIdInMesh2D, DataArrayIdType *&cellIdInMesh1D)
1792 if(!mesh2D || !mesh1D)
1793 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1794 mesh2D->checkFullyDefined();
1795 mesh1D->checkFullyDefined();
1796 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1797 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1798 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1799 // Step 1: compute all edge intersections (new nodes)
1800 std::vector< std::vector<mcIdType> > intersectEdge1, colinear2, subDiv2;
1801 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1802 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1804 // Build desc connectivity
1805 DataArrayIdType *desc1(DataArrayIdType::New()),*descIndx1(DataArrayIdType::New()),*revDesc1(DataArrayIdType::New()),*revDescIndx1(DataArrayIdType::New());
1806 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1807 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1808 std::map<mcIdType,mcIdType> mergedNodes;
1809 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1810 // use mergeNodes to fix intersectEdge1
1811 for(std::vector< std::vector<mcIdType> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1813 std::size_t n((*it0).size()/2);
1814 mcIdType eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1815 std::map<mcIdType,mcIdType>::const_iterator it1;
1816 it1=mergedNodes.find(eltStart);
1817 if(it1!=mergedNodes.end())
1818 (*it0)[0]=(*it1).second;
1819 it1=mergedNodes.find(eltEnd);
1820 if(it1!=mergedNodes.end())
1821 (*it0)[2*n-1]=(*it1).second;
1824 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1825 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2);
1826 // Step 2: re-order newly created nodes according to the ordering found in m2
1827 std::vector< std::vector<mcIdType> > intersectEdge2;
1828 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1830 // Step 3: compute splitMesh1D
1831 MCAuto<DataArrayIdType> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1832 MCAuto<DataArrayIdType> ret2(DataArrayIdType::New()); ret2->alloc(0,1);
1833 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1834 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1836 // ### Colinearity fix :
1837 // if a node in ret1 has been merged with an already existing node in mesh2D,
1838 // we might end up with edges in ret1 that are colinear with some edges in mesh2D
1839 // even if none of the edges of the two original meshes were colinear.
1840 // this procedure detects such edges and adds them in idsInRet1Colinear/idsInDescMesh2DForIdsInRetColinear
1841 // a- find edges in ret1 that are in contact with 2 cells
1842 MCAuto<DataArrayDouble> centerOfMassRet1(ret1->computeCellCenterOfMass());
1843 MCAuto<DataArrayIdType> cells,cellsIndex;
1844 mesh2D->getCellsContainingPoints(centerOfMassRet1->begin(),centerOfMassRet1->getNumberOfTuples(),eps,cells,cellsIndex);
1845 MCAuto<DataArrayIdType> cellsIndex2(DataArrayIdType::New()); cellsIndex2->alloc(0,1);
1846 if (cellsIndex->getNumberOfTuples() > 1)
1847 cellsIndex2 = cellsIndex->deltaShiftIndex();
1848 MCAuto<DataArrayIdType> idsInRet1With2Contacts(cellsIndex2->findIdsEqual(2));
1850 MCAuto<DataArrayIdType> realIdsInDesc2D(desc1->deepCopy());
1851 realIdsInDesc2D->abs(); realIdsInDesc2D->applyLin(1,-1);
1852 const mcIdType *cRet1(ret1->getNodalConnectivity()->begin()),*ciRet1(ret1->getNodalConnectivityIndex()->begin());
1853 for(const mcIdType *it=idsInRet1With2Contacts->begin();it!=idsInRet1With2Contacts->end();it++)
1855 // b- find the edge that the 2 cells in m1Desc have in common:
1856 // this is the edge which is colinear with the one in ret1
1857 const mcIdType* cellId1 = cells->begin() + cellsIndex->begin()[*it];
1858 const mcIdType* cellId2 = cells->begin() + cellsIndex->begin()[*it]+1;
1860 std::set<mcIdType> s1(realIdsInDesc2D->begin()+dd2->begin()[*cellId1], realIdsInDesc2D->begin()+dd2->begin()[*cellId1+1]);
1861 std::set<mcIdType> s2(realIdsInDesc2D->begin()+dd2->begin()[*cellId2], realIdsInDesc2D->begin()+dd2->begin()[*cellId2+1]);
1863 std::vector<mcIdType> commonEdgeId;
1864 std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(), std::back_inserter(commonEdgeId));
1866 // c- find correct orientation for commonEdgeId
1867 const mcIdType* firstNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+1;
1868 const mcIdType* secondNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+2;
1869 mcIdType commonEdgeIdCorrectlyOriented;
1870 if(IsColinearOfACellOf(intersectEdge1, commonEdgeId, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1, commonEdgeIdCorrectlyOriented))
1872 idsInRet1Colinear->pushBackSilent(*it);
1873 idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeIdCorrectlyOriented);
1876 // ### End colinearity fix
1878 MCAuto<DataArrayIdType> ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<mcIdType>::max()); ret3->rearrange(2);
1879 MCAuto<DataArrayIdType> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1880 // deal with cells in mesh2D that are not cut but only some of their edges are
1881 MCAuto<DataArrayIdType> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1882 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1883 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1885 MCAuto<DataArrayIdType> out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells
1886 if(!idsInDesc2DToBeRefined->empty())
1888 DataArrayIdType *out0(0),*outi0(0);
1889 DataArrayIdType::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1890 MCAuto<DataArrayIdType> outi0s(outi0);
1892 out0s=out0s->buildUnique();
1896 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1897 MCAuto<DataArrayDouble> baryRet1(centerOfMassRet1->selectByTupleId(idsInRet1NotColinear->begin(), idsInRet1NotColinear->end()));
1898 DataArrayIdType *out(0),*outi(0);
1899 DataArrayIdType::ExtractFromIndexedArrays(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end(),cells,cellsIndex,out,outi);
1900 MCAuto<DataArrayIdType> elts(out),eltsIndex(outi);
1902 MCAuto<DataArrayIdType> eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1);
1903 if (eltsIndex->getNumberOfTuples() > 1)
1904 eltsIndex2 = eltsIndex->deltaShiftIndex();
1905 MCAuto<DataArrayIdType> eltsIndex3(eltsIndex2->findIdsEqual(1));
1906 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1907 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1908 MCAuto<DataArrayIdType> cellsToBeModified(elts->buildUnique());
1909 MCAuto<DataArrayIdType> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1910 if((DataArrayIdType *)out0s)
1911 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1912 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1913 // OK all is ready to insert in ret2 mesh
1914 if(!untouchedCells->empty())
1915 {// the most easy part, cells in mesh2D not impacted at all
1916 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1917 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1918 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1920 if((DataArrayIdType *)out0s)
1921 {// here dealing with cells in out0s but not in cellsToBeModified
1922 MCAuto<DataArrayIdType> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1923 const mcIdType *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1924 for(const mcIdType *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1926 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1927 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1929 mcIdType offset(ret2->getNumberOfTuples());
1930 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1931 MCAuto<DataArrayIdType> partOfRet3(DataArrayIdType::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1932 partOfRet3->fillWithValue(std::numeric_limits<mcIdType>::max()); partOfRet3->rearrange(2);
1933 mcIdType kk(0),*ret3ptr(partOfRet3->getPointer());
1934 for(const mcIdType *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1936 mcIdType faceId(std::abs(*it)-1);
1937 for(const mcIdType *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1939 mcIdType tmp(fewModifiedCells->findIdFirstEqual(*it2));
1942 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1943 ret3ptr[2*kk]=tmp+offset;
1944 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1945 ret3ptr[2*kk+1]=tmp+offset;
1948 {//the current edge is shared by a 2D cell that will be split just after
1949 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1950 ret3ptr[2*kk]=-(*it2+1);
1951 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1952 ret3ptr[2*kk+1]=-(*it2+1);
1956 m1Desc->setCoords(ret1->getCoords());
1957 ret1NonCol->setCoords(ret1->getCoords());
1958 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1959 if(!outMesh2DSplit.empty())
1961 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1962 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1963 (*itt)->setCoords(da);
1966 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1967 for(const mcIdType *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1969 MCAuto<DataArrayIdType> idsNonColPerCell(elts->findIdsEqual(*it));
1970 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1971 MCAuto<DataArrayIdType> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1972 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1973 MCAuto<DataArrayIdType> partOfRet3;
1974 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));
1975 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1976 outMesh2DSplit.push_back(splitOfOneCell);
1977 for(mcIdType i=0;i<splitOfOneCell->getNumberOfCells();i++)
1978 ret2->pushBackSilent(*it);
1981 std::size_t nbOfMeshes(outMesh2DSplit.size());
1982 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1983 for(std::size_t i=0;i<nbOfMeshes;i++)
1984 tmp[i]=outMesh2DSplit[i];
1986 ret1->getCoords()->setInfoOnComponents(compNames);
1987 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1988 // To finish - filter ret3 - std::numeric_limits<mcIdType>::max() -> -1 - negate values must be resolved.
1990 MCAuto<DataArrayIdType> edgesToDealWith(ret3->findIdsStrictlyNegative());
1991 for(const mcIdType *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1993 mcIdType old2DCellId(-ret3->getIJ(*it,0)-1);
1994 MCAuto<DataArrayIdType> candidates(ret2->findIdsEqual(old2DCellId));
1995 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
1997 ret3->changeValue(std::numeric_limits<mcIdType>::max(),-1);
2000 splitMesh1D=ret1.retn();
2001 splitMesh2D=ret2D.retn();
2002 cellIdInMesh2D=ret2.retn();
2003 cellIdInMesh1D=ret3.retn();
2007 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
2008 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2009 * 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
2010 * 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).
2011 * 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.
2013 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
2014 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
2016 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
2017 * This method expects that all nodes in \a this are not closer than \a eps.
2018 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2020 * \param [in] eps the relative error to detect merged edges.
2021 * \return DataArrayIdType * - 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
2022 * that the user is expected to deal with.
2024 * \throw If \a this is not coherent.
2025 * \throw If \a this has not spaceDim equal to 2.
2026 * \throw If \a this has not meshDim equal to 2.
2027 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
2029 DataArrayIdType *MEDCouplingUMesh::conformize2D(double eps)
2031 static const int SPACEDIM=2;
2032 checkConsistencyLight();
2033 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2034 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2035 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),descIndx1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescIndx1(DataArrayIdType::New());
2036 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
2037 const mcIdType *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
2038 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2039 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
2040 mcIdType nCell=getNumberOfCells(),nDescCell=mDesc->getNumberOfCells();
2041 std::vector< std::vector<mcIdType> > intersectEdge(nDescCell),overlapEdge(nDescCell);
2042 std::vector<double> addCoo;
2043 BBTree<SPACEDIM,mcIdType> myTree(bbox,0,0,nDescCell,-eps);
2044 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2045 for(mcIdType i=0;i<nDescCell;i++)
2047 std::vector<mcIdType> candidates;
2048 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
2049 for(std::vector<mcIdType>::const_iterator it=candidates.begin();it!=candidates.end();it++)
2050 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
2052 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
2053 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
2054 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
2055 INTERP_KERNEL::MergePoints merge;
2056 INTERP_KERNEL::QuadraticPolygon c1,c2;
2057 e1->intersectWith(e2,merge,c1,c2);
2058 e1->decrRef(); e2->decrRef();
2059 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
2060 overlapEdge[i].push_back(*it);
2061 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
2062 overlapEdge[*it].push_back(i);
2065 // splitting done. sort intersect point in intersectEdge.
2066 std::vector< std::vector<mcIdType> > middle(nDescCell);
2067 mcIdType nbOf2DCellsToBeSplit(0);
2068 bool middleNeedsToBeUsed(false);
2069 std::vector<bool> cells2DToTreat(nDescCell,false);
2070 for(mcIdType i=0;i<nDescCell;i++)
2072 std::vector<mcIdType>& isect(intersectEdge[i]);
2073 std::size_t sz(isect.size());
2076 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
2077 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
2078 e->sortSubNodesAbs(coords,isect);
2083 mcIdType idx0(rdi[i]),idx1(rdi[i+1]);
2085 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
2086 if(!cells2DToTreat[rd[idx0]])
2088 cells2DToTreat[rd[idx0]]=true;
2089 nbOf2DCellsToBeSplit++;
2091 // try to reuse at most eventual 'middle' of SEG3
2092 std::vector<mcIdType>& mid(middle[i]);
2093 mid.resize(sz+1,-1);
2094 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
2096 middleNeedsToBeUsed=true;
2097 const std::vector<mcIdType>& candidates(overlapEdge[i]);
2098 std::vector<mcIdType> trueCandidates;
2099 for(std::vector<mcIdType>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
2100 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
2101 trueCandidates.push_back(*itc);
2102 mcIdType stNode(c[ci[i]+1]),endNode(isect[0]);
2103 for(std::size_t j=0;j<sz+1;j++)
2105 for(std::vector<mcIdType>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
2107 mcIdType tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
2108 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
2109 { mid[j]=*itc; break; }
2112 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
2117 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()),notRet(DataArrayIdType::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
2118 if(nbOf2DCellsToBeSplit==0)
2121 mcIdType *retPtr(ret->getPointer());
2122 for(mcIdType i=0;i<nCell;i++)
2123 if(cells2DToTreat[i])
2126 MCAuto<DataArrayIdType> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
2127 DataArrayIdType *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
2128 DataArrayIdType::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
2129 DataArrayIdType::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
2130 if(middleNeedsToBeUsed)
2131 { DataArrayIdType::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
2132 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
2133 mcIdType nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
2134 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.
2135 setPartOfMySelf(ret->begin(),ret->end(),*modif);
2137 bool areNodesMerged; mcIdType newNbOfNodes;
2138 if(nbOfNodesCreated!=0)
2139 MCAuto<DataArrayIdType> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2145 * 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.
2146 * 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).
2147 * 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
2148 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2149 * 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
2150 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2152 * 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
2153 * using new instance, idem for coordinates.
2155 * 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.
2157 * \return DataArrayIdType * - The list of cellIds in \a this that have at least one edge colinearized.
2159 * \throw If \a this is not coherent.
2160 * \throw If \a this has not spaceDim equal to 2.
2161 * \throw If \a this has not meshDim equal to 2.
2163 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2165 DataArrayIdType *MEDCouplingUMesh::colinearize2D(double eps)
2167 return internalColinearize2D(eps, false);
2171 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2172 * 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
2173 * merged, contrary to colinearize2D().
2175 * \sa MEDCouplingUMesh::colinearize2D
2177 DataArrayIdType *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2179 return internalColinearize2D(eps, true);
2184 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2186 DataArrayIdType *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2188 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
2189 checkConsistencyLight();
2190 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2191 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2192 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2193 mcIdType nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2194 const mcIdType *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2195 MCAuto<DataArrayIdType> newc(DataArrayIdType::New()),newci(DataArrayIdType::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2196 std::map<mcIdType, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2199 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2200 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),descI1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescI1(DataArrayIdType::New());
2201 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2202 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
2203 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2204 MCAuto<DataArrayIdType> dsi(revDescI2->deltaShiftIndex());
2205 MCAuto<DataArrayIdType> ids(dsi->findIdsGreaterThan(2));
2206 const mcIdType * cPtr(mDesc0D->getNodalConnectivity()->begin());
2207 for(auto it = ids->begin(); it != ids->end(); it++)
2208 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2211 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2212 const double *coords(_coords->begin());
2213 mcIdType *newciptr(newci->getPointer());
2214 for(mcIdType i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2216 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2217 ret->pushBackSilent(i);
2218 newciptr[1]=newc->getNumberOfTuples();
2223 if(!appendedCoords->empty())
2225 appendedCoords->rearrange(2);
2226 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2228 setCoords(newCoords);
2231 setConnectivity(newc,newci,true);
2237 * c, cI describe a wire mesh in 3D space, in local numbering
2238 * startNode, endNode in global numbering
2239 *\return true if the segment is indeed split
2241 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, mcIdType startNode, mcIdType endNode,
2242 const mcIdType * c, const mcIdType * cI, const mcIdType *idsBg, const mcIdType *endBg,
2243 std::vector<mcIdType> & pointIds, std::vector<mcIdType> & hitSegs)
2245 using namespace std;
2247 const int SPACEDIM=3;
2248 typedef pair<double, mcIdType> PairDI;
2250 for (const mcIdType * it = idsBg; it != endBg; ++it)
2252 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2253 mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2];
2254 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2255 x.insert(make_pair(coo[end*SPACEDIM], end));
2258 vector<PairDI> xx(x.begin(), x.end());
2259 sort(xx.begin(),xx.end());
2260 pointIds.reserve(xx.size());
2262 // Keep what is inside [startNode, endNode]:
2264 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2266 const mcIdType idx = (*it).second;
2269 if (idx == startNode) go = 1;
2270 if (idx == endNode) go = 2;
2271 if (go) pointIds.push_back(idx);
2274 pointIds.push_back(idx);
2275 if (idx == endNode || idx == startNode)
2279 // vector<mcIdType> pointIds2(pointIds.size()+2);
2280 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2281 // pointIds2[0] = startNode;
2282 // pointIds2[pointIds2.size()-1] = endNode;
2285 reverse(pointIds.begin(), pointIds.end());
2287 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2288 for (const mcIdType * it = idsBg; it != endBg; ++it)
2290 mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2];
2291 vector<mcIdType>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2292 if (itStart == pointIds.end()) continue;
2293 vector<mcIdType>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2294 if (itEnd == pointIds.end()) continue;
2295 if (abs(distance(itEnd, itStart)) != 1) continue;
2296 hitSegs.push_back(*it); // segment is undivided.
2299 return (pointIds.size() > 2); // something else apart start and end node
2302 void MEDCouplingUMesh::ReplaceEdgeInFace(const mcIdType * sIdxConn, const mcIdType * sIdxConnE, mcIdType startNode, mcIdType endNode,
2303 const std::vector<mcIdType>& insidePoints, std::vector<mcIdType>& modifiedFace)
2305 using namespace std;
2306 size_t dst = distance(sIdxConn, sIdxConnE);
2307 modifiedFace.reserve(dst + insidePoints.size()-2);
2308 modifiedFace.resize(dst);
2309 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2311 vector<mcIdType>::iterator shortEnd = modifiedFace.begin()+dst;
2312 vector<mcIdType>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2313 if (startPos == shortEnd)
2314 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2315 vector<mcIdType>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2316 if (endPos == shortEnd)
2317 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2318 size_t d = distance(startPos, endPos);
2319 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2320 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2322 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2329 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2330 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2331 * This method performs a conformization of \b this.
2333 * Only polyhedron cells are supported. You can call convertAllToPoly()
2335 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2336 * This method expects that all nodes in \a this are not closer than \a eps.
2337 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2339 * \b WARNING this method only works for 'partition-like' non-conformities. When joining adjacent faces, the set of all small faces must fit exactly into the
2340 * neighboring bigger face. No real face intersection is actually computed.
2342 * \param [in] eps the relative error to detect merged edges.
2343 * \return DataArrayIdType * - 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
2344 * that the user is expected to deal with.
2346 * \throw If \a this is not coherent.
2347 * \throw If \a this has not spaceDim equal to 3.
2348 * \throw If \a this has not meshDim equal to 3.
2349 * \throw If the 'partition-like' condition described above is not respected.
2350 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2352 DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps)
2354 using namespace std;
2356 static const int SPACEDIM=3;
2357 checkConsistencyLight();
2358 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2359 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2360 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2361 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2363 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2364 const double * coo(_coords->begin());
2365 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
2368 /*************************
2370 *************************/
2371 MCAuto<DataArrayIdType> descDNU(DataArrayIdType::New()),descIDNU(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
2372 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2373 const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2374 const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2375 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2378 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2379 const double *bbox(bboxArr->begin()); getCoords()->begin();
2380 mcIdType nDescCell=mDesc->getNumberOfCells();
2381 BBTree<SPACEDIM,mcIdType> myTree(bbox,0,0,nDescCell,-eps);
2382 // Surfaces - handle biggest first
2383 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2384 DataArrayDouble * surfs = surfF->getArray();
2386 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2387 DataArrayDouble * normals = normalsF->getArray();
2388 const double * normalsP = normals->getConstPointer();
2390 // Sort faces by decreasing surface:
2391 vector< pair<double,mcIdType> > S;
2392 for(mcIdType i=0;i < surfs->getNumberOfTuples();i++)
2394 pair<double,mcIdType> p = make_pair(surfs->begin()[i], i);
2397 sort(S.rbegin(),S.rend()); // reverse sort
2398 vector<bool> hit(nDescCell);
2399 fill(hit.begin(), hit.end(), false);
2400 vector<mcIdType> hitPoly; // the final result: which 3D cells have been modified.
2402 for( vector<pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
2404 mcIdType faceIdx = (*it).second;
2405 if (hit[faceIdx]) continue;
2407 vector<mcIdType> candidates, cands2;
2408 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2409 // Keep only candidates whose normal matches the normal of current face
2410 for(vector<mcIdType>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2412 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2413 if (*it2 != faceIdx && col)
2414 cands2.push_back(*it2);
2419 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2420 INTERP_KERNEL::TranslationRotationMatrix rotation;
2421 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2422 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2423 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2425 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2426 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2427 double * cooPartRef(mPartRef->_coords->getPointer());
2428 double * cooPartCand(mPartCand->_coords->getPointer());
2429 for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2430 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2431 for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2432 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2434 // Localize faces in 2D thanks to barycenters
2435 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2436 vector<std::size_t> compo; compo.push_back(2);
2437 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2438 MCAuto<DataArrayIdType> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2439 if (!idsGoodPlane->getNumberOfTuples())
2442 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2444 compo[0] = 0; compo.push_back(1);
2445 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2446 mPartRef->changeSpaceDimension(2,0.0);
2447 MCAuto<DataArrayIdType> cc(DataArrayIdType::New()), ccI(DataArrayIdType::New());
2448 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2450 if (!cc->getNumberOfTuples())
2452 MCAuto<DataArrayIdType> dsi = ccI->deltaShiftIndex();
2455 MCAuto<DataArrayIdType> tmp = dsi->findIdsInRange(0, 2);
2456 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2459 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2460 throw INTERP_KERNEL::Exception(oss.str());
2464 MCAuto<DataArrayIdType> ids = dsi->findIdsEqual(1);
2466 if (!ids->getNumberOfTuples())
2469 double checkSurf=0.0;
2470 const mcIdType * idsGoodPlaneP(idsGoodPlane->begin());
2471 for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++)
2473 mcIdType faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2474 hit[faceIdx2] = true;
2475 checkSurf += surfs->begin()[faceIdx2];
2477 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2480 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2481 throw INTERP_KERNEL::Exception(oss.str());
2484 // For all polyhedrons using this face, replace connectivity:
2485 vector<mcIdType> polyIndices, packsIds, facePack;
2486 for (mcIdType ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2487 polyIndices.push_back(revDescP[ii]);
2488 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2490 // Current face connectivity
2491 const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2492 const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2493 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2494 // Deletion of old faces
2496 for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2498 if (packsIds[jj] == -1)
2499 // The below should never happen - if a face is used several times, with a different layout of the nodes
2500 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2501 // faces which are actually used only once, by a single cell. This is different for edges below.
2502 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2504 connSla->deletePack(*it2, packsIds[jj]);
2506 // Insertion of new faces:
2507 for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++)
2509 // Build pack from the face to insert:
2510 mcIdType faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2511 mcIdType facePack2Sz;
2512 const mcIdType * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2513 // Insert it in all hit polyhedrons:
2514 for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2515 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2520 // Set back modified connectivity
2521 MCAuto<DataArrayIdType> cAuto; cAuto.takeRef(_nodal_connec);
2522 MCAuto<DataArrayIdType> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2523 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2526 /************************
2528 ************************/
2529 // Now we have a face-conform mesh.
2531 // Recompute descending
2532 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
2533 // Rebuild desc connectivity with orientation this time!!
2534 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2535 const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2536 const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2537 const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2538 MCAuto<DataArrayIdType> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2539 MCAuto<DataArrayIdType> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2540 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2541 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
2542 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2543 // std::cout << "writing!\n";
2544 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2545 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2546 const mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2547 const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2548 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2549 const double *bbox2(bboxArr->begin());
2550 mcIdType nDesc2Cell=mDesc2->getNumberOfCells();
2551 BBTree<SPACEDIM,mcIdType> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2553 // Edges - handle longest first
2554 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2555 DataArrayDouble * lens = lenF->getArray();
2557 // Sort edges by decreasing length:
2558 vector<pair<double,mcIdType> > S;
2559 for(mcIdType i=0;i < lens->getNumberOfTuples();i++)
2561 pair<double,mcIdType> p = make_pair(lens->getIJ(i, 0), i);
2564 sort(S.rbegin(),S.rend()); // reverse sort
2566 vector<bool> hit(nDesc2Cell);
2567 fill(hit.begin(), hit.end(), false);
2569 for( vector<pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
2571 mcIdType eIdx = (*it).second;
2575 vector<mcIdType> candidates, cands2;
2576 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2577 // Keep only candidates colinear with current edge
2579 mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2580 for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2581 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2582 for(vector<mcIdType>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2585 mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2586 for (mcIdType i3=0; i3 < 3; i3++)
2587 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2588 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2589 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2590 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2592 cands2.push_back(*it2);
2594 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2597 // Now rotate edges to bring them on Ox
2598 mcIdType startNode = cDesc2[cIDesc2[eIdx]+1];
2599 mcIdType endNode = cDesc2[cIDesc2[eIdx]+2];
2600 INTERP_KERNEL::TranslationRotationMatrix rotation;
2601 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2602 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2603 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2604 MCAuto<DataArrayIdType> nodeMap(mPartCand->zipCoordsTraducer());
2605 mcIdType nbElemsNotM1;
2607 MCAuto<DataArrayIdType> tmp(nodeMap->findIdsNotEqual(-1));
2608 nbElemsNotM1 = tmp->getNbOfElems();
2610 MCAuto<DataArrayIdType> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2611 double * cooPartRef(mPartRef->_coords->getPointer());
2612 double * cooPartCand(mPartCand->_coords->getPointer());
2613 for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2614 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2615 for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2616 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2619 // Eliminate all edges for which y or z is not null
2620 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2621 vector<std::size_t> compo; compo.push_back(1);
2622 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2624 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2625 MCAuto<DataArrayIdType> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2626 MCAuto<DataArrayIdType> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2627 MCAuto<DataArrayIdType> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2628 if (!idsGoodLine->getNumberOfTuples())
2631 // Now the ordering along the Ox axis:
2632 std::vector<mcIdType> insidePoints, hitSegs;
2633 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2634 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2635 idsGoodLine->begin(), idsGoodLine->end(),
2636 /*out*/insidePoints, hitSegs);
2637 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2638 for (vector<mcIdType>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2639 hit[cands2[*its]] = true;
2641 if (!isSplit) // current segment remains in one piece
2644 // Get original node IDs in global coords array
2645 for (std::vector<mcIdType>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2646 *iit = nodeMapInv->begin()[*iit];
2648 vector<mcIdType> polyIndices, packsIds, facePack;
2649 // For each face implying this edge
2650 for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2652 mcIdType faceIdx = revDescP2[ii];
2653 // each cell where this face is involved connectivity will be modified:
2654 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2656 // Current face connectivity
2657 const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2658 const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2660 vector<mcIdType> modifiedFace;
2661 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2662 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2663 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2667 // Rebuild 3D connectivity from descending:
2668 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2669 MCAuto<DataArrayIdType> superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2670 MCAuto<DataArrayIdType> idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0);
2671 MCAuto<DataArrayIdType> vals(DataArrayIdType::New()); vals->alloc(0);
2672 newConn->set3(superIdx, idx, vals);
2673 mcIdType nbCells=getNumberOfCells();
2674 for(mcIdType ii = 0; ii < nbCells; ii++)
2675 for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2677 mcIdType sz, faceIdx = abs(descP[jj])-1;
2678 bool orient = descP[jj]>0;
2679 const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2681 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2684 vector<mcIdType> rev(sz-1);
2685 for (mcIdType kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2686 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2690 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2693 ret = ret->buildUniqueNotSorted();