1 // Copyright (C) 2007-2020 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++)
559 MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<mcIdType> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<mcIdType,mcIdType>& mergedNodes, const std::vector< std::vector<mcIdType> >& colinear2, const std::vector< std::vector<mcIdType> >& intersectEdge1,
560 MCAuto<DataArrayIdType>& idsInRetColinear, MCAuto<DataArrayIdType>& idsInMesh1DForIdsInRetColinear)
562 idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1);
563 idsInMesh1DForIdsInRetColinear=DataArrayIdType::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
564 mcIdType nCells=mesh1D->getNumberOfCells();
565 if(nCells!=ToIdType(intersectEdge2.size()))
566 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
567 const DataArrayDouble *coo2(mesh1D->getCoords());
568 const mcIdType *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
569 const double *coo2Ptr(coo2->begin());
570 mcIdType offset1(coords1->getNumberOfTuples());
571 mcIdType offset2(offset1+coo2->getNumberOfTuples());
572 mcIdType offset3(offset2+ToIdType(addCoo.size())/2);
573 std::vector<double> addCooQuad;
574 MCAuto<DataArrayIdType> cOut(DataArrayIdType::New()),ciOut(DataArrayIdType::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
575 mcIdType tmp[4],cicnt(0),kk(0);
576 for(mcIdType i=0;i<nCells;i++)
578 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
579 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
580 const std::vector<mcIdType>& subEdges(intersectEdge2[i]);
581 mcIdType nbSubEdge=ToIdType(subEdges.size()/2);
582 for(mcIdType j=0;j<nbSubEdge;j++,kk++)
584 MCAuto<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
585 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
586 INTERP_KERNEL::Edge *e2Ptr(e2);
587 std::map<mcIdType,mcIdType>::const_iterator itm;
588 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
590 tmp[0]=INTERP_KERNEL::NORM_SEG3;
591 itm=mergedNodes.find(subEdges[2*j]);
592 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
593 itm=mergedNodes.find(subEdges[2*j+1]);
594 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
595 tmp[3]=offset3+ToIdType(addCooQuad.size()/2);
597 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
599 cOut->insertAtTheEnd(tmp,tmp+4);
600 ciOut->pushBackSilent(cicnt);
604 tmp[0]=INTERP_KERNEL::NORM_SEG2;
605 itm=mergedNodes.find(subEdges[2*j]);
606 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
607 itm=mergedNodes.find(subEdges[2*j+1]);
608 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
610 cOut->insertAtTheEnd(tmp,tmp+3);
611 ciOut->pushBackSilent(cicnt);
614 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
616 idsInRetColinear->pushBackSilent(kk);
617 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
622 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
623 ret->setConnectivity(cOut,ciOut,true);
624 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
625 arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2);
626 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,addCooQuad.size()/2,2);
627 std::vector<const DataArrayDouble *> coordss(4);
628 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
629 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
634 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
636 std::vector<mcIdType> allEdges;
637 for(const mcIdType *it2(descBg);it2!=descEnd;it2++)
639 const std::vector<mcIdType>& edge1(intersectEdge1[std::abs(*it2)-1]);
641 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
643 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
645 std::size_t nb(allEdges.size());
647 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
648 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
649 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
650 ret->setCoords(coords);
651 ret->allocateCells(1);
652 std::vector<mcIdType> connOut(nbOfEdgesOf2DCellSplit);
653 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
654 connOut[kk]=allEdges[2*kk];
655 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(connOut.size()),&connOut[0]);
659 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
661 const mcIdType *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
662 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
664 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
665 if(sz!=std::distance(descBg,descEnd))
666 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
667 INTERP_KERNEL::AutoPtr<mcIdType> tmpPtr(new mcIdType[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
668 std::vector<mcIdType> allEdges,centers;
669 const double *coordsPtr(coords->begin());
670 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
671 mcIdType offset(coords->getNumberOfTuples());
672 for(const mcIdType *it2(descBg);it2!=descEnd;it2++,ii++)
674 INTERP_KERNEL::NormalizedCellType typeOfSon;
675 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
676 const std::vector<mcIdType>& edge1(intersectEdge1[std::abs(*it2)-1]);
678 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
680 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
682 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
684 {//the current edge has been subsplit -> create corresponding centers.
685 mcIdType nbOfCentersToAppend=ToIdType(edge1.size()/2);
686 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
687 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
688 std::vector<mcIdType>::const_iterator it3(allEdges.end()-edge1.size());
689 for(mcIdType k=0;k<nbOfCentersToAppend;k++)
692 const double *aa(coordsPtr+2*(*it3++));
693 const double *bb(coordsPtr+2*(*it3++));
694 ee->getMiddleOfPoints(aa,bb,tmpp);
695 addCoo->insertAtTheEnd(tmpp,tmpp+2);
696 centers.push_back(offset+k);
700 std::size_t nb(allEdges.size());
702 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
703 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
704 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
706 ret->setCoords(coords);
709 addCoo->rearrange(2);
710 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
711 ret->setCoords(addCoo);
713 ret->allocateCells(1);
714 std::vector<mcIdType> connOut(nbOfEdgesOf2DCellSplit);
715 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
716 connOut[kk]=allEdges[2*kk];
717 connOut.insert(connOut.end(),centers.begin(),centers.end());
718 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(connOut.size()),&connOut[0]);
723 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
726 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
728 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1)
730 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
731 if(!cm.isQuadratic())
732 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
734 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
737 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<mcIdType>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
740 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
742 const INTERP_KERNEL::Edge *ee(*it);
743 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
747 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(conn.size()),&conn[0]);
750 const double *coo(mesh2D->getCoords()->begin());
751 std::size_t sz(conn.size());
752 std::vector<double> addCoo;
753 std::vector<mcIdType> conn2(conn);
754 mcIdType offset(mesh2D->getNumberOfNodes());
755 for(std::size_t i=0;i<sz;i++)
758 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
759 addCoo.insert(addCoo.end(),tmp,tmp+2);
760 conn2.push_back(offset+ToIdType(i));
762 mesh2D->getCoords()->rearrange(1);
763 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
764 mesh2D->getCoords()->rearrange(2);
765 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(conn2.size()),&conn2[0]);
770 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
772 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
773 * a set of edges defined in \a splitMesh1D.
775 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<mcIdType>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
776 std::vector< std::vector<mcIdType> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
778 std::size_t nb(edge1Bis.size()/2);
779 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
780 mcIdType iEnd=splitMesh1D->getNumberOfCells();
782 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
784 const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
785 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
786 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
789 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
790 out0.resize(1); out1.resize(1);
791 std::vector<mcIdType>& connOut(out0[0]);
792 connOut.resize(nbOfEdgesOf2DCellSplit);
793 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
794 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
795 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
797 connOut[kk]=edge1Bis[2*kk];
798 edgesPtr[kk]=edge1BisPtr[2*kk];
803 // [i,iEnd[ contains the
804 out0.resize(2); out1.resize(2);
805 std::vector<mcIdType>& connOutLeft(out0[0]);
806 std::vector<mcIdType>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
807 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
808 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
809 for(std::size_t k=ii;k<jj+1;k++)
810 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
811 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
812 for(mcIdType ik=0;ik<iEnd;ik++)
814 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
815 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
818 for(mcIdType ik=iEnd-1;ik>=0;ik--)
819 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
820 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
821 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
822 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
823 for(mcIdType ik=0;ik<iEnd;ik++)
824 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
825 eright.insert(eright.end(),ees.begin(),ees.end());
833 CellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
835 std::vector<mcIdType> _edges;
836 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
839 CellInfo::CellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
841 std::size_t nbe(edges.size());
842 std::vector<mcIdType> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
843 for(std::size_t i=0;i<nbe;i++)
845 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
846 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
848 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
849 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
850 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
856 EdgeInfo(mcIdType istart, mcIdType iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
857 EdgeInfo(mcIdType istart, mcIdType iend, mcIdType pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
858 bool isInMyRange(mcIdType pos) const { return pos>=_istart && pos<_iend; }
859 void somethingHappendAt(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
860 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, mcIdType offset, mcIdType neighbors[2]) const;
864 MCAuto<MEDCouplingUMesh> _mesh;
865 MCAuto<INTERP_KERNEL::Edge> _edge;
866 mcIdType _left; // index (local numbering) of the left 2D cell bordering the edge '_edge'
867 mcIdType _right; // same as above, right side.
871 * Update indices of left and right 2D cell bordering the current edge.
873 void EdgeInfo::somethingHappendAt(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
875 const MEDCouplingUMesh *mesh(_mesh);
881 { _left++; _right++; return ; }
882 if (_right > pos && _left != pos)
883 { _right++; return ; }
886 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
887 if((isLeft && isRight) || (!isLeft && !isRight))
888 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
899 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
900 if((isLeft && isRight) || (!isLeft && !isRight))
901 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
916 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, mcIdType offset, mcIdType neighbors[2]) const
918 const MEDCouplingUMesh *mesh(_mesh);
921 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
924 {// not fully splitting cell case
925 if(mesh2D->getNumberOfCells()==1)
926 {//little optimization. 1 cell no need to find in which cell mesh is !
927 neighbors[0]=offset; neighbors[1]=offset;
932 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
933 mcIdType cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
935 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
936 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
941 class VectorOfCellInfo
944 VectorOfCellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
945 std::size_t size() const { return _pool.size(); }
946 mcIdType getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
947 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);
948 const std::vector<mcIdType>& getConnOf(mcIdType pos) const { return get(pos)._edges; }
949 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(mcIdType pos) const { return get(pos)._edges_ptr; }
950 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
951 void feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const;
953 mcIdType getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const;
954 void updateEdgeInfo(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
955 const CellInfo& get(mcIdType pos) const;
956 CellInfo& get(mcIdType pos);
958 std::vector<CellInfo> _pool; // for a newly created 2D cell, the list of edges ToIdType( and edges ptr constiuing it
959 MCAuto<MEDCouplingUMesh> _ze_mesh; // the aggregated mesh
960 std::vector<EdgeInfo> _edge_info; // for each new edge added when cuting the 2D cell, the information on left and right bordering 2D cell
963 VectorOfCellInfo::VectorOfCellInfo(const std::vector<mcIdType>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
965 _pool[0]._edges=edges;
966 _pool[0]._edges_ptr=edgesPtr;
969 mcIdType VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
972 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
975 const MEDCouplingUMesh *zeMesh(_ze_mesh);
977 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
978 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
979 return zeMesh->getCellContainingPoint(barys->begin(),eps);
982 void VectorOfCellInfo::setMeshAt(mcIdType pos, const MCAuto<MEDCouplingUMesh>& mesh, mcIdType istart, mcIdType iend,
983 const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<mcIdType> >& edges,
984 const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
986 get(pos);//to check pos
987 bool isFast(pos==0 && _pool.size()==1);
988 std::size_t sz(edges.size());
989 // dealing with edges
991 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
993 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
995 std::vector<CellInfo> pool(_pool.size()-1+sz);
996 for(mcIdType i=0;i<pos;i++)
998 for(std::size_t j=0;j<sz;j++)
999 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
1000 for(std::size_t i=pos+1;i<_pool.size();i++)
1001 pool[i+sz-1]=_pool[i];
1005 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
1013 std::vector< MCAuto<MEDCouplingUMesh> > ms;
1016 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
1020 if(pos<_ze_mesh->getNumberOfCells()-1)
1022 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
1025 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
1026 for(std::size_t j=0;j<ms2.size();j++)
1028 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
1031 void VectorOfCellInfo::feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const
1033 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
1036 mcIdType VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const
1039 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
1041 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
1043 if((*it).isInMyRange(pos))
1046 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
1049 void VectorOfCellInfo::updateEdgeInfo(mcIdType pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
1051 get(pos);//to perform the sanity check;
1052 if(_edge_info.empty())
1054 std::size_t sz(_edge_info.size()-1);
1055 for(std::size_t i=0;i<sz;i++)
1056 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
1059 const CellInfo& VectorOfCellInfo::get(mcIdType pos) const
1061 if(pos<0 || pos>=ToIdType(_pool.size()))
1062 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1066 CellInfo& VectorOfCellInfo::get(mcIdType pos)
1068 if(pos<0 || pos>=ToIdType(_pool.size()))
1069 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1075 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1076 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1078 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1080 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1082 * \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.
1084 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<mcIdType>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, mcIdType offset,
1085 MCAuto<DataArrayIdType>& idsLeftRight)
1087 mcIdType nbCellsInSplitMesh1D=splitMesh1D->getNumberOfCells();
1088 if(nbCellsInSplitMesh1D==0)
1089 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1090 const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1091 std::size_t nb(allEdges.size()),jj;
1093 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1094 std::vector<mcIdType> edge1Bis(nb*2);
1095 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1096 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1097 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1098 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1099 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1101 idsLeftRight=DataArrayIdType::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1102 mcIdType *idsLeftRightPtr(idsLeftRight->getPointer());
1103 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1105 // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1106 // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1107 // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1108 // of the connectivity.
1109 MCAuto <DataArrayIdType> renumb(DataArrayIdType::New());
1110 renumb->alloc(nbCellsInSplitMesh1D,1);
1111 const mcIdType * renumbP(renumb->begin());
1113 mcIdType i, first=cSplitPtr[1];
1114 // Follow 1D line backward as long as it is connected:
1115 for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1116 first=cSplitPtr[ciSplitPtr[i]+1];
1117 if (i < nbCellsInSplitMesh1D-1)
1119 // Build circular permutation to shift consecutive edges together
1121 renumb->applyModulus(nbCellsInSplitMesh1D);
1122 splitMesh1D->renumberCells(renumbP, false);
1123 cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1124 ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1129 // The 1D first piece is used to intersect the 2D cell resulting in max two 2D cells.
1130 // The next 1D piece is localised (getPositionOf()) into this previous cut.
1131 // The result of the next intersection replaces the former single 2D cell that has been cut in the
1132 // pool. The neighbourhood information detained by pool._edge_info is also updated so that left and right
1133 // adjacent 2D cell of a 1D piece is kept up to date.
1134 // And so on and so forth.
1135 for(mcIdType iStart=0;iStart<nbCellsInSplitMesh1D;)
1136 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1137 mcIdType iEnd(iStart);
1138 for(;iEnd<nbCellsInSplitMesh1D;)
1140 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1146 if(iEnd<nbCellsInSplitMesh1D)
1149 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1150 mcIdType pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1152 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1153 retTmp->setCoords(splitMesh1D->getCoords());
1154 retTmp->allocateCells();
1156 std::vector< std::vector<mcIdType> > out0;
1157 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1159 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1160 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1161 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1162 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1166 for(mcIdType mm=0;mm<nbCellsInSplitMesh1D;mm++)
1167 pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1169 return pool.getZeMesh().retn();
1173 * splitMesh1D is an input parameter but might have its cells renumbered.
1175 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, mcIdType cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1176 const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector<mcIdType> >& intersectEdge1, mcIdType offset,
1177 MCAuto<DataArrayIdType>& idsLeftRight)
1179 const mcIdType *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1181 std::vector<mcIdType> allEdges;
1182 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1183 for(const mcIdType *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1185 mcIdType edgeId(std::abs(*it)-1);
1186 std::map< MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
1187 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1188 const std::vector<mcIdType>& edge1(intersectEdge1[edgeId]);
1190 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1192 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1193 std::size_t sz(edge1.size());
1194 for(std::size_t cnt=0;cnt<sz;cnt++)
1195 allEdgesPtr.push_back(ee);
1198 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1201 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const mcIdType *conn1, const INTERP_KERNEL::CellModel& typ2, const mcIdType *conn2, double eps)
1203 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1204 {//easy case comparison not
1205 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1207 else if(typ1.isQuadratic() && typ2.isQuadratic())
1209 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1212 if(conn1[2]==conn2[2])
1214 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1215 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1219 {//only one is quadratic
1220 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1223 const double *a(0),*bb(0),*be(0);
1224 if(typ1.isQuadratic())
1226 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1230 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1232 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1233 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1239 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1240 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1242 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1244 mcIdType FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const mcIdType *candidatesIn2DBg, const mcIdType *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, mcIdType cellIdInMesh1DSplitRelative, double eps)
1246 if(candidatesIn2DEnd==candidatesIn2DBg)
1247 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1248 const double *coo(mesh2DSplit->getCoords()->begin());
1249 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1250 return *candidatesIn2DBg;
1251 mcIdType edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1252 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1253 if(cellIdInMesh1DSplitRelative<0)
1254 cur1D->changeOrientationOfCells();
1255 const mcIdType *c1D(cur1D->getNodalConnectivity()->begin());
1256 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1257 for(const mcIdType *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1259 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1260 const mcIdType *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1261 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1262 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1263 INTERP_KERNEL::AutoPtr<mcIdType> tmpPtr(new mcIdType[ci[1]-ci[0]]);
1264 for(unsigned it2=0;it2<sz;it2++)
1266 INTERP_KERNEL::NormalizedCellType typeOfSon;
1267 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1268 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1269 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1273 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1277 * \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.
1278 * 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.
1279 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1280 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1281 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1282 * \param [out] addCoo - nodes to be append at the end
1283 * \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.
1285 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1286 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)
1288 static const int SPACEDIM=2;
1289 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1290 const mcIdType *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1291 // Build BB tree of all edges in the tool mesh (second mesh)
1292 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1293 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1294 mcIdType nDescCell1=m1Desc->getNumberOfCells(),nDescCell2=m2Desc->getNumberOfCells();
1295 intersectEdge1.resize(nDescCell1);
1296 colinear2.resize(nDescCell2);
1297 subDiv2.resize(nDescCell2);
1298 BBTree<SPACEDIM,mcIdType> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1299 BBTreePts<SPACEDIM,mcIdType> treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps);
1301 std::vector<mcIdType> candidates1(1);
1302 mcIdType offset1(m1Desc->getNumberOfNodes());
1303 mcIdType offset2(offset1+m2Desc->getNumberOfNodes());
1304 for(mcIdType i=0;i<nDescCell1;i++) // for all edges in the first mesh
1306 std::vector<mcIdType> candidates2; // edges of mesh2 candidate for intersection
1307 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1308 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1310 std::map<INTERP_KERNEL::Node *,mcIdType> map1,map2;
1311 std::map<mcIdType, INTERP_KERNEL::Node *> revMap2;
1312 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1313 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1315 for (auto& kv : map2)
1316 revMap2[kv.second] = kv.first;
1318 // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged.
1319 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2);
1320 // 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
1321 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1322 std::set<INTERP_KERNEL::Node *> nodes;
1323 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1324 std::size_t szz(nodes.size());
1325 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1326 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1327 for(std::size_t iii=0;iii<szz;iii++,itt++)
1328 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1329 // end of protection
1330 // Performs edge cutting:
1331 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1336 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1337 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1343 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1344 * It builds the descending connectivity of the two meshes, and then using a binary tree
1345 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1346 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1348 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1349 std::vector< std::vector<mcIdType> >& intersectEdge1, std::vector< std::vector<mcIdType> >& colinear2, std::vector< std::vector<mcIdType> >& subDiv2,
1350 MEDCouplingUMesh *& m1Desc, DataArrayIdType *&desc1, DataArrayIdType *&descIndx1, DataArrayIdType *&revDesc1, DataArrayIdType *&revDescIndx1,
1351 std::vector<double>& addCoo,
1352 MEDCouplingUMesh *& m2Desc, DataArrayIdType *&desc2, DataArrayIdType *&descIndx2, DataArrayIdType *&revDesc2, DataArrayIdType *&revDescIndx2)
1354 // Build desc connectivity
1355 desc1=DataArrayIdType::New(); descIndx1=DataArrayIdType::New(); revDesc1=DataArrayIdType::New(); revDescIndx1=DataArrayIdType::New();
1356 desc2=DataArrayIdType::New();
1357 descIndx2=DataArrayIdType::New();
1358 revDesc2=DataArrayIdType::New();
1359 revDescIndx2=DataArrayIdType::New();
1360 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1361 MCAuto<DataArrayIdType> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1362 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1363 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1364 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1365 std::map<mcIdType,mcIdType> notUsedMap;
1366 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1367 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1368 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1372 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1373 * (newly created) nodes corresponding to the edge intersections.
1375 * @param[out] cr connectivity of the resulting mesh
1376 * @param[out] crI connectivity of the resulting mesh
1377 * @param[out] cNb1 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1378 * @param[out] cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1379 * TODO: describe input parameters
1381 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const mcIdType *desc1, const mcIdType *descIndx1,
1382 const std::vector<std::vector<mcIdType> >& intesctEdges1, const std::vector< std::vector<mcIdType> >& colinear2,
1383 const MEDCouplingUMesh *m2, const mcIdType *desc2, const mcIdType *descIndx2, const std::vector<std::vector<mcIdType> >& intesctEdges2,
1384 const std::vector<double>& addCoords,
1385 std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& cr, std::vector<mcIdType>& crI, std::vector<mcIdType>& cNb1, std::vector<mcIdType>& cNb2)
1387 static const int SPACEDIM=2;
1388 const double *coo1(m1->getCoords()->begin());
1389 const mcIdType *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1390 mcIdType offset1(m1->getNumberOfNodes());
1391 const double *coo2(m2->getCoords()->begin());
1392 const mcIdType *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1393 mcIdType offset2(offset1+m2->getNumberOfNodes());
1394 mcIdType offset3(offset2+ToIdType(addCoords.size())/2);
1395 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1396 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1397 // Here a BBTree on 2D-cells, not on segments:
1398 BBTree<SPACEDIM,mcIdType> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1399 mcIdType ncell1=m1->getNumberOfCells();
1401 for(mcIdType i=0;i<ncell1;i++)
1403 std::vector<mcIdType> candidates2;
1404 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1405 std::map<INTERP_KERNEL::Node *,mcIdType> mapp;
1406 std::map<mcIdType,INTERP_KERNEL::Node *> mappRev;
1407 INTERP_KERNEL::QuadraticPolygon pol1;
1408 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1409 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1410 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1411 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1412 // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1413 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1414 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1416 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
1417 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1418 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1419 for(it1.first();!it1.finished();it1.next())
1420 edges1.insert(it1.current()->getPtr());
1422 std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1423 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1425 // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1426 // Again all the additional intersecting nodes are there.
1427 for(std::vector<mcIdType>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1429 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1430 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1431 // Complete mapping with elements coming from the current cell it2 in mesh2:
1432 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1433 // pol2 is the new QP in the final merged result.
1434 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1435 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1437 // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1438 for (auto &p: pol2s)
1439 p.cleanDegeneratedConsecutiveEdges();
1440 edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges.
1443 // Now rebuild intersected cells from all this:
1444 for(std::vector<mcIdType>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1446 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1447 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1448 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1449 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1451 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1452 // by m2 but that we still want to keep in the final result.
1457 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1459 catch(INTERP_KERNEL::Exception& e)
1461 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();
1462 throw INTERP_KERNEL::Exception(oss.str());
1465 for(std::map<mcIdType,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1466 (*it).second->decrRef();
1470 void InsertNodeInConnIfNecessary(mcIdType nodeIdToInsert, std::vector<mcIdType>& conn, const double *coords, double eps)
1472 std::vector<mcIdType>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1475 std::size_t sz(conn.size());
1476 std::size_t found(std::numeric_limits<std::size_t>::max());
1477 for(std::size_t i=0;i<sz;i++)
1479 mcIdType pt0(conn[i]),pt1(conn[(i+1)%sz]);
1480 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]};
1481 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1482 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1483 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1485 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];
1486 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]);
1488 if(dotTest>eps && dotTest<1.-eps)
1494 if(found==std::numeric_limits<std::size_t>::max())
1495 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1496 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1499 void SplitIntoToPart(const std::vector<mcIdType>& conn, mcIdType pt0, mcIdType pt1, std::vector<mcIdType>& part0, std::vector<mcIdType>& part1)
1501 std::size_t sz(conn.size());
1502 std::vector<mcIdType> *curPart(&part0);
1503 for(std::size_t i=0;i<sz;i++)
1505 mcIdType nextt(conn[(i+1)%sz]);
1506 (*curPart).push_back(nextt);
1507 if(nextt==pt0 || nextt==pt1)
1513 (*curPart).push_back(nextt);
1519 * 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.
1521 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1522 const mcIdType *desc, const mcIdType *descIndx, const double *coords, double eps,
1523 std::vector<std::vector<mcIdType> >& res) const
1525 checkFullyDefined();
1526 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1527 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1528 const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1529 mcIdType nbOfCells=getNumberOfCells();
1531 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1532 for(mcIdType i=0;i<nbOfCells;i++)
1534 mcIdType offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1535 for(int j=0;j<nbOfFaces;j++)
1537 const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1538 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1539 mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1540 INTERP_KERNEL::AutoPtr<mcIdType> tmp(new mcIdType[sz]);
1541 INTERP_KERNEL::NormalizedCellType cmsId;
1542 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1543 std::vector<mcIdType> elt((mcIdType *)tmp,(mcIdType *)tmp+nbOfNodesSon);
1544 if(p.first!=-1 && p.second!=-1)
1548 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1549 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1550 std::vector<mcIdType> elt1,elt2;
1551 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1552 res.push_back(elt1);
1553 res.push_back(elt2);
1565 * 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).
1567 * \sa MEDCouplingUMesh::split2DCells
1569 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI)
1571 checkConnectivityFullyDefined();
1572 mcIdType ncells=getNumberOfCells();
1573 mcIdType lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1574 MCAuto<DataArrayIdType> c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach);
1575 const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1576 mcIdType *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1577 mcIdType prevPosOfCi(ciPtr[0]);
1578 for(mcIdType i=0;i<ncells;i++,ciPtr++,descIPtr++)
1580 mcIdType offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1581 *cPtr++=ToIdType(INTERP_KERNEL::NORM_POLYGON); *cPtr++=oldConn[prevPosOfCi+1];
1582 for(mcIdType j=0;j<sz;j++)
1584 mcIdType offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1585 for(mcIdType k=0;k<sz2;k++)
1586 *cPtr++=subPtr[offset2+k];
1588 *cPtr++=oldConn[prevPosOfCi+j+2];
1591 prevPosOfCi=ciPtr[1];
1592 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1595 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1596 _nodal_connec->decrRef();
1597 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1602 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1604 * \return mcIdType - the number of new nodes created.
1605 * \sa MEDCouplingUMesh::split2DCells
1607 mcIdType MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayIdType *desc, const DataArrayIdType *descI, const DataArrayIdType *subNodesInSeg, const DataArrayIdType *subNodesInSegI, const DataArrayIdType *mid, const DataArrayIdType *midI)
1609 checkConsistencyLight();
1610 mcIdType ncells=getNumberOfCells();
1611 mcIdType lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples());
1612 mcIdType nodesCnt(getNumberOfNodes());
1613 MCAuto<DataArrayIdType> c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach);
1614 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1615 const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1616 const mcIdType *midPtr(mid->begin()),*midIPtr(midI->begin());
1617 const double *oldCoordsPtr(getCoords()->begin());
1618 mcIdType *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1619 mcIdType prevPosOfCi(ciPtr[0]);
1620 for(mcIdType i=0;i<ncells;i++,ciPtr++,descIPtr++)
1622 mcIdType offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1623 for(mcIdType j=0;j<sz;j++)
1624 { mcIdType sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1625 *cPtr++=ToIdType(INTERP_KERNEL::NORM_QPOLYG); cPtr[0]=oldConn[prevPosOfCi+1];
1626 for(mcIdType j=0;j<sz;j++)//loop over subedges of oldConn
1628 mcIdType offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1632 cPtr[1]=oldConn[prevPosOfCi+2+j];
1633 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1636 std::vector<INTERP_KERNEL::Node *> ns(3);
1637 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1638 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1639 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1640 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1641 for(mcIdType k=0;k<sz2;k++)//loop over subsplit of current subedge
1643 cPtr[1]=subPtr[offset2+k];
1644 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1646 mcIdType tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1649 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1651 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1652 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1655 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1656 _nodal_connec->decrRef();
1657 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1658 addCoo->rearrange(2);
1659 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1661 return addCoo->getNumberOfTuples();
1668 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1669 * returns a result mesh constituted by polygons.
1670 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1671 * all nodes from m2.
1672 * The meshes should be in 2D space. In
1673 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1675 * \b WARNING: the two meshes should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and
1676 * orientCorrectly2DCells() can be used for this.
1677 * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells)
1678 * \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
1679 * 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)
1680 * \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
1681 * 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)
1682 * \param [in] eps - precision used to detect coincident mesh entities.
1683 * \param [out] cellNb1 - a new instance of DataArrayIdType holding for each result
1684 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1685 * this array using decrRef() as it is no more needed.
1686 * \param [out] cellNb2 - a new instance of DataArrayIdType holding for each result
1687 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1688 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1689 * any cell of \a m2. The caller is to delete this array using decrRef() as
1690 * it is no more needed.
1691 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1692 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1693 * is no more needed.
1694 * \throw If the coordinates array is not set in any of the meshes.
1695 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1696 * \throw If any of the meshes is not a 2D mesh in 2D space.
1698 * \sa conformize2D, mergeNodes
1700 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1701 double eps, DataArrayIdType *&cellNb1, DataArrayIdType *&cellNb2)
1704 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1705 m1->checkFullyDefined();
1706 m2->checkFullyDefined();
1707 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1708 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1709 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1711 // Step 1: compute all edge intersections (new nodes)
1712 std::vector< std::vector<mcIdType> > intersectEdge1, colinear2, subDiv2;
1713 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1714 DataArrayIdType *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1715 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1716 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1717 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1718 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1719 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1720 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1721 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1723 // Step 2: re-order newly created nodes according to the ordering found in m2
1724 std::vector< std::vector<mcIdType> > intersectEdge2;
1725 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1726 subDiv2.clear(); dd5=0; dd6=0;
1729 std::vector<mcIdType> cr,crI; //no DataArrayIdType because interface with Geometric2D
1730 std::vector<mcIdType> cNb1,cNb2; //no DataArrayIdType because interface with Geometric2D
1731 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1732 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1734 // Step 4: Prepare final result:
1735 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1736 addCooDa->alloc(addCoo.size()/2,2);
1737 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1738 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1739 addCoordsQuadraticDa->alloc(addCoordsQuadratic.size()/2,2);
1740 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1741 std::vector<const DataArrayDouble *> coordss(4);
1742 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1743 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1744 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1745 MCAuto<DataArrayIdType> conn(DataArrayIdType::New()); conn->alloc(cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1746 MCAuto<DataArrayIdType> connI(DataArrayIdType::New()); connI->alloc(crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1747 MCAuto<DataArrayIdType> c1(DataArrayIdType::New()); c1->alloc(cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1748 MCAuto<DataArrayIdType> c2(DataArrayIdType::New()); c2->alloc(cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1749 ret->setConnectivity(conn,connI,true);
1750 ret->setCoords(coo);
1751 cellNb1=c1.retn(); cellNb2=c2.retn();
1756 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1757 * 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
1758 * and finally, in case of quadratic polygon the centers of edges new nodes.
1759 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1761 * \b WARNING: the 2D mesh should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and
1762 * orientCorrectly2DCells() can be used for this.
1763 * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells)
1765 * \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
1766 * 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)
1767 * \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
1768 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1769 * \param [in] eps - precision used to perform intersections and localization operations.
1770 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1771 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1772 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1773 * 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.
1774 * \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
1775 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1776 * 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.
1778 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1780 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayIdType *&cellIdInMesh2D, DataArrayIdType *&cellIdInMesh1D)
1782 if(!mesh2D || !mesh1D)
1783 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1784 mesh2D->checkFullyDefined();
1785 mesh1D->checkFullyDefined();
1786 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1787 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1788 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1789 // Step 1: compute all edge intersections (new nodes)
1790 std::vector< std::vector<mcIdType> > intersectEdge1, colinear2, subDiv2;
1791 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1792 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1794 // Build desc connectivity
1795 DataArrayIdType *desc1(DataArrayIdType::New()),*descIndx1(DataArrayIdType::New()),*revDesc1(DataArrayIdType::New()),*revDescIndx1(DataArrayIdType::New());
1796 MCAuto<DataArrayIdType> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1797 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1798 std::map<mcIdType,mcIdType> mergedNodes;
1799 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1800 // use mergeNodes to fix intersectEdge1
1801 for(std::vector< std::vector<mcIdType> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1803 std::size_t n((*it0).size()/2);
1804 mcIdType eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1805 std::map<mcIdType,mcIdType>::const_iterator it1;
1806 it1=mergedNodes.find(eltStart);
1807 if(it1!=mergedNodes.end())
1808 (*it0)[0]=(*it1).second;
1809 it1=mergedNodes.find(eltEnd);
1810 if(it1!=mergedNodes.end())
1811 (*it0)[2*n-1]=(*it1).second;
1814 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1815 addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2);
1816 // Step 2: re-order newly created nodes according to the ordering found in m2
1817 std::vector< std::vector<mcIdType> > intersectEdge2;
1818 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1820 // Step 3: compute splitMesh1D
1821 MCAuto<DataArrayIdType> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1822 MCAuto<DataArrayIdType> ret2(DataArrayIdType::New()); ret2->alloc(0,1);
1823 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1824 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1825 MCAuto<DataArrayIdType> ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<mcIdType>::max()); ret3->rearrange(2);
1826 MCAuto<DataArrayIdType> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1827 // deal with cells in mesh2D that are not cut but only some of their edges are
1828 MCAuto<DataArrayIdType> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1829 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1830 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1831 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
1832 if(!idsInDesc2DToBeRefined->empty())
1834 DataArrayIdType *out0(0),*outi0(0);
1835 DataArrayIdType::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1836 MCAuto<DataArrayIdType> outi0s(outi0);
1838 out0s=out0s->buildUnique();
1842 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1843 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1844 MCAuto<DataArrayIdType> elts,eltsIndex;
1845 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1846 MCAuto<DataArrayIdType> eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1);
1847 if (eltsIndex->getNumberOfTuples() > 1)
1848 eltsIndex2 = eltsIndex->deltaShiftIndex();
1849 MCAuto<DataArrayIdType> eltsIndex3(eltsIndex2->findIdsEqual(1));
1850 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1851 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1852 MCAuto<DataArrayIdType> cellsToBeModified(elts->buildUnique());
1853 MCAuto<DataArrayIdType> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1854 if((DataArrayIdType *)out0s)
1855 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1856 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1857 // OK all is ready to insert in ret2 mesh
1858 if(!untouchedCells->empty())
1859 {// the most easy part, cells in mesh2D not impacted at all
1860 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1861 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1862 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1864 if((DataArrayIdType *)out0s)
1865 {// here dealing with cells in out0s but not in cellsToBeModified
1866 MCAuto<DataArrayIdType> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1867 const mcIdType *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1868 for(const mcIdType *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1870 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1871 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1873 mcIdType offset(ret2->getNumberOfTuples());
1874 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1875 MCAuto<DataArrayIdType> partOfRet3(DataArrayIdType::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1876 partOfRet3->fillWithValue(std::numeric_limits<mcIdType>::max()); partOfRet3->rearrange(2);
1877 mcIdType kk(0),*ret3ptr(partOfRet3->getPointer());
1878 for(const mcIdType *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1880 mcIdType faceId(std::abs(*it)-1);
1881 for(const mcIdType *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1883 mcIdType tmp(fewModifiedCells->findIdFirstEqual(*it2));
1886 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1887 ret3ptr[2*kk]=tmp+offset;
1888 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1889 ret3ptr[2*kk+1]=tmp+offset;
1892 {//the current edge is shared by a 2D cell that will be split just after
1893 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1894 ret3ptr[2*kk]=-(*it2+1);
1895 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1896 ret3ptr[2*kk+1]=-(*it2+1);
1900 m1Desc->setCoords(ret1->getCoords());
1901 ret1NonCol->setCoords(ret1->getCoords());
1902 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1903 if(!outMesh2DSplit.empty())
1905 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1906 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1907 (*itt)->setCoords(da);
1910 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1911 for(const mcIdType *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1913 MCAuto<DataArrayIdType> idsNonColPerCell(elts->findIdsEqual(*it));
1914 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1915 MCAuto<DataArrayIdType> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1916 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1917 MCAuto<DataArrayIdType> partOfRet3;
1918 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));
1919 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1920 outMesh2DSplit.push_back(splitOfOneCell);
1921 for(mcIdType i=0;i<splitOfOneCell->getNumberOfCells();i++)
1922 ret2->pushBackSilent(*it);
1925 std::size_t nbOfMeshes(outMesh2DSplit.size());
1926 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1927 for(std::size_t i=0;i<nbOfMeshes;i++)
1928 tmp[i]=outMesh2DSplit[i];
1930 ret1->getCoords()->setInfoOnComponents(compNames);
1931 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1932 // To finish - filter ret3 - std::numeric_limits<mcIdType>::max() -> -1 - negate values must be resolved.
1934 MCAuto<DataArrayIdType> edgesToDealWith(ret3->findIdsStrictlyNegative());
1935 for(const mcIdType *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1937 mcIdType old2DCellId(-ret3->getIJ(*it,0)-1);
1938 MCAuto<DataArrayIdType> candidates(ret2->findIdsEqual(old2DCellId));
1939 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
1941 ret3->changeValue(std::numeric_limits<mcIdType>::max(),-1);
1944 splitMesh1D=ret1.retn();
1945 splitMesh2D=ret2D.retn();
1946 cellIdInMesh2D=ret2.retn();
1947 cellIdInMesh1D=ret3.retn();
1951 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1952 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1953 * 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
1954 * 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).
1955 * 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.
1957 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1958 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1960 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1961 * This method expects that all nodes in \a this are not closer than \a eps.
1962 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1964 * \param [in] eps the relative error to detect merged edges.
1965 * \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
1966 * that the user is expected to deal with.
1968 * \throw If \a this is not coherent.
1969 * \throw If \a this has not spaceDim equal to 2.
1970 * \throw If \a this has not meshDim equal to 2.
1971 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1973 DataArrayIdType *MEDCouplingUMesh::conformize2D(double eps)
1975 static const int SPACEDIM=2;
1976 checkConsistencyLight();
1977 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1978 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1979 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),descIndx1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescIndx1(DataArrayIdType::New());
1980 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1981 const mcIdType *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1982 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1983 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1984 mcIdType nCell=getNumberOfCells(),nDescCell=mDesc->getNumberOfCells();
1985 std::vector< std::vector<mcIdType> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1986 std::vector<double> addCoo;
1987 BBTree<SPACEDIM,mcIdType> myTree(bbox,0,0,nDescCell,-eps);
1988 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1989 for(mcIdType i=0;i<nDescCell;i++)
1991 std::vector<mcIdType> candidates;
1992 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1993 for(std::vector<mcIdType>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1994 if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice
1996 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
1997 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1998 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1999 INTERP_KERNEL::MergePoints merge;
2000 INTERP_KERNEL::QuadraticPolygon c1,c2;
2001 e1->intersectWith(e2,merge,c1,c2);
2002 e1->decrRef(); e2->decrRef();
2003 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
2004 overlapEdge[i].push_back(*it);
2005 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
2006 overlapEdge[*it].push_back(i);
2009 // splitting done. sort intersect point in intersectEdge.
2010 std::vector< std::vector<mcIdType> > middle(nDescCell);
2011 mcIdType nbOf2DCellsToBeSplit(0);
2012 bool middleNeedsToBeUsed(false);
2013 std::vector<bool> cells2DToTreat(nDescCell,false);
2014 for(mcIdType i=0;i<nDescCell;i++)
2016 std::vector<mcIdType>& isect(intersectEdge[i]);
2017 std::size_t sz(isect.size());
2020 std::map<MCAuto<INTERP_KERNEL::Node>,mcIdType> m;
2021 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
2022 e->sortSubNodesAbs(coords,isect);
2027 mcIdType idx0(rdi[i]),idx1(rdi[i+1]);
2029 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
2030 if(!cells2DToTreat[rd[idx0]])
2032 cells2DToTreat[rd[idx0]]=true;
2033 nbOf2DCellsToBeSplit++;
2035 // try to reuse at most eventual 'middle' of SEG3
2036 std::vector<mcIdType>& mid(middle[i]);
2037 mid.resize(sz+1,-1);
2038 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
2040 middleNeedsToBeUsed=true;
2041 const std::vector<mcIdType>& candidates(overlapEdge[i]);
2042 std::vector<mcIdType> trueCandidates;
2043 for(std::vector<mcIdType>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
2044 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
2045 trueCandidates.push_back(*itc);
2046 mcIdType stNode(c[ci[i]+1]),endNode(isect[0]);
2047 for(std::size_t j=0;j<sz+1;j++)
2049 for(std::vector<mcIdType>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
2051 mcIdType tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
2052 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
2053 { mid[j]=*itc; break; }
2056 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
2061 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()),notRet(DataArrayIdType::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
2062 if(nbOf2DCellsToBeSplit==0)
2065 mcIdType *retPtr(ret->getPointer());
2066 for(mcIdType i=0;i<nCell;i++)
2067 if(cells2DToTreat[i])
2070 MCAuto<DataArrayIdType> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
2071 DataArrayIdType *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
2072 DataArrayIdType::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
2073 DataArrayIdType::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
2074 if(middleNeedsToBeUsed)
2075 { DataArrayIdType::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
2076 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
2077 mcIdType nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
2078 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.
2079 setPartOfMySelf(ret->begin(),ret->end(),*modif);
2081 bool areNodesMerged; mcIdType newNbOfNodes;
2082 if(nbOfNodesCreated!=0)
2083 MCAuto<DataArrayIdType> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2089 * 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.
2090 * 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).
2091 * 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
2092 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2093 * 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
2094 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2096 * 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
2097 * using new instance, idem for coordinates.
2099 * 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.
2101 * \return DataArrayIdType * - The list of cellIds in \a this that have at least one edge colinearized.
2103 * \throw If \a this is not coherent.
2104 * \throw If \a this has not spaceDim equal to 2.
2105 * \throw If \a this has not meshDim equal to 2.
2107 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2109 DataArrayIdType *MEDCouplingUMesh::colinearize2D(double eps)
2111 return internalColinearize2D(eps, false);
2115 * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2116 * 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
2117 * merged, contrary to colinearize2D().
2119 * \sa MEDCouplingUMesh::colinearize2D
2121 DataArrayIdType *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2123 return internalColinearize2D(eps, true);
2128 * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2130 DataArrayIdType *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2132 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
2133 checkConsistencyLight();
2134 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2135 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2136 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2137 mcIdType nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2138 const mcIdType *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2139 MCAuto<DataArrayIdType> newc(DataArrayIdType::New()),newci(DataArrayIdType::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2140 std::map<mcIdType, bool> forbiddenPoints; // list of points that can not be removed (or it will break conformity)
2143 // A point that is used by more than 2 edges can not be removed without breaking conformity:
2144 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),descI1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescI1(DataArrayIdType::New());
2145 MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2146 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
2147 MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2148 MCAuto<DataArrayIdType> dsi(revDescI2->deltaShiftIndex());
2149 MCAuto<DataArrayIdType> ids(dsi->findIdsGreaterThan(2));
2150 const mcIdType * cPtr(mDesc0D->getNodalConnectivity()->begin());
2151 for(auto it = ids->begin(); it != ids->end(); it++)
2152 forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2155 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2156 const double *coords(_coords->begin());
2157 mcIdType *newciptr(newci->getPointer());
2158 for(mcIdType i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2160 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2161 ret->pushBackSilent(i);
2162 newciptr[1]=newc->getNumberOfTuples();
2167 if(!appendedCoords->empty())
2169 appendedCoords->rearrange(2);
2170 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2172 setCoords(newCoords);
2175 setConnectivity(newc,newci,true);
2181 * c, cI describe a wire mesh in 3D space, in local numbering
2182 * startNode, endNode in global numbering
2183 *\return true if the segment is indeed split
2185 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, mcIdType startNode, mcIdType endNode,
2186 const mcIdType * c, const mcIdType * cI, const mcIdType *idsBg, const mcIdType *endBg,
2187 std::vector<mcIdType> & pointIds, std::vector<mcIdType> & hitSegs)
2189 using namespace std;
2191 const int SPACEDIM=3;
2192 typedef pair<double, mcIdType> PairDI;
2194 for (const mcIdType * it = idsBg; it != endBg; ++it)
2196 assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2197 mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2];
2198 x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
2199 x.insert(make_pair(coo[end*SPACEDIM], end));
2202 vector<PairDI> xx(x.begin(), x.end());
2203 sort(xx.begin(),xx.end());
2204 pointIds.reserve(xx.size());
2206 // Keep what is inside [startNode, endNode]:
2208 for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2210 const mcIdType idx = (*it).second;
2213 if (idx == startNode) go = 1;
2214 if (idx == endNode) go = 2;
2215 if (go) pointIds.push_back(idx);
2218 pointIds.push_back(idx);
2219 if (idx == endNode || idx == startNode)
2223 // vector<mcIdType> pointIds2(pointIds.size()+2);
2224 // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2225 // pointIds2[0] = startNode;
2226 // pointIds2[pointIds2.size()-1] = endNode;
2229 reverse(pointIds.begin(), pointIds.end());
2231 // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2232 for (const mcIdType * it = idsBg; it != endBg; ++it)
2234 mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2];
2235 vector<mcIdType>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2236 if (itStart == pointIds.end()) continue;
2237 vector<mcIdType>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2238 if (itEnd == pointIds.end()) continue;
2239 if (abs(distance(itEnd, itStart)) != 1) continue;
2240 hitSegs.push_back(*it); // segment is undivided.
2243 return (pointIds.size() > 2); // something else apart start and end node
2246 void MEDCouplingUMesh::ReplaceEdgeInFace(const mcIdType * sIdxConn, const mcIdType * sIdxConnE, mcIdType startNode, mcIdType endNode,
2247 const std::vector<mcIdType>& insidePoints, std::vector<mcIdType>& modifiedFace)
2249 using namespace std;
2250 size_t dst = distance(sIdxConn, sIdxConnE);
2251 modifiedFace.reserve(dst + insidePoints.size()-2);
2252 modifiedFace.resize(dst);
2253 copy(sIdxConn, sIdxConnE, modifiedFace.data());
2255 vector<mcIdType>::iterator shortEnd = modifiedFace.begin()+dst;
2256 vector<mcIdType>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2257 if (startPos == shortEnd)
2258 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2259 vector<mcIdType>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2260 if (endPos == shortEnd)
2261 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2262 size_t d = distance(startPos, endPos);
2263 if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2264 modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted.
2266 modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2273 * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2274 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2275 * This method performs a conformization of \b this.
2277 * Only polyhedron cells are supported. You can call convertAllToPoly()
2279 * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2280 * This method expects that all nodes in \a this are not closer than \a eps.
2281 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2283 * \param [in] eps the relative error to detect merged edges.
2284 * \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
2285 * that the user is expected to deal with.
2287 * \throw If \a this is not coherent.
2288 * \throw If \a this has not spaceDim equal to 3.
2289 * \throw If \a this has not meshDim equal to 3.
2290 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2292 DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps)
2294 using namespace std;
2296 static const int SPACEDIM=3;
2297 checkConsistencyLight();
2298 if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2299 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2300 if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2301 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2303 MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2304 const double * coo(_coords->begin());
2305 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
2308 /*************************
2310 *************************/
2311 MCAuto<DataArrayIdType> descDNU(DataArrayIdType::New()),descIDNU(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
2312 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2313 const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2314 const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2315 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2318 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2319 const double *bbox(bboxArr->begin()); getCoords()->begin();
2320 mcIdType nDescCell=mDesc->getNumberOfCells();
2321 BBTree<SPACEDIM,mcIdType> myTree(bbox,0,0,nDescCell,-eps);
2322 // Surfaces - handle biggest first
2323 MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2324 DataArrayDouble * surfs = surfF->getArray();
2326 MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2327 DataArrayDouble * normals = normalsF->getArray();
2328 const double * normalsP = normals->getConstPointer();
2330 // Sort faces by decreasing surface:
2331 vector< pair<double,mcIdType> > S;
2332 for(mcIdType i=0;i < surfs->getNumberOfTuples();i++)
2334 pair<double,mcIdType> p = make_pair(surfs->begin()[i], i);
2337 sort(S.rbegin(),S.rend()); // reverse sort
2338 vector<bool> hit(nDescCell);
2339 fill(hit.begin(), hit.end(), false);
2340 vector<mcIdType> hitPoly; // the final result: which 3D cells have been modified.
2342 for( vector<pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
2344 mcIdType faceIdx = (*it).second;
2345 if (hit[faceIdx]) continue;
2347 vector<mcIdType> candidates, cands2;
2348 myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2349 // Keep only candidates whose normal matches the normal of current face
2350 for(vector<mcIdType>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2352 bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2353 if (*it2 != faceIdx && col)
2354 cands2.push_back(*it2);
2359 // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2360 INTERP_KERNEL::TranslationRotationMatrix rotation;
2361 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2362 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2363 coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2365 MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false)); // false=zipCoords is called
2366 MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2367 double * cooPartRef(mPartRef->_coords->getPointer());
2368 double * cooPartCand(mPartCand->_coords->getPointer());
2369 for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2370 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2371 for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2372 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2374 // Localize faces in 2D thanks to barycenters
2375 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2376 vector<std::size_t> compo; compo.push_back(2);
2377 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2378 MCAuto<DataArrayIdType> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2379 if (!idsGoodPlane->getNumberOfTuples())
2382 baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2384 compo[0] = 0; compo.push_back(1);
2385 MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2386 mPartRef->changeSpaceDimension(2,0.0);
2387 MCAuto<DataArrayIdType> cc(DataArrayIdType::New()), ccI(DataArrayIdType::New());
2388 mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2390 if (!cc->getNumberOfTuples())
2392 MCAuto<DataArrayIdType> dsi = ccI->deltaShiftIndex();
2395 MCAuto<DataArrayIdType> tmp = dsi->findIdsInRange(0, 2);
2396 if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2399 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2400 throw INTERP_KERNEL::Exception(oss.str());
2404 MCAuto<DataArrayIdType> ids = dsi->findIdsEqual(1);
2406 if (!ids->getNumberOfTuples())
2409 double checkSurf=0.0;
2410 const mcIdType * idsGoodPlaneP(idsGoodPlane->begin());
2411 for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++)
2413 mcIdType faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2414 hit[faceIdx2] = true;
2415 checkSurf += surfs->begin()[faceIdx2];
2417 if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2420 oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2421 throw INTERP_KERNEL::Exception(oss.str());
2424 // For all polyhedrons using this face, replace connectivity:
2425 vector<mcIdType> polyIndices, packsIds, facePack;
2426 for (mcIdType ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2427 polyIndices.push_back(revDescP[ii]);
2428 ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2430 // Current face connectivity
2431 const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2432 const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2433 connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2434 // Deletion of old faces
2436 for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2438 if (packsIds[jj] == -1)
2439 // The below should never happen - if a face is used several times, with a different layout of the nodes
2440 // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2441 // faces which are actually used only once, by a single cell. This is different for edges below.
2442 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2444 connSla->deletePack(*it2, packsIds[jj]);
2446 // Insertion of new faces:
2447 for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++)
2449 // Build pack from the face to insert:
2450 mcIdType faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2451 mcIdType facePack2Sz;
2452 const mcIdType * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2453 // Insert it in all hit polyhedrons:
2454 for (vector<mcIdType>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2455 connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type
2460 // Set back modified connectivity
2461 MCAuto<DataArrayIdType> cAuto; cAuto.takeRef(_nodal_connec);
2462 MCAuto<DataArrayIdType> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2463 connSla->convertToPolyhedronConn(cAuto, cIAuto);
2466 /************************
2468 ************************/
2469 // Now we have a face-conform mesh.
2471 // Recompute descending
2472 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New());
2473 // Rebuild desc connectivity with orientation this time!!
2474 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2475 const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2476 const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2477 const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2478 MCAuto<DataArrayIdType> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2479 MCAuto<DataArrayIdType> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2480 MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2481 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New());
2482 MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2483 // std::cout << "writing!\n";
2484 // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2485 // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2486 const mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2487 const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2488 MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2489 const double *bbox2(bboxArr->begin());
2490 mcIdType nDesc2Cell=mDesc2->getNumberOfCells();
2491 BBTree<SPACEDIM,mcIdType> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2493 // Edges - handle longest first
2494 MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2495 DataArrayDouble * lens = lenF->getArray();
2497 // Sort edges by decreasing length:
2498 vector<pair<double,mcIdType> > S;
2499 for(mcIdType i=0;i < lens->getNumberOfTuples();i++)
2501 pair<double,mcIdType> p = make_pair(lens->getIJ(i, 0), i);
2504 sort(S.rbegin(),S.rend()); // reverse sort
2506 vector<bool> hit(nDesc2Cell);
2507 fill(hit.begin(), hit.end(), false);
2509 for( vector<pair<double,mcIdType> >::const_iterator it = S.begin(); it != S.end(); it++)
2511 mcIdType eIdx = (*it).second;
2515 vector<mcIdType> candidates, cands2;
2516 myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2517 // Keep only candidates colinear with current edge
2519 mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2520 for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
2521 vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2522 for(vector<mcIdType>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2525 mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2526 for (mcIdType i3=0; i3 < 3; i3++)
2527 vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2528 bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2529 // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2530 // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2532 cands2.push_back(*it2);
2534 if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
2537 // Now rotate edges to bring them on Ox
2538 mcIdType startNode = cDesc2[cIDesc2[eIdx]+1];
2539 mcIdType endNode = cDesc2[cIDesc2[eIdx]+2];
2540 INTERP_KERNEL::TranslationRotationMatrix rotation;
2541 INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2542 MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
2543 MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2544 MCAuto<DataArrayIdType> nodeMap(mPartCand->zipCoordsTraducer());
2545 mcIdType nbElemsNotM1;
2547 MCAuto<DataArrayIdType> tmp(nodeMap->findIdsNotEqual(-1));
2548 nbElemsNotM1 = tmp->getNbOfElems();
2550 MCAuto<DataArrayIdType> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2551 double * cooPartRef(mPartRef->_coords->getPointer());
2552 double * cooPartCand(mPartCand->_coords->getPointer());
2553 for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2554 rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2555 for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2556 rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2559 // Eliminate all edges for which y or z is not null
2560 MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2561 vector<std::size_t> compo; compo.push_back(1);
2562 MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2564 MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2565 MCAuto<DataArrayIdType> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2566 MCAuto<DataArrayIdType> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2567 MCAuto<DataArrayIdType> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2568 if (!idsGoodLine->getNumberOfTuples())
2571 // Now the ordering along the Ox axis:
2572 std::vector<mcIdType> insidePoints, hitSegs;
2573 bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2574 mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2575 idsGoodLine->begin(), idsGoodLine->end(),
2576 /*out*/insidePoints, hitSegs);
2577 // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2578 for (vector<mcIdType>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2579 hit[cands2[*its]] = true;
2581 if (!isSplit) // current segment remains in one piece
2584 // Get original node IDs in global coords array
2585 for (std::vector<mcIdType>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2586 *iit = nodeMapInv->begin()[*iit];
2588 vector<mcIdType> polyIndices, packsIds, facePack;
2589 // For each face implying this edge
2590 for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2592 mcIdType faceIdx = revDescP2[ii];
2593 // each cell where this face is involved connectivity will be modified:
2594 ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2596 // Current face connectivity
2597 const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2598 const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2600 vector<mcIdType> modifiedFace;
2601 ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2602 modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2603 connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2607 // Rebuild 3D connectivity from descending:
2608 MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2609 MCAuto<DataArrayIdType> superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2610 MCAuto<DataArrayIdType> idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0);
2611 MCAuto<DataArrayIdType> vals(DataArrayIdType::New()); vals->alloc(0);
2612 newConn->set3(superIdx, idx, vals);
2613 mcIdType nbCells=getNumberOfCells();
2614 for(mcIdType ii = 0; ii < nbCells; ii++)
2615 for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2617 mcIdType sz, faceIdx = abs(descP[jj])-1;
2618 bool orient = descP[jj]>0;
2619 const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2621 newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type
2624 vector<mcIdType> rev(sz-1);
2625 for (mcIdType kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2626 newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2630 newConn->convertToPolyhedronConn(cAuto, cIAuto);
2633 ret = ret->buildUniqueNotSorted();