1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingUMesh.hxx"
22 #include "MEDCoupling1GTUMesh.hxx"
23 #include "CellModel.hxx"
24 #include "VolSurfUser.txx"
25 #include "InterpolationUtils.hxx"
26 #include "PointLocatorAlgos.txx"
28 #include "BBTreeDst.txx"
29 #include "DirectedBoundingBox.hxx"
30 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
31 #include "InterpKernelAutoPtr.hxx"
32 #include "InterpKernelGeo2DNode.hxx"
33 #include "InterpKernelGeo2DEdgeLin.hxx"
34 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
35 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
44 using namespace MEDCoupling;
48 int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
54 int ret(nodesCnter++);
56 e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
57 addCoo.insertAtTheEnd(newPt,newPt+2);
62 int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
68 int ret(nodesCnter++);
70 e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt);
71 addCoo.insertAtTheEnd(newPt,newPt+2);
77 void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
80 int trueStart(start>=0?start:nbOfEdges+start);
81 tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
82 newConnOfCell->insertAtTheEnd(tmp,tmp+3);
87 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
88 InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2);
89 middles.push_back(tmp3+offset);
92 middles.push_back(connBg[trueStart+nbOfEdges]);
96 void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
98 int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
99 newConnOfCell->pushBackSilent(tmpEnd);
104 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
105 InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
106 middles.push_back(tmp3+offset);
109 middles.push_back(connBg[start+nbOfEdges]);
113 void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
115 // only the quadratic point to deal with:
120 int tmpSrt(connBg[start]),tmpEnd(connBg[stp]);
121 int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
122 InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
123 middles.push_back(tmp3+offset);
126 middles.push_back(connBg[start+nbOfEdges]);
130 void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
132 MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
133 std::map<MCAuto<INTERP_KERNEL::Node>,int>::const_iterator it(m.find(nTmp));
135 throw INTERP_KERNEL::Exception("Internal error in remapping !");
137 if(v==forbVal0 || v==forbVal1)
139 if(std::find(isect.begin(),isect.end(),v)==isect.end())
143 bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
148 bool presenceOfOn(false);
149 for(int i=0;i<sz;i++)
151 INTERP_KERNEL::ElementaryEdge *e(c[i]);
152 if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
154 IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
155 IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
161 namespace MEDCoupling
164 INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map< MCAuto<INTERP_KERNEL::Node>,int>& m)
166 INTERP_KERNEL::Edge *ret(0);
167 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]));
168 m[n0]=bg[0]; m[n1]=bg[1];
171 case INTERP_KERNEL::NORM_SEG2:
173 ret=new INTERP_KERNEL::EdgeLin(n0,n1);
176 case INTERP_KERNEL::NORM_SEG3:
178 INTERP_KERNEL::Node *n2(new INTERP_KERNEL::Node(coords2D[2*bg[2]],coords2D[2*bg[2]+1])); m[n2]=bg[2];
179 INTERP_KERNEL::EdgeLin *e1(new INTERP_KERNEL::EdgeLin(n0,n2)),*e2(new INTERP_KERNEL::EdgeLin(n2,n1));
180 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
181 // is the SEG3 degenerated, and thus can be reduced to a SEG2?
182 bool colinearity(inters.areColinears());
183 delete e1; delete e2;
185 { ret=new INTERP_KERNEL::EdgeLin(n0,n1); }
187 { ret=new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1); }
191 throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge2 : Expecting a mesh with spaceDim==2 and meshDim==1 !");
196 INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
198 INTERP_KERNEL::Edge *ret=0;
201 case INTERP_KERNEL::NORM_SEG2:
203 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
206 case INTERP_KERNEL::NORM_SEG3:
208 INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
209 INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
210 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
211 // is the SEG3 degenerated, and thus can be reduced to a SEG2?
212 bool colinearity=inters.areColinears();
213 delete e1; delete e2;
215 ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
217 ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
218 mapp2[bg[2]].second=false;
222 throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
228 * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from
229 * the global mesh 'mDesc'.
230 * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2.
231 * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
233 INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
234 std::map<INTERP_KERNEL::Node *,int>& mapp)
237 std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3.
238 const double *coo=mDesc->getCoords()->getConstPointer();
239 const int *c=mDesc->getNodalConnectivity()->getConstPointer();
240 const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
242 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
243 s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
244 for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
246 INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
247 mapp2[*it2]=std::pair<INTERP_KERNEL::Node *,bool>(n,true);
249 INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
250 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
252 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
253 ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
255 for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
257 if((*it2).second.second)
258 mapp[(*it2).second.first]=(*it2).first;
259 ((*it2).second.first)->decrRef();
264 INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
268 int locId=nodeId-offset2;
269 return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
273 int locId=nodeId-offset1;
274 return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
276 return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
280 * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI).
282 void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
283 const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
284 /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
286 for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
288 int eltId1=abs(*desc1)-1;
289 for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
291 std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
292 if(it==mappRev.end())
294 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
306 * 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 ) .
307 * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
309 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
311 std::size_t sz(std::distance(connBg,connEnd));
312 if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
313 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
315 INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
316 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
317 unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
318 unsigned nbOfHit(0); // number of fusions operated
319 int posBaseElt(0),posEndElt(0),nbOfTurn(0);
320 const unsigned int maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3; // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least
321 INTERP_KERNEL::NormalizedCellType typeOfSon;
322 std::vector<int> middles;
324 for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
326 cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
327 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
328 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
329 posEndElt = posBaseElt+1;
331 // Look backward first: are the final edges of the cells colinear with the first ones?
332 // This initializes posBaseElt.
335 for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell wih one single edge
337 cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
338 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
339 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
340 bool isColinear=eint->areColinears();
354 const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward
355 for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++) // 2nd condition is to avoid ending with a cell wih one single edge
357 cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #j's connectivity
358 INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
359 INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
360 bool isColinear(eint->areColinears());
372 //push [posBaseElt,posEndElt) in newConnOfCell using e
373 // The if clauses below are (volontary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its begining!
375 // at the begining of the connectivity (insert type)
376 EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
377 else if((nbOfHit+nbOfTurn) != (nbs-1))
379 EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
380 if ((nbOfHit+nbOfTurn) == (nbs-1))
381 // at the end (only quad points to deal with)
382 EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
383 posBaseElt=posEndElt;
387 newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
393 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
395 if(candidates.empty())
397 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
399 const std::vector<int>& pool(intersectEdge1[*it]);
400 int tmp[2]; tmp[0]=start; tmp[1]=stop;
401 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
406 tmp[0]=stop; tmp[1]=start;
407 if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
417 * This method performs the 2nd step of Partition of 2D mesh.
418 * This method has 4 inputs :
419 * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
420 * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
421 * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted.
422 * 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'
423 * Nodes end up lying consecutively on a cutted edge.
424 * \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.
425 * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
426 * \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.
427 * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
428 * \param[out] intersectEdge the same content as subDiv, but correclty oriented.
430 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
431 const std::vector<double>& addCoo,
432 const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
434 int offset1=m1->getNumberOfNodes();
435 int ncell=m2->getNumberOfCells();
436 const int *c=m2->getNodalConnectivity()->begin();
437 const int *cI=m2->getNodalConnectivityIndex()->begin();
438 const double *coo=m2->getCoords()->begin();
439 const double *cooBis=m1->getCoords()->begin();
440 int offset2=offset1+m2->getNumberOfNodes();
441 intersectEdge.resize(ncell);
442 for(int i=0;i<ncell;i++,cI++)
444 const std::vector<int>& divs=subDiv[i];
445 int nnode=cI[1]-cI[0]-1;
446 std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
447 std::map<INTERP_KERNEL::Node *, int> mapp22;
448 for(int j=0;j<nnode;j++)
450 INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
451 int nnid=c[(*cI)+j+1];
452 mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
453 mapp22[nn]=nnid+offset1;
455 INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
456 for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
457 ((*it).second.first)->decrRef();
458 std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
459 std::map<INTERP_KERNEL::Node *,int> mapp3;
460 for(std::size_t j=0;j<divs.size();j++)
463 INTERP_KERNEL::Node *tmp=0;
465 tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
467 tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
469 tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
473 e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
474 for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
480 MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<int> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<int,int>& mergedNodes, const std::vector< std::vector<int> >& colinear2, const std::vector< std::vector<int> >& intersectEdge1,
481 MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
483 idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
484 idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
485 int nCells(mesh1D->getNumberOfCells());
486 if(nCells!=(int)intersectEdge2.size())
487 throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
488 const DataArrayDouble *coo2(mesh1D->getCoords());
489 const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
490 const double *coo2Ptr(coo2->begin());
491 int offset1(coords1->getNumberOfTuples());
492 int offset2(offset1+coo2->getNumberOfTuples());
493 int offset3(offset2+addCoo.size()/2);
494 std::vector<double> addCooQuad;
495 MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
496 int tmp[4],cicnt(0),kk(0);
497 for(int i=0;i<nCells;i++)
499 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
500 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
501 const std::vector<int>& subEdges(intersectEdge2[i]);
502 int nbSubEdge(subEdges.size()/2);
503 for(int j=0;j<nbSubEdge;j++,kk++)
505 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));
506 MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
507 INTERP_KERNEL::Edge *e2Ptr(e2);
508 std::map<int,int>::const_iterator itm;
509 if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
511 tmp[0]=INTERP_KERNEL::NORM_SEG3;
512 itm=mergedNodes.find(subEdges[2*j]);
513 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
514 itm=mergedNodes.find(subEdges[2*j+1]);
515 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
516 tmp[3]=offset3+(int)addCooQuad.size()/2;
518 e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
520 cOut->insertAtTheEnd(tmp,tmp+4);
521 ciOut->pushBackSilent(cicnt);
525 tmp[0]=INTERP_KERNEL::NORM_SEG2;
526 itm=mergedNodes.find(subEdges[2*j]);
527 tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
528 itm=mergedNodes.find(subEdges[2*j+1]);
529 tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
531 cOut->insertAtTheEnd(tmp,tmp+3);
532 ciOut->pushBackSilent(cicnt);
535 if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
537 idsInRetColinear->pushBackSilent(kk);
538 idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
543 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
544 ret->setConnectivity(cOut,ciOut,true);
545 MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
546 arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
547 MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
548 std::vector<const DataArrayDouble *> coordss(4);
549 coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
550 MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
555 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
557 std::vector<int> allEdges;
558 for(const int *it2(descBg);it2!=descEnd;it2++)
560 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
562 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
564 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
566 std::size_t nb(allEdges.size());
568 throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
569 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
570 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
571 ret->setCoords(coords);
572 ret->allocateCells(1);
573 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
574 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
575 connOut[kk]=allEdges[2*kk];
576 ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
580 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
582 const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
583 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
585 unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
586 if(sz!=std::distance(descBg,descEnd))
587 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
588 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
589 std::vector<int> allEdges,centers;
590 const double *coordsPtr(coords->begin());
591 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
592 int offset(coords->getNumberOfTuples());
593 for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
595 INTERP_KERNEL::NormalizedCellType typeOfSon;
596 cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
597 const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
599 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
601 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
603 centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
605 {//the current edge has been subsplit -> create corresponding centers.
606 std::size_t nbOfCentersToAppend(edge1.size()/2);
607 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
608 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
609 std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
610 for(std::size_t k=0;k<nbOfCentersToAppend;k++)
613 const double *aa(coordsPtr+2*(*it3++));
614 const double *bb(coordsPtr+2*(*it3++));
615 ee->getMiddleOfPoints(aa,bb,tmpp);
616 addCoo->insertAtTheEnd(tmpp,tmpp+2);
617 centers.push_back(offset+k);
621 std::size_t nb(allEdges.size());
623 throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
624 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
625 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
627 ret->setCoords(coords);
630 addCoo->rearrange(2);
631 addCoo=DataArrayDouble::Aggregate(coords,addCoo);
632 ret->setCoords(addCoo);
634 ret->allocateCells(1);
635 std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
636 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
637 connOut[kk]=allEdges[2*kk];
638 connOut.insert(connOut.end(),centers.begin(),centers.end());
639 ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
644 * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
647 * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
649 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
651 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
652 if(!cm.isQuadratic())
653 return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
655 return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
658 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
661 for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
663 const INTERP_KERNEL::Edge *ee(*it);
664 if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
668 mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
671 const double *coo(mesh2D->getCoords()->begin());
672 std::size_t sz(conn.size());
673 std::vector<double> addCoo;
674 std::vector<int> conn2(conn);
675 int offset(mesh2D->getNumberOfNodes());
676 for(std::size_t i=0;i<sz;i++)
679 edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
680 addCoo.insert(addCoo.end(),tmp,tmp+2);
681 conn2.push_back(offset+(int)i);
683 mesh2D->getCoords()->rearrange(1);
684 mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
685 mesh2D->getCoords()->rearrange(2);
686 mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
691 * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
693 * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
694 * a set of edges defined in \a splitMesh1D.
696 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
697 std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
699 std::size_t nb(edge1Bis.size()/2);
700 std::size_t nbOfEdgesOf2DCellSplit(nb/2);
701 int iEnd(splitMesh1D->getNumberOfCells());
703 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
705 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
706 for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
707 for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
710 {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
711 out0.resize(1); out1.resize(1);
712 std::vector<int>& connOut(out0[0]);
713 connOut.resize(nbOfEdgesOf2DCellSplit);
714 std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
715 edgesPtr.resize(nbOfEdgesOf2DCellSplit);
716 for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
718 connOut[kk]=edge1Bis[2*kk];
719 edgesPtr[kk]=edge1BisPtr[2*kk];
724 // [i,iEnd[ contains the
725 out0.resize(2); out1.resize(2);
726 std::vector<int>& connOutLeft(out0[0]);
727 std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
728 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
729 std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
730 for(std::size_t k=ii;k<jj+1;k++)
731 { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
732 std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
733 for(int ik=0;ik<iEnd;ik++)
735 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
736 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
739 for(int ik=iEnd-1;ik>=0;ik--)
740 connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
741 for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
742 { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
743 eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
744 for(int ik=0;ik<iEnd;ik++)
745 connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
746 eright.insert(eright.end(),ees.begin(),ees.end());
754 CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
756 std::vector<int> _edges;
757 std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
760 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
762 std::size_t nbe(edges.size());
763 std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
764 for(std::size_t i=0;i<nbe;i++)
766 edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
767 edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
769 _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
770 std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
771 std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
777 EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
778 EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
779 bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
780 void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
781 void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
785 MCAuto<MEDCouplingUMesh> _mesh;
786 MCAuto<INTERP_KERNEL::Edge> _edge;
791 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
793 const MEDCouplingUMesh *mesh(_mesh);
799 { _left++; _right++; return ; }
802 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
803 if((isLeft && isRight) || (!isLeft && !isRight))
804 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
815 bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
816 if((isLeft && isRight) || (!isLeft && !isRight))
817 throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
832 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
834 const MEDCouplingUMesh *mesh(_mesh);
837 neighbors[0]=offset+_left; neighbors[1]=offset+_right;
840 {// not fully splitting cell case
841 if(mesh2D->getNumberOfCells()==1)
842 {//little optimization. 1 cell no need to find in which cell mesh is !
843 neighbors[0]=offset; neighbors[1]=offset;
848 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
849 int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
851 throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
852 neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
857 class VectorOfCellInfo
860 VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
861 std::size_t size() const { return _pool.size(); }
862 int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
863 void setMeshAt(int pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs);
864 const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
865 const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
866 MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
867 void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
869 int getZePosOfEdgeGivenItsGlobalId(int pos) const;
870 void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
871 const CellInfo& get(int pos) const;
872 CellInfo& get(int pos);
874 std::vector<CellInfo> _pool;
875 MCAuto<MEDCouplingUMesh> _ze_mesh;
876 std::vector<EdgeInfo> _edge_info;
879 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
881 _pool[0]._edges=edges;
882 _pool[0]._edges_ptr=edgesPtr;
885 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
888 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
891 const MEDCouplingUMesh *zeMesh(_ze_mesh);
893 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
894 MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
895 return zeMesh->getCellContainingPoint(barys->begin(),eps);
898 void VectorOfCellInfo::setMeshAt(int pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
900 get(pos);//to check pos
901 bool isFast(pos==0 && _pool.size()==1);
902 std::size_t sz(edges.size());
903 // dealing with edges
905 _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
907 _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
909 std::vector<CellInfo> pool(_pool.size()-1+sz);
910 for(int i=0;i<pos;i++)
912 for(std::size_t j=0;j<sz;j++)
913 pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
914 for(int i=pos+1;i<(int)_pool.size();i++)
915 pool[i+sz-1]=_pool[i];
919 updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
927 std::vector< MCAuto<MEDCouplingUMesh> > ms;
930 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
934 if(pos<_ze_mesh->getNumberOfCells()-1)
936 MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
939 std::vector< const MEDCouplingUMesh *> ms2(ms.size());
940 for(std::size_t j=0;j<ms2.size();j++)
942 _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
945 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
947 _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
950 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
953 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
955 for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
957 if((*it).isInMyRange(pos))
960 throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
963 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
966 if(_edge_info.empty())
968 std::size_t sz(_edge_info.size()-1);
969 for(std::size_t i=0;i<sz;i++)
970 _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
973 const CellInfo& VectorOfCellInfo::get(int pos) const
975 if(pos<0 || pos>=(int)_pool.size())
976 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
980 CellInfo& VectorOfCellInfo::get(int pos)
982 if(pos<0 || pos>=(int)_pool.size())
983 throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
989 * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
990 * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
992 * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
994 * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
996 * \param [in] allEdges a list of pairs (beginNode, endNode). Linked with \a allEdgesPtr to get the equation of edge.
998 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
999 MCAuto<DataArrayInt>& idsLeftRight)
1001 int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1002 if(nbCellsInSplitMesh1D==0)
1003 throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1004 const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1005 std::size_t nb(allEdges.size()),jj;
1007 throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1008 std::vector<int> edge1Bis(nb*2);
1009 std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1010 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1011 std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1012 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1013 std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1015 idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1016 int *idsLeftRightPtr(idsLeftRight->getPointer());
1017 VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1018 for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1019 {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1021 for(;iEnd<nbCellsInSplitMesh1D;)
1023 for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1029 if(iEnd<nbCellsInSplitMesh1D)
1032 MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1033 int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1035 MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1036 retTmp->setCoords(splitMesh1D->getCoords());
1037 retTmp->allocateCells();
1039 std::vector< std::vector<int> > out0;
1040 std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1042 BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1043 for(std::size_t cnt=0;cnt<out0.size();cnt++)
1044 AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1045 pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1049 for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1050 pool.feedEdgeInfoAt(eps,mm,offset,idsLeftRightPtr+2*mm);
1051 return pool.getZeMesh().retn();
1054 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, const MEDCouplingUMesh *splitMesh1D,
1055 const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1056 MCAuto<DataArrayInt>& idsLeftRight)
1058 const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1060 std::vector<int> allEdges;
1061 std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1062 for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1064 int edgeId(std::abs(*it)-1);
1065 std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1066 MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1067 const std::vector<int>& edge1(intersectEdge1[edgeId]);
1069 allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1071 allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1072 std::size_t sz(edge1.size());
1073 for(std::size_t cnt=0;cnt<sz;cnt++)
1074 allEdgesPtr.push_back(ee);
1077 return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1080 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1082 if(!typ1.isQuadratic() && !typ2.isQuadratic())
1083 {//easy case comparison not
1084 return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1086 else if(typ1.isQuadratic() && typ2.isQuadratic())
1088 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1091 if(conn1[2]==conn2[2])
1093 const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1094 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1098 {//only one is quadratic
1099 bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1102 const double *a(0),*bb(0),*be(0);
1103 if(typ1.isQuadratic())
1105 a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1109 a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1111 double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1112 double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1118 * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1119 * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1121 * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1123 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1125 if(candidatesIn2DEnd==candidatesIn2DBg)
1126 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1127 const double *coo(mesh2DSplit->getCoords()->begin());
1128 if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1129 return *candidatesIn2DBg;
1130 int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1131 MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1132 if(cellIdInMesh1DSplitRelative<0)
1133 cur1D->changeOrientationOfCells();
1134 const int *c1D(cur1D->getNodalConnectivity()->begin());
1135 const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1136 for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1138 MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1139 const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1140 const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1141 unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1142 INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1143 for(unsigned it2=0;it2<sz;it2++)
1145 INTERP_KERNEL::NormalizedCellType typeOfSon;
1146 cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1147 const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1148 if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1152 throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1156 * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of int given resp start and stop.
1157 * 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.
1158 * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1159 * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1160 * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1161 * \param [out] addCoo - nodes to be append at the end
1162 * \param [out] mergedNodes - gives all pair of nodes of \a m2Desc that have same location than some nodes in \a m1Desc. key is id in \a m2Desc offseted and value is id in \a m1Desc.
1164 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1165 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2, std::vector<double>& addCoo, std::map<int,int>& mergedNodes)
1167 static const int SPACEDIM=2;
1168 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1169 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1170 const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1171 // Build BB tree of all edges in the tool mesh (second mesh)
1172 MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree());
1173 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1174 int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1175 intersectEdge1.resize(nDescCell1);
1176 colinear2.resize(nDescCell2);
1177 subDiv2.resize(nDescCell2);
1178 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1180 std::vector<int> candidates1(1);
1181 int offset1(m1Desc->getNumberOfNodes());
1182 int offset2(offset1+m2Desc->getNumberOfNodes());
1183 for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
1185 std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1186 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1187 if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1189 std::map<INTERP_KERNEL::Node *,int> map1,map2;
1190 // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1191 INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1193 INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1194 // 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
1195 // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1196 std::set<INTERP_KERNEL::Node *> nodes;
1197 pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1198 std::size_t szz(nodes.size());
1199 std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1200 std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1201 for(std::size_t iii=0;iii<szz;iii++,itt++)
1202 { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1203 // end of protection
1204 // Performs egde cutting:
1205 pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1210 // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1211 intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1217 * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1218 * It builds the descending connectivity of the two meshes, and then using a binary tree
1219 * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1220 * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1222 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1223 std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1224 MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1225 std::vector<double>& addCoo,
1226 MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1228 // Build desc connectivity
1229 desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1230 desc2=DataArrayInt::New();
1231 descIndx2=DataArrayInt::New();
1232 revDesc2=DataArrayInt::New();
1233 revDescIndx2=DataArrayInt::New();
1234 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1235 MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1236 m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1237 m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1238 MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1239 std::map<int,int> notUsedMap;
1240 Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1241 m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1242 m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1246 * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1247 * (newly created) nodes corresponding to the edge intersections.
1249 * @param[out] cr, crI connectivity of the resulting mesh
1250 * @param[out] cNb1, cNb2 correspondance arrays giving for the merged mesh the initial cells IDs in m1 / m2
1251 * TODO: describe input parameters
1253 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1254 const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1255 const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1256 const std::vector<double>& addCoords,
1257 std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1259 static const int SPACEDIM=2;
1260 const double *coo1(m1->getCoords()->begin());
1261 const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1262 int offset1(m1->getNumberOfNodes());
1263 const double *coo2(m2->getCoords()->begin());
1264 const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1265 int offset2(offset1+m2->getNumberOfNodes());
1266 int offset3(offset2+((int)addCoords.size())/2);
1267 MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
1268 const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1269 // Here a BBTree on 2D-cells, not on segments:
1270 BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1271 int ncell1(m1->getNumberOfCells());
1273 for(int i=0;i<ncell1;i++)
1275 std::vector<int> candidates2;
1276 myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1277 std::map<INTERP_KERNEL::Node *,int> mapp;
1278 std::map<int,INTERP_KERNEL::Node *> mappRev;
1279 INTERP_KERNEL::QuadraticPolygon pol1;
1280 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1281 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1282 // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1283 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1284 // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
1285 pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1286 desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1288 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
1289 std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1290 INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1291 for(it1.first();!it1.finished();it1.next())
1292 edges1.insert(it1.current()->getPtr());
1294 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1295 std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1297 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1299 INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1300 const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1301 // Complete mapping with elements coming from the current cell it2 in mesh2:
1302 MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1303 // pol2 is the new QP in the final merged result.
1304 pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1305 pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1308 for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1310 INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1311 pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1312 //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1313 pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1315 // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1316 // by m2 but that we still want to keep in the final result.
1321 INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1323 catch(INTERP_KERNEL::Exception& e)
1325 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();
1326 throw INTERP_KERNEL::Exception(oss.str());
1329 for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1330 (*it).second->decrRef();
1334 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1336 std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1339 std::size_t sz(conn.size());
1340 std::size_t found(std::numeric_limits<std::size_t>::max());
1341 for(std::size_t i=0;i<sz;i++)
1343 int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1344 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]};
1345 double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1346 std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1347 std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1349 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];
1350 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]);
1352 if(dotTest>eps && dotTest<1.-eps)
1358 if(found==std::numeric_limits<std::size_t>::max())
1359 throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1360 conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1363 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1365 std::size_t sz(conn.size());
1366 std::vector<int> *curPart(&part0);
1367 for(std::size_t i=0;i<sz;i++)
1369 int nextt(conn[(i+1)%sz]);
1370 (*curPart).push_back(nextt);
1371 if(nextt==pt0 || nextt==pt1)
1377 (*curPart).push_back(nextt);
1383 * 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.
1385 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1386 const int *desc, const int *descIndx, const double *coords, double eps,
1387 std::vector<std::vector<int> >& res) const
1389 checkFullyDefined();
1390 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1391 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1392 const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1393 int nbOfCells(getNumberOfCells());
1395 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1396 for(int i=0;i<nbOfCells;i++)
1398 int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset),start(-1),end(-1);
1399 for(int j=0;j<nbOfFaces;j++)
1401 const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1402 const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1403 int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1404 INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1405 INTERP_KERNEL::NormalizedCellType cmsId;
1406 unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1407 std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1408 if(p.first!=-1 && p.second!=-1)
1412 InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1413 InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1414 std::vector<int> elt1,elt2;
1415 SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1416 res.push_back(elt1);
1417 res.push_back(elt2);
1429 * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
1431 * \sa MEDCouplingUMesh::split2DCells
1433 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1435 checkConnectivityFullyDefined();
1436 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1437 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1438 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1439 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1440 int prevPosOfCi(ciPtr[0]);
1441 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1443 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1444 *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1445 for(int j=0;j<sz;j++)
1447 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1448 for(int k=0;k<sz2;k++)
1449 *cPtr++=subPtr[offset2+k];
1451 *cPtr++=oldConn[prevPosOfCi+j+2];
1454 prevPosOfCi=ciPtr[1];
1455 ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1458 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1459 _nodal_connec->decrRef();
1460 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1465 * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
1467 * \return int - the number of new nodes created.
1468 * \sa MEDCouplingUMesh::split2DCells
1470 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1472 checkConsistencyLight();
1473 int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1474 MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1475 MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1476 const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1477 const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1478 const double *oldCoordsPtr(getCoords()->begin());
1479 int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1480 int prevPosOfCi(ciPtr[0]);
1481 for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1483 int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1484 for(int j=0;j<sz;j++)
1485 { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1486 *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1487 for(int j=0;j<sz;j++)//loop over subedges of oldConn
1489 int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1493 cPtr[1]=oldConn[prevPosOfCi+2+j];
1494 cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1497 std::vector<INTERP_KERNEL::Node *> ns(3);
1498 ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1499 ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1500 ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1501 MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1502 for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1504 cPtr[1]=subPtr[offset2+k];
1505 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1507 int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1510 cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1512 prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1513 ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1516 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1517 _nodal_connec->decrRef();
1518 _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1519 addCoo->rearrange(2);
1520 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1522 return addCoo->getNumberOfTuples();
1529 * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1530 * returns a result mesh constituted by polygons.
1531 * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1532 * all nodes from m2.
1533 * The meshes should be in 2D space. In
1534 * addition, returns two arrays mapping cells of the result mesh to cells of the input
1536 * \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
1537 * 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)
1538 * \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
1539 * 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)
1540 * \param [in] eps - precision used to detect coincident mesh entities.
1541 * \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1542 * cell an id of the cell of \a m1 it comes from. The caller is to delete
1543 * this array using decrRef() as it is no more needed.
1544 * \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1545 * cell an id of the cell of \a m2 it comes from. -1 value means that a
1546 * result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1547 * any cell of \a m2. The caller is to delete this array using decrRef() as
1548 * it is no more needed.
1549 * \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1550 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1551 * is no more needed.
1552 * \throw If the coordinates array is not set in any of the meshes.
1553 * \throw If the nodal connectivity of cells is not defined in any of the meshes.
1554 * \throw If any of the meshes is not a 2D mesh in 2D space.
1556 * \sa conformize2D, mergeNodes
1558 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1559 double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1562 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1563 m1->checkFullyDefined();
1564 m2->checkFullyDefined();
1565 if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1566 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!");
1568 // Step 1: compute all edge intersections (new nodes)
1569 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1570 MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1571 DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1572 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1573 IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1574 m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1575 addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1576 revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1577 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1578 MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1580 // Step 2: re-order newly created nodes according to the ordering found in m2
1581 std::vector< std::vector<int> > intersectEdge2;
1582 BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1583 subDiv2.clear(); dd5=0; dd6=0;
1586 std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1587 std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1588 BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1589 /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1591 // Step 4: Prepare final result:
1592 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1593 addCooDa->alloc((int)(addCoo.size())/2,2);
1594 std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1595 MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1596 addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1597 std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1598 std::vector<const DataArrayDouble *> coordss(4);
1599 coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1600 MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1601 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1602 MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1603 MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1604 MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1605 MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1606 ret->setConnectivity(conn,connI,true);
1607 ret->setCoords(coo);
1608 cellNb1=c1.retn(); cellNb2=c2.retn();
1613 * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1614 * 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
1615 * and finaly, in case of quadratic polygon the centers of edges new nodes.
1616 * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1618 * \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
1619 * 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)
1620 * \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
1621 * you can invoke orderConsecutiveCells1D on \a mesh1D.
1622 * \param [in] eps - precision used to perform intersections and localization operations.
1623 * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1624 * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1625 * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1626 * 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.
1627 * \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
1628 * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1629 * 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.
1631 * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1633 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1635 if(!mesh2D || !mesh1D)
1636 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1637 mesh2D->checkFullyDefined();
1638 mesh1D->checkFullyDefined();
1639 const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1640 if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1641 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1642 // Step 1: compute all edge intersections (new nodes)
1643 std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1644 std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
1645 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1646 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1648 // Build desc connectivity
1649 DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1650 MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1651 MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1652 std::map<int,int> mergedNodes;
1653 Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1654 // use mergeNodes to fix intersectEdge1
1655 for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1657 std::size_t n((*it0).size()/2);
1658 int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1659 std::map<int,int>::const_iterator it1;
1660 it1=mergedNodes.find(eltStart);
1661 if(it1!=mergedNodes.end())
1662 (*it0)[0]=(*it1).second;
1663 it1=mergedNodes.find(eltEnd);
1664 if(it1!=mergedNodes.end())
1665 (*it0)[2*n-1]=(*it1).second;
1668 MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1669 addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
1670 // Step 2: re-order newly created nodes according to the ordering found in m2
1671 std::vector< std::vector<int> > intersectEdge2;
1672 BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1674 // Step 3: compute splitMesh1D
1675 MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1676 MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1677 MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1678 idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1679 MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1680 MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1681 // deal with cells in mesh2D that are not cut but only some of their edges are
1682 MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1683 idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1684 idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1685 MCAuto<DataArrayInt> out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells
1686 if(!idsInDesc2DToBeRefined->empty())
1688 DataArrayInt *out0(0),*outi0(0);
1689 MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1690 MCAuto<DataArrayInt> outi0s(outi0);
1692 out0s=out0s->buildUnique();
1696 MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1697 MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1698 MCAuto<DataArrayInt> elts,eltsIndex;
1699 mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1700 MCAuto<DataArrayInt> eltsIndex2(eltsIndex->deltaShiftIndex());
1701 MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1702 if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1703 throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1704 MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1705 MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1706 if((DataArrayInt *)out0s)
1707 untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1708 std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1709 // OK all is ready to insert in ret2 mesh
1710 if(!untouchedCells->empty())
1711 {// the most easy part, cells in mesh2D not impacted at all
1712 outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1713 outMesh2DSplit.back()->setCoords(ret1->getCoords());
1714 ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1716 if((DataArrayInt *)out0s)
1717 {// here dealing with cells in out0s but not in cellsToBeModified
1718 MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1719 const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1720 for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1722 outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1723 ret1->setCoords(outMesh2DSplit.back()->getCoords());
1725 int offset(ret2->getNumberOfTuples());
1726 ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1727 MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1728 partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1729 int kk(0),*ret3ptr(partOfRet3->getPointer());
1730 for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1732 int faceId(std::abs(*it)-1);
1733 for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1735 int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1738 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1739 ret3ptr[2*kk]=tmp+offset;
1740 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1741 ret3ptr[2*kk+1]=tmp+offset;
1744 {//the current edge is shared by a 2D cell that will be split just after
1745 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1746 ret3ptr[2*kk]=-(*it2+1);
1747 if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1748 ret3ptr[2*kk+1]=-(*it2+1);
1752 m1Desc->setCoords(ret1->getCoords());
1753 ret1NonCol->setCoords(ret1->getCoords());
1754 ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1755 if(!outMesh2DSplit.empty())
1757 DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1758 for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1759 (*itt)->setCoords(da);
1762 cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1763 for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1765 MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1766 idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1767 MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1768 MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1769 MCAuto<DataArrayInt> partOfRet3;
1770 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));
1771 ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1772 outMesh2DSplit.push_back(splitOfOneCell);
1773 for(int i=0;i<splitOfOneCell->getNumberOfCells();i++)
1774 ret2->pushBackSilent(*it);
1777 std::size_t nbOfMeshes(outMesh2DSplit.size());
1778 std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1779 for(std::size_t i=0;i<nbOfMeshes;i++)
1780 tmp[i]=outMesh2DSplit[i];
1782 ret1->getCoords()->setInfoOnComponents(compNames);
1783 MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1784 // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1786 MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStricltyNegative());
1787 for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1789 int old2DCellId(-ret3->getIJ(*it,0)-1);
1790 MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1791 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
1793 ret3->changeValue(std::numeric_limits<int>::max(),-1);
1796 splitMesh1D=ret1.retn();
1797 splitMesh2D=ret2D.retn();
1798 cellIdInMesh2D=ret2.retn();
1799 cellIdInMesh1D=ret3.retn();
1803 * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1804 * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1805 * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this this method
1806 * 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).
1807 * 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 invokation of this method.
1809 * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1810 * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1812 * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1813 * This method expects that all nodes in \a this are not closer than \a eps.
1814 * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1816 * \param [in] eps the relative error to detect merged edges.
1817 * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array
1818 * that the user is expected to deal with.
1820 * \throw If \a this is not coherent.
1821 * \throw If \a this has not spaceDim equal to 2.
1822 * \throw If \a this has not meshDim equal to 2.
1823 * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1825 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1827 static const int SPACEDIM=2;
1828 checkConsistencyLight();
1829 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1830 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1831 MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1832 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1833 const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1834 MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
1835 const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1836 int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1837 std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1838 std::vector<double> addCoo;
1839 BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1840 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1841 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1842 for(int i=0;i<nDescCell;i++)
1844 std::vector<int> candidates;
1845 myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1846 for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1849 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1850 INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1851 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1852 INTERP_KERNEL::MergePoints merge;
1853 INTERP_KERNEL::QuadraticPolygon c1,c2;
1854 e1->intersectWith(e2,merge,c1,c2);
1855 e1->decrRef(); e2->decrRef();
1856 if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1857 overlapEdge[i].push_back(*it);
1858 if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1859 overlapEdge[*it].push_back(i);
1862 // splitting done. sort intersect point in intersectEdge.
1863 std::vector< std::vector<int> > middle(nDescCell);
1864 int nbOf2DCellsToBeSplit(0);
1865 bool middleNeedsToBeUsed(false);
1866 std::vector<bool> cells2DToTreat(nDescCell,false);
1867 for(int i=0;i<nDescCell;i++)
1869 std::vector<int>& isect(intersectEdge[i]);
1870 int sz((int)isect.size());
1873 std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1874 INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1875 e->sortSubNodesAbs(coords,isect);
1880 int idx0(rdi[i]),idx1(rdi[i+1]);
1882 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1883 if(!cells2DToTreat[rd[idx0]])
1885 cells2DToTreat[rd[idx0]]=true;
1886 nbOf2DCellsToBeSplit++;
1888 // try to reuse at most eventual 'middle' of SEG3
1889 std::vector<int>& mid(middle[i]);
1890 mid.resize(sz+1,-1);
1891 if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1893 middleNeedsToBeUsed=true;
1894 const std::vector<int>& candidates(overlapEdge[i]);
1895 std::vector<int> trueCandidates;
1896 for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1897 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1898 trueCandidates.push_back(*itc);
1899 int stNode(c[ci[i]+1]),endNode(isect[0]);
1900 for(int j=0;j<sz+1;j++)
1902 for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1904 int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1905 if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1906 { mid[j]=*itc; break; }
1909 endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1914 MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1915 if(nbOf2DCellsToBeSplit==0)
1918 int *retPtr(ret->getPointer());
1919 for(int i=0;i<nCell;i++)
1920 if(cells2DToTreat[i])
1923 MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1924 DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1925 MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1926 DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1927 if(middleNeedsToBeUsed)
1928 { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1929 MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1930 int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1931 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.
1932 setPartOfMySelf(ret->begin(),ret->end(),*modif);
1934 bool areNodesMerged; int newNbOfNodes;
1935 if(nbOfNodesCreated!=0)
1936 MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
1942 * 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.
1943 * 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).
1944 * 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
1945 * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
1946 * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automaticaly. For 2D "repaired" quadratic cells
1947 * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
1949 * 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
1950 * using new instance, idem for coordinates.
1952 * 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.
1954 * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
1956 * \throw If \a this is not coherent.
1957 * \throw If \a this has not spaceDim equal to 2.
1958 * \throw If \a this has not meshDim equal to 2.
1960 * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
1962 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
1964 MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1965 checkConsistencyLight();
1966 if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1967 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1968 INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1969 INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1970 int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
1971 const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
1972 MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
1973 MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
1974 const double *coords(_coords->begin());
1975 int *newciptr(newci->getPointer());
1976 for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
1978 if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
1979 ret->pushBackSilent(i);
1980 newciptr[1]=newc->getNumberOfTuples();
1985 if(!appendedCoords->empty())
1987 appendedCoords->rearrange(2);
1988 MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
1990 setCoords(newCoords);
1993 setConnectivity(newc,newci,true);