1 // Copyright (C) 2007-2020 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingUMesh.hxx"
22 #include "MEDCouplingCMesh.hxx"
23 #include "MEDCoupling1GTUMesh.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingSkyLineArray.hxx"
26 #include "CellModel.hxx"
27 #include "VolSurfUser.txx"
28 #include "InterpolationUtils.hxx"
29 #include "PointLocatorAlgos.txx"
31 #include "BBTreeDst.txx"
32 #include "SplitterTetra.hxx"
33 #include "DiameterCalculator.hxx"
34 #include "DirectedBoundingBox.hxx"
35 #include "InterpKernelMatrixTools.hxx"
36 #include "InterpKernelMeshQuality.hxx"
37 #include "InterpKernelCellSimplify.hxx"
38 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
39 #include "InterpKernelAutoPtr.hxx"
40 #include "InterpKernelGeo2DNode.hxx"
41 #include "InterpKernelGeo2DEdgeLin.hxx"
42 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
43 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
44 #include "MEDCouplingUMesh_internal.hxx"
53 using namespace MEDCoupling;
56 * This method checks that all arrays are set. If yes nothing done if no an exception is thrown.
58 void MEDCouplingUMesh::checkFullyDefined() const
60 if(!_nodal_connec_index || !_nodal_connec || !_coords)
61 throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity and coordinates set in unstructured mesh.");
65 * This method checks that all connectivity arrays are set. If yes nothing done if no an exception is thrown.
67 void MEDCouplingUMesh::checkConnectivityFullyDefined() const
69 if(!_nodal_connec_index || !_nodal_connec)
70 throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity set in unstructured mesh.");
73 void MEDCouplingUMesh::reprConnectivityOfThisLL(std::ostringstream& stream) const
75 if(_nodal_connec!=0 && _nodal_connec_index!=0)
77 mcIdType nbOfCells=getNumberOfCells();
78 const mcIdType *c=_nodal_connec->getConstPointer();
79 const mcIdType *ci=_nodal_connec_index->getConstPointer();
80 for(mcIdType i=0;i<nbOfCells;i++)
82 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[i]]);
83 stream << "Cell #" << i << " " << cm.getRepr() << " : ";
84 std::copy(c+ci[i]+1,c+ci[i+1],std::ostream_iterator<int>(stream," "));
89 stream << "Connectivity not defined !\n";
94 * This method implements policy 0 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
96 DataArrayIdType *MEDCouplingUMesh::simplexizePol0()
98 checkConnectivityFullyDefined();
99 if(getMeshDimension()!=2)
100 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
101 mcIdType nbOfCells=getNumberOfCells();
102 MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
103 mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
104 ret->alloc(nbOfCells+nbOfCutCells,1);
105 if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
106 mcIdType *retPt=ret->getPointer();
107 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
108 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
109 newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
110 newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
111 mcIdType *pt=newConn->getPointer();
112 mcIdType *ptI=newConnI->getPointer();
114 const mcIdType *oldc=_nodal_connec->begin();
115 const mcIdType *ci=_nodal_connec_index->begin();
116 for(mcIdType i=0;i<nbOfCells;i++,ci++)
118 if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
120 const mcIdType tmp[8]={ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+3],
121 ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+1],oldc[ci[0]+3],oldc[ci[0]+4]};
122 pt=std::copy(tmp,tmp+8,pt);
131 pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
132 ptI[1]=ptI[0]+ci[1]-ci[0];
137 _nodal_connec->decrRef();
138 _nodal_connec=newConn.retn();
139 _nodal_connec_index->decrRef();
140 _nodal_connec_index=newConnI.retn();
147 * This method implements policy 1 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
149 DataArrayIdType *MEDCouplingUMesh::simplexizePol1()
151 checkConnectivityFullyDefined();
152 if(getMeshDimension()!=2)
153 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
154 mcIdType nbOfCells=getNumberOfCells();
155 MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
156 mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
157 ret->alloc(nbOfCells+nbOfCutCells,1);
158 if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
159 mcIdType *retPt=ret->getPointer();
160 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
161 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
162 newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
163 newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
164 mcIdType *pt=newConn->getPointer();
165 mcIdType *ptI=newConnI->getPointer();
167 const mcIdType *oldc=_nodal_connec->begin();
168 const mcIdType *ci=_nodal_connec_index->begin();
169 for(mcIdType i=0;i<nbOfCells;i++,ci++)
171 if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
173 const mcIdType tmp[8]={ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+4],
174 ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+2],oldc[ci[0]+3],oldc[ci[0]+4]};
175 pt=std::copy(tmp,tmp+8,pt);
184 pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
185 ptI[1]=ptI[0]+ci[1]-ci[0];
190 _nodal_connec->decrRef();
191 _nodal_connec=newConn.retn();
192 _nodal_connec_index->decrRef();
193 _nodal_connec_index=newConnI.retn();
200 * This method implements policy INTERP_KERNEL::PLANAR_FACE_5 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
202 DataArrayIdType *MEDCouplingUMesh::simplexizePlanarFace5()
204 checkConnectivityFullyDefined();
205 if(getMeshDimension()!=3)
206 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePlanarFace5 : this policy is only available for mesh with meshdim == 3 !");
207 mcIdType nbOfCells=getNumberOfCells();
208 MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
209 mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
210 ret->alloc(nbOfCells+4*nbOfCutCells,1);
211 if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
212 mcIdType *retPt=ret->getPointer();
213 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
214 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
215 newConnI->alloc(nbOfCells+4*nbOfCutCells+1,1);
216 newConn->alloc(getNodalConnectivityArrayLen()+16*nbOfCutCells,1);//21
217 mcIdType *pt=newConn->getPointer();
218 mcIdType *ptI=newConnI->getPointer();
220 const mcIdType *oldc=_nodal_connec->begin();
221 const mcIdType *ci=_nodal_connec_index->begin();
222 for(mcIdType i=0;i<nbOfCells;i++,ci++)
224 if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_HEXA8)
226 for(int j=0;j<5;j++,pt+=5,ptI++)
228 pt[0]=ToIdType(INTERP_KERNEL::NORM_TETRA4);
229 pt[1]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+0]+1]; pt[2]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+1]+1]; pt[3]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+2]+1]; pt[4]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+3]+1];
236 pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
237 ptI[1]=ptI[0]+ci[1]-ci[0];
242 _nodal_connec->decrRef();
243 _nodal_connec=newConn.retn();
244 _nodal_connec_index->decrRef();
245 _nodal_connec_index=newConnI.retn();
252 * This method implements policy INTERP_KERNEL::PLANAR_FACE_6 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
254 DataArrayIdType *MEDCouplingUMesh::simplexizePlanarFace6()
256 checkConnectivityFullyDefined();
257 if(getMeshDimension()!=3)
258 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePlanarFace6 : this policy is only available for mesh with meshdim == 3 !");
259 mcIdType nbOfCells=getNumberOfCells();
260 MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
261 mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
262 ret->alloc(nbOfCells+5*nbOfCutCells,1);
263 if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
264 mcIdType *retPt=ret->getPointer();
265 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
266 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
267 newConnI->alloc(nbOfCells+5*nbOfCutCells+1,1);
268 newConn->alloc(getNodalConnectivityArrayLen()+21*nbOfCutCells,1);
269 mcIdType *pt=newConn->getPointer();
270 mcIdType *ptI=newConnI->getPointer();
272 const mcIdType *oldc=_nodal_connec->begin();
273 const mcIdType *ci=_nodal_connec_index->begin();
274 for(mcIdType i=0;i<nbOfCells;i++,ci++)
276 if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_HEXA8)
278 for(int j=0;j<6;j++,pt+=5,ptI++)
280 pt[0]=ToIdType(INTERP_KERNEL::NORM_TETRA4);
281 pt[1]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+0]+1]; pt[2]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+1]+1]; pt[3]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+2]+1]; pt[4]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+3]+1];
288 pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
289 ptI[1]=ptI[0]+ci[1]-ci[0];
294 _nodal_connec->decrRef();
295 _nodal_connec=newConn.retn();
296 _nodal_connec_index->decrRef();
297 _nodal_connec_index=newConnI.retn();
304 * Tessellates \a this 2D mesh by dividing not straight edges of quadratic faces,
305 * so that the number of cells remains the same. Quadratic faces are converted to
306 * polygons. This method works only for 2D meshes in
307 * 2D space. If no cells are quadratic (INTERP_KERNEL::NORM_QUAD8,
308 * INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QPOLYG ), \a this mesh remains unchanged.
309 * \warning This method can lead to a huge amount of nodes if \a eps is very low.
310 * \param [in] eps - specifies the maximal angle (in radians) between 2 sub-edges of
311 * a polylinized edge constituting the input polygon.
312 * \throw If the coordinates array is not set.
313 * \throw If the nodal connectivity of cells is not defined.
314 * \throw If \a this->getMeshDimension() != 2.
315 * \throw If \a this->getSpaceDimension() != 2.
317 void MEDCouplingUMesh::tessellate2DInternal(double eps)
320 if(getMeshDimension()!=2 || getSpaceDimension()!=2)
321 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DInternal works on umeshes with meshdim equal to 2 and spaceDim equal to 2 too!");
322 double epsa=fabs(eps);
323 if(epsa<std::numeric_limits<double>::min())
324 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DInternal : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
325 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),descIndx1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescIndx1(DataArrayIdType::New());
326 MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
327 revDesc1=0; revDescIndx1=0;
328 mDesc->tessellate2D(eps);
329 subDivide2DMesh(mDesc->_nodal_connec->begin(),mDesc->_nodal_connec_index->begin(),desc1->begin(),descIndx1->begin());
330 setCoords(mDesc->getCoords());
334 * Tessellates \a this 1D mesh in 2D space by dividing not straight quadratic edges.
335 * \warning This method can lead to a huge amount of nodes if \a eps is very low.
336 * \param [in] eps - specifies the maximal angle (in radian) between 2 sub-edges of
337 * a sub-divided edge.
338 * \throw If the coordinates array is not set.
339 * \throw If the nodal connectivity of cells is not defined.
340 * \throw If \a this->getMeshDimension() != 1.
341 * \throw If \a this->getSpaceDimension() != 2.
343 void MEDCouplingUMesh::tessellate2DCurveInternal(double eps)
346 if(getMeshDimension()!=1 || getSpaceDimension()!=2)
347 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurveInternal works on umeshes with meshdim equal to 1 and spaceDim equal to 2 too!");
348 double epsa=fabs(eps);
349 if(epsa<std::numeric_limits<double>::min())
350 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurveInternal : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
351 INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1.e-10); // RAII
352 mcIdType nbCells=getNumberOfCells();
353 mcIdType nbNodes=getNumberOfNodes();
354 const mcIdType *conn=_nodal_connec->begin();
355 const mcIdType *connI=_nodal_connec_index->begin();
356 const double *coords=_coords->begin();
357 std::vector<double> addCoo;
358 std::vector<mcIdType> newConn;//no direct DataArrayIdType because interface with Geometric2D
359 MCAuto<DataArrayIdType> newConnI(DataArrayIdType::New());
360 newConnI->alloc(nbCells+1,1);
361 mcIdType *newConnIPtr=newConnI->getPointer();
364 INTERP_KERNEL::Node *tmp2[3];
365 std::set<INTERP_KERNEL::NormalizedCellType> types;
366 for(mcIdType i=0;i<nbCells;i++,newConnIPtr++)
368 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
370 {//assert(connI[i+1]-connI[i]-1==3)
371 tmp1[0]=conn[connI[i]+1+0]; tmp1[1]=conn[connI[i]+1+1]; tmp1[2]=conn[connI[i]+1+2];
372 tmp2[0]=new INTERP_KERNEL::Node(coords[2*tmp1[0]],coords[2*tmp1[0]+1]);
373 tmp2[1]=new INTERP_KERNEL::Node(coords[2*tmp1[1]],coords[2*tmp1[1]+1]);
374 tmp2[2]=new INTERP_KERNEL::Node(coords[2*tmp1[2]],coords[2*tmp1[2]+1]);
375 INTERP_KERNEL::EdgeArcCircle *eac=INTERP_KERNEL::EdgeArcCircle::BuildFromNodes(tmp2[0],tmp2[2],tmp2[1]);
378 eac->tesselate(tmp1,nbNodes,epsa,newConn,addCoo);
379 types.insert((INTERP_KERNEL::NormalizedCellType)newConn[newConnIPtr[0]]);
381 newConnIPtr[1]=ToIdType(newConn.size());
385 types.insert(INTERP_KERNEL::NORM_SEG2);
386 newConn.push_back(INTERP_KERNEL::NORM_SEG2);
387 newConn.insert(newConn.end(),conn+connI[i]+1,conn+connI[i]+3);
388 newConnIPtr[1]=newConnIPtr[0]+3;
393 types.insert((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
394 newConn.insert(newConn.end(),conn+connI[i],conn+connI[i+1]);
395 newConnIPtr[1]=newConnIPtr[0]+3;
398 if(addCoo.empty() && ToIdType(newConn.size())==_nodal_connec->getNumberOfTuples())//nothing happens during tessellation : no update needed
401 DataArrayIdType::SetArrayIn(newConnI,_nodal_connec_index);
402 MCAuto<DataArrayIdType> newConnArr=DataArrayIdType::New();
403 newConnArr->alloc(newConn.size(),1);
404 std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer());
405 DataArrayIdType::SetArrayIn(newConnArr,_nodal_connec);
406 MCAuto<DataArrayDouble> newCoords=DataArrayDouble::New();
407 newCoords->alloc(nbNodes+addCoo.size()/2,2);
408 double *work=std::copy(_coords->begin(),_coords->end(),newCoords->getPointer());
409 std::copy(addCoo.begin(),addCoo.end(),work);
410 DataArrayDouble::SetArrayIn(newCoords,_coords);
416 * This private method is used to subdivide edges of a mesh with meshdim==2. If \a this has no a meshdim equal to 2 an exception will be thrown.
417 * This method completely ignores coordinates.
418 * \param nodeSubdived is the nodal connectivity of subdivision of edges
419 * \param nodeIndxSubdived is the nodal connectivity index of subdivision of edges
420 * \param desc is descending connectivity in format specified in MEDCouplingUMesh::buildDescendingConnectivity2
421 * \param descIndex is descending connectivity index in format specified in MEDCouplingUMesh::buildDescendingConnectivity2
423 void MEDCouplingUMesh::subDivide2DMesh(const mcIdType *nodeSubdived, const mcIdType *nodeIndxSubdived, const mcIdType *desc, const mcIdType *descIndex)
426 if(getMeshDimension()!=2)
427 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : works only on umesh with meshdim==2 !");
428 mcIdType nbOfCells=getNumberOfCells();
429 mcIdType *connI=_nodal_connec_index->getPointer();
430 mcIdType newConnLgth=0;
431 for(mcIdType i=0;i<nbOfCells;i++,connI++)
433 mcIdType offset=descIndex[i];
434 mcIdType nbOfEdges=descIndex[i+1]-offset;
436 bool ddirect=desc[offset+nbOfEdges-1]>0;
437 mcIdType eedgeId=std::abs(desc[offset+nbOfEdges-1])-1;
438 mcIdType ref=ddirect?nodeSubdived[nodeIndxSubdived[eedgeId+1]-1]:nodeSubdived[nodeIndxSubdived[eedgeId]+1];
439 for(mcIdType j=0;j<nbOfEdges;j++)
441 bool direct=desc[offset+j]>0;
442 mcIdType edgeId=std::abs(desc[offset+j])-1;
443 if(!INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodeSubdived[nodeIndxSubdived[edgeId]]).isQuadratic())
445 mcIdType id1=nodeSubdived[nodeIndxSubdived[edgeId]+1];
446 mcIdType id2=nodeSubdived[nodeIndxSubdived[edgeId+1]-1];
447 mcIdType ref2=direct?id1:id2;
450 mcIdType nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
451 newConnLgth+=nbOfSubNodes-1;
456 std::ostringstream oss; oss << "MEDCouplingUMesh::subDivide2DMesh : On polygon #" << i << " edgeid #" << j << " subedges mismatch : end subedge k!=start subedge k+1 !";
457 throw INTERP_KERNEL::Exception(oss.str());
462 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : this method only subdivides into linear edges !");
465 newConnLgth++;//+1 is for cell type
466 connI[1]=newConnLgth;
469 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
470 newConn->alloc(newConnLgth,1);
471 mcIdType *work=newConn->getPointer();
472 for(mcIdType i=0;i<nbOfCells;i++)
474 *work++=INTERP_KERNEL::NORM_POLYGON;
475 mcIdType offset=descIndex[i];
476 mcIdType nbOfEdges=descIndex[i+1]-offset;
477 for(mcIdType j=0;j<nbOfEdges;j++)
479 bool direct=desc[offset+j]>0;
480 mcIdType edgeId=std::abs(desc[offset+j])-1;
482 work=std::copy(nodeSubdived+nodeIndxSubdived[edgeId]+1,nodeSubdived+nodeIndxSubdived[edgeId+1]-1,work);
485 mcIdType nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
486 std::reverse_iterator<const mcIdType *> it(nodeSubdived+nodeIndxSubdived[edgeId+1]);
487 work=std::copy(it,it+nbOfSubNodes-1,work);
491 DataArrayIdType::SetArrayIn(newConn,_nodal_connec);
494 _types.insert(INTERP_KERNEL::NORM_POLYGON);
498 * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
499 * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
500 * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
501 * If \a fullyIn is true only cells whose ids are \b fully contained in [ \a begin,\a end ) tab will be kept.
503 * \param [in] begin input start of array of node ids.
504 * \param [in] end input end of array of node ids.
505 * \param [in] fullyIn input that specifies if all node ids must be in [ \a begin,\a end ) array to consider cell to be in.
506 * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
508 void MEDCouplingUMesh::fillCellIdsToKeepFromNodeIds(const mcIdType *begin, const mcIdType *end, bool fullyIn, DataArrayIdType *&cellIdsKeptArr) const
510 MCAuto<DataArrayIdType> cellIdsKept=DataArrayIdType::New(); cellIdsKept->alloc(0,1);
511 checkConnectivityFullyDefined();
513 mcIdType sz=getNodalConnectivity()->getMaxValue(tmp); sz=std::max(sz,ToIdType(0))+1;
514 std::vector<bool> fastFinder(sz,false);
515 for(const mcIdType *work=begin;work!=end;work++)
516 if(*work>=0 && *work<sz)
517 fastFinder[*work]=true;
518 mcIdType nbOfCells=getNumberOfCells();
519 const mcIdType *conn=getNodalConnectivity()->getConstPointer();
520 const mcIdType *connIndex=getNodalConnectivityIndex()->getConstPointer();
521 for(mcIdType i=0;i<nbOfCells;i++)
523 mcIdType ref=0,nbOfHit=0;
524 for(const mcIdType *work2=conn+connIndex[i]+1;work2!=conn+connIndex[i+1];work2++)
528 if(fastFinder[*work2])
531 if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
532 cellIdsKept->pushBackSilent(i);
534 cellIdsKeptArr=cellIdsKept.retn();
538 * This method works on a 3D curve linear mesh that is to say (meshDim==1 and spaceDim==3).
539 * If it is not the case an exception will be thrown.
540 * This method is non const because the coordinate of \a this can be appended with some new points issued from
541 * intersection of plane defined by ('origin','vec').
542 * This method has one in/out parameter : 'cut3DCurve'.
543 * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
544 * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
545 * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
546 * This method will throw an exception if \a this contains a non linear segment.
548 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<mcIdType>& cut3DCurve)
551 if(getMeshDimension()!=1 || getSpaceDimension()!=3)
552 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
553 mcIdType ncells=getNumberOfCells();
554 mcIdType nnodes=getNumberOfNodes();
555 double vec2[3],vec3[3],vec4[3];
556 double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
558 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
559 vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
560 const mcIdType *conn=_nodal_connec->getConstPointer();
561 const mcIdType *connI=_nodal_connec_index->getConstPointer();
562 const double *coo=_coords->getConstPointer();
563 std::vector<double> addCoo;
564 for(mcIdType i=0;i<ncells;i++)
566 if(conn[connI[i]]==ToIdType(INTERP_KERNEL::NORM_SEG2))
568 if(cut3DCurve[i]==-2)
570 mcIdType st=conn[connI[i]+1],endd=conn[connI[i]+2];
571 vec3[0]=coo[3*endd]-coo[3*st]; vec3[1]=coo[3*endd+1]-coo[3*st+1]; vec3[2]=coo[3*endd+2]-coo[3*st+2];
572 double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
573 double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
574 if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
576 const double *st2=coo+3*st;
577 vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
578 double pos=-(vec4[0]*vec2[0]+vec4[1]*vec2[1]+vec4[2]*vec2[2])/((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2]));
579 if(pos>eps && pos<1-eps)
581 mcIdType nNode=ToIdType(addCoo.size())/3;
582 vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
583 addCoo.insert(addCoo.end(),vec4,vec4+3);
584 cut3DCurve[i]=nnodes+nNode;
590 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
594 mcIdType newNbOfNodes=nnodes+ToIdType(addCoo.size())/3;
595 MCAuto<DataArrayDouble> coo2=DataArrayDouble::New();
596 coo2->alloc(newNbOfNodes,3);
597 double *tmp=coo2->getPointer();
598 tmp=std::copy(_coords->begin(),_coords->end(),tmp);
599 std::copy(addCoo.begin(),addCoo.end(),tmp);
600 DataArrayDouble::SetArrayIn(coo2,_coords);
605 * This method incarnates the policy 0 for MEDCouplingUMesh::buildExtrudedMesh method.
606 * \param mesh1D is the input 1D mesh used for translation computation.
607 * \return newCoords new coords filled by this method.
609 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
611 mcIdType oldNbOfNodes=getNumberOfNodes();
612 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
613 std::size_t spaceDim=getSpaceDimension();
614 DataArrayDouble *ret=DataArrayDouble::New();
615 std::vector<bool> isQuads;
616 mcIdType nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
617 ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
618 double *retPtr=ret->getPointer();
619 const double *coords=getCoords()->getConstPointer();
620 double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
621 std::vector<mcIdType> v;
622 std::vector<double> c;
626 for(mcIdType i=0;i<nbOf1DCells;i++)
629 mesh1D->getNodeIdsOfCell(i,v);
631 mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
632 mesh1D->getCoordinatesOfNode(v[0],c);
633 std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
634 for(mcIdType j=0;j<oldNbOfNodes;j++)
635 work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
639 mesh1D->getCoordinatesOfNode(v[1],c);
640 mesh1D->getCoordinatesOfNode(v[0],c);
641 std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
642 for(mcIdType j=0;j<oldNbOfNodes;j++)
643 work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
646 ret->copyStringInfoFrom(*getCoords());
651 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
652 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
653 * \return newCoords new coords filled by this method.
655 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
657 if(mesh1D->getSpaceDimension()==2)
658 return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
659 if(mesh1D->getSpaceDimension()==3)
660 return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
661 throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
665 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
666 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
667 * \return newCoords new coords filled by this method.
669 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
672 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
673 mcIdType oldNbOfNodes=getNumberOfNodes();
674 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
676 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
677 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
678 mcIdType nbOfLevsInVec=nbOf1DCells+1;
679 ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
680 double *retPtr=ret->getPointer();
681 retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
682 MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
683 MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
684 tmp->setCoords(tmp2);
685 const double *coo1D=mesh1D->getCoords()->getConstPointer();
686 const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
687 const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
688 for(mcIdType i=1;i<nbOfLevsInVec;i++)
690 const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
691 const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
692 const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
693 const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
695 double tmp3[2],radius,alpha,alpha0;
696 const double *p0=i+1<nbOfLevsInVec?begin:third;
697 const double *p1=i+1<nbOfLevsInVec?end:begin;
698 const double *p2=i+1<nbOfLevsInVec?third:end;
699 INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
700 double cosangle=i+1<nbOfLevsInVec?(p0[0]-tmp3[0])*(p1[0]-tmp3[0])+(p0[1]-tmp3[1])*(p1[1]-tmp3[1]):(p2[0]-tmp3[0])*(p1[0]-tmp3[0])+(p2[1]-tmp3[1])*(p1[1]-tmp3[1]);
701 double angle=acos(cosangle/(radius*radius));
702 tmp->rotate(end,0,angle);
703 retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
709 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
710 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
711 * \return newCoords new coords filled by this method.
713 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
716 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
717 mcIdType oldNbOfNodes=getNumberOfNodes();
718 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
720 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
721 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
722 mcIdType nbOfLevsInVec=nbOf1DCells+1;
723 ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
724 double *retPtr=ret->getPointer();
725 retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
726 MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
727 MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
728 tmp->setCoords(tmp2);
729 const double *coo1D=mesh1D->getCoords()->getConstPointer();
730 const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
731 const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
732 for(mcIdType i=1;i<nbOfLevsInVec;i++)
734 const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
735 const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
736 const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
737 const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
739 double tmp3[2],radius,alpha,alpha0;
740 const double *p0=i+1<nbOfLevsInVec?begin:third;
741 const double *p1=i+1<nbOfLevsInVec?end:begin;
742 const double *p2=i+1<nbOfLevsInVec?third:end;
744 (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
745 (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
746 (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
748 double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
751 vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
752 double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
753 double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
755 double c2=cos(asin(s2));
757 {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
758 {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
759 {-vec2[1]*s2, vec2[0]*s2, c2}
761 double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]};
762 double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]};
763 double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]};
764 INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
765 double cosangle=i+1<nbOfLevsInVec?(p0r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p0r[1]-tmp3[1])*(p1r[1]-tmp3[1]):(p2r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p2r[1]-tmp3[1])*(p1r[1]-tmp3[1]);
766 double angle=acos(cosangle/(radius*radius));
767 tmp->rotate(end,vecPlane,angle);
769 retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
775 * This method is private because not easy to use for end user. This method is const contrary to
776 * MEDCouplingUMesh::buildExtrudedMesh method because this->_coords are expected to contain
777 * the coords sorted slice by slice.
778 * \param isQuad specifies presence of quadratic cells.
780 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(mcIdType nbOfNodesOf1Lev, bool isQuad) const
782 mcIdType nbOf1DCells(getNumberOfNodes()/nbOfNodesOf1Lev-1);
783 mcIdType nbOf2DCells=getNumberOfCells();
784 mcIdType nbOf3DCells(nbOf2DCells*nbOf1DCells);
785 MEDCouplingUMesh *ret(MEDCouplingUMesh::New("Extruded",getMeshDimension()+1));
786 const mcIdType *conn(_nodal_connec->begin()),*connI(_nodal_connec_index->begin());
787 MCAuto<DataArrayIdType> newConn(DataArrayIdType::New()),newConnI(DataArrayIdType::New());
788 newConnI->alloc(nbOf3DCells+1,1);
789 mcIdType *newConnIPtr(newConnI->getPointer());
791 std::vector<mcIdType> newc;
792 for(mcIdType j=0;j<nbOf2DCells;j++)
794 AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
795 *newConnIPtr++=ToIdType(newc.size());
797 newConn->alloc(newc.size()*nbOf1DCells,1);
798 mcIdType *newConnPtr(newConn->getPointer());
799 mcIdType deltaPerLev(isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev);
800 newConnIPtr=newConnI->getPointer();
801 for(mcIdType iz=0;iz<nbOf1DCells;iz++)
804 std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind2nd(std::plus<mcIdType>(),newConnIPtr[iz*nbOf2DCells]));
805 const mcIdType *posOfTypeOfCell(newConnIPtr);
806 for(std::vector<mcIdType>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
808 mcIdType icell(ToIdType(iter-newc.begin()));//std::distance unfortunately cannot been called here in C++98
809 if(icell!=*posOfTypeOfCell)
812 *newConnPtr=(*iter)+iz*deltaPerLev;
823 ret->setConnectivity(newConn,newConnI,true);
824 ret->setCoords(getCoords());
830 * This method find in candidate pool defined by 'candidates' the cells equal following the polycy 'compType'.
831 * If any true is returned and the results will be put at the end of 'result' output parameter. If not false is returned
832 * and result remains unchanged.
833 * The semantic of 'compType' is specified in MEDCouplingPointSet::zipConnectivityTraducer method.
834 * If in 'candidates' pool -1 value is considered as an empty value.
835 * WARNING this method returns only ONE set of result !
837 bool MEDCouplingUMesh::AreCellsEqualInPool(const std::vector<mcIdType>& candidates, int compType, const mcIdType *conn, const mcIdType *connI, DataArrayIdType *result)
839 if(candidates.size()<1)
842 std::vector<mcIdType>::const_iterator iter=candidates.begin();
843 mcIdType start=(*iter++);
844 for(;iter!=candidates.end();iter++)
846 int status=AreCellsEqual(conn,connI,start,*iter,compType);
851 result->pushBackSilent(start);
855 result->pushBackSilent(*iter);
857 result->pushBackSilent(status==2?(*iter+1):-(*iter+1));
864 * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelf.
865 * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
866 * The return newly allocated mesh will share the same coordinates as \a this.
868 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const mcIdType *begin, const mcIdType *end) const
870 checkConnectivityFullyDefined();
871 mcIdType ncell=getNumberOfCells();
872 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
873 ret->_mesh_dim=_mesh_dim;
874 ret->setCoords(_coords);
875 std::size_t nbOfElemsRet=std::distance(begin,end);
876 mcIdType *connIndexRet=(mcIdType *)malloc((nbOfElemsRet+1)*sizeof(mcIdType));
878 const mcIdType *conn=_nodal_connec->getConstPointer();
879 const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
880 mcIdType newNbring=0;
881 for(const mcIdType *work=begin;work!=end;work++,newNbring++)
883 if(*work>=0 && *work<ncell)
884 connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
888 std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
889 throw INTERP_KERNEL::Exception(oss.str());
892 mcIdType *connRet=(mcIdType *)malloc(connIndexRet[nbOfElemsRet]*sizeof(mcIdType));
893 mcIdType *connRetWork=connRet;
894 std::set<INTERP_KERNEL::NormalizedCellType> types;
895 for(const mcIdType *work=begin;work!=end;work++)
897 types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*work]]);
898 connRetWork=std::copy(conn+connIndex[*work],conn+connIndex[*work+1],connRetWork);
900 MCAuto<DataArrayIdType> connRetArr=DataArrayIdType::New();
901 connRetArr->useArray(connRet,true,DeallocType::C_DEALLOC,connIndexRet[nbOfElemsRet],1);
902 MCAuto<DataArrayIdType> connIndexRetArr=DataArrayIdType::New();
903 connIndexRetArr->useArray(connIndexRet,true,DeallocType::C_DEALLOC,nbOfElemsRet+1,1);
904 ret->setConnectivity(connRetArr,connIndexRetArr,false);
906 ret->copyTinyInfoFrom(this);
911 * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelfSlice.
912 * CellIds are given using range specified by a start an end and step.
914 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice(mcIdType start, mcIdType end, mcIdType step) const
917 mcIdType ncell=getNumberOfCells();
918 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
919 ret->_mesh_dim=_mesh_dim;
920 ret->setCoords(_coords);
921 mcIdType newNbOfCells=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : ");
922 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(newNbOfCells+1,1);
923 mcIdType *newConnIPtr=newConnI->getPointer(); *newConnIPtr=0;
925 const mcIdType *conn=_nodal_connec->getConstPointer();
926 const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
927 for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
929 if(work>=0 && work<ncell)
931 newConnIPtr[1]=newConnIPtr[0]+connIndex[work+1]-connIndex[work];
935 std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : On pos #" << i << " input cell id =" << work << " should be in [0," << ncell << ") !";
936 throw INTERP_KERNEL::Exception(oss.str());
939 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(newConnIPtr[0],1);
940 mcIdType *newConnPtr=newConn->getPointer();
941 std::set<INTERP_KERNEL::NormalizedCellType> types;
943 for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
945 types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[work]]);
946 newConnPtr=std::copy(conn+connIndex[work],conn+connIndex[work+1],newConnPtr);
948 ret->setConnectivity(newConn,newConnI,false);
950 ret->copyTinyInfoFrom(this);
955 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
960 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
966 if(cm.getOrientationStatus(nb,conn1,conn2))
975 * Implements \a conversionType 0 for meshes with meshDim = 1, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
976 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
977 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
979 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic1D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
981 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
982 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
983 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
984 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
985 mcIdType nbOfCells=getNumberOfCells();
986 mcIdType nbOfNodes=getNumberOfNodes();
987 const mcIdType *cPtr=_nodal_connec->begin();
988 const mcIdType *icPtr=_nodal_connec_index->begin();
989 mcIdType lastVal=0,offset=nbOfNodes;
990 for(mcIdType i=0;i<nbOfCells;i++,icPtr++)
992 INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
993 if(type==INTERP_KERNEL::NORM_SEG2)
995 types.insert(INTERP_KERNEL::NORM_SEG3);
996 newConn->pushBackSilent(ToIdType(INTERP_KERNEL::NORM_SEG3));
997 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[0]+3);
998 newConn->pushBackSilent(offset++);
1000 newConnI->pushBackSilent(lastVal);
1001 ret->pushBackSilent(i);
1006 lastVal+=(icPtr[1]-icPtr[0]);
1007 newConnI->pushBackSilent(lastVal);
1008 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1011 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1012 coords=DataArrayDouble::Aggregate(getCoords(),tmp); conn=newConn.retn(); connI=newConnI.retn();
1016 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2DAnd3D0(const MEDCouplingUMesh *m1D, const DataArrayIdType *desc, const DataArrayIdType *descI, DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1018 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1019 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1020 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1022 const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1023 DataArrayIdType *conn1D=0,*conn1DI=0;
1024 std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1025 DataArrayDouble *coordsTmp=0;
1026 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1027 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1028 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1029 const mcIdType *c1DPtr=conn1D->begin();
1030 const mcIdType *c1DIPtr=conn1DI->begin();
1031 mcIdType nbOfCells=getNumberOfCells();
1032 const mcIdType *cPtr=_nodal_connec->begin();
1033 const mcIdType *icPtr=_nodal_connec_index->begin();
1035 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1037 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1038 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1039 if(!cm.isQuadratic())
1041 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType();
1042 types.insert(typ2); newConn->pushBackSilent(typ2);
1043 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1044 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1045 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1046 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0]);
1047 newConnI->pushBackSilent(lastVal);
1048 ret->pushBackSilent(i);
1053 lastVal+=(icPtr[1]-icPtr[0]);
1054 newConnI->pushBackSilent(lastVal);
1055 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1058 conn=newConn.retn(); connI=newConnI.retn(); coords=coordsTmpSafe.retn();
1063 * Implements \a conversionType 0 for meshes with meshDim = 2, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1064 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1065 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1067 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1069 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1070 MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1071 return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1074 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1076 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1077 MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1079 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1080 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1081 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1083 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1084 const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1085 DataArrayIdType *conn1D=0,*conn1DI=0;
1086 std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1087 DataArrayDouble *coordsTmp=0;
1088 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1089 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1090 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1091 const mcIdType *c1DPtr=conn1D->begin();
1092 const mcIdType *c1DIPtr=conn1DI->begin();
1093 mcIdType nbOfCells=getNumberOfCells();
1094 const mcIdType *cPtr=_nodal_connec->begin();
1095 const mcIdType *icPtr=_nodal_connec_index->begin();
1097 mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1098 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1100 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1101 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1102 if(!cm.isQuadratic())
1104 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1105 types.insert(typ2); newConn->pushBackSilent(typ2);
1106 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1107 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1108 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1109 newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1110 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+1;
1111 newConnI->pushBackSilent(lastVal);
1112 ret->pushBackSilent(i);
1117 lastVal+=(icPtr[1]-icPtr[0]);
1118 newConnI->pushBackSilent(lastVal);
1119 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1122 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1123 coords=DataArrayDouble::Aggregate(coordsTmpSafe,tmp); conn=newConn.retn(); connI=newConnI.retn();
1128 * Implements \a conversionType 0 for meshes with meshDim = 3, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1129 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1130 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1132 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1134 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1135 MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1136 return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1139 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1141 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),desc2I(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1142 MCAuto<MEDCouplingUMesh> m2D=buildDescendingConnectivityGen<MinusOneSonsGeneratorBiQuadratic>(desc2,desc2I,tmp2,tmp3,MEDCouplingFastNbrer); tmp2=0; tmp3=0;
1143 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),desc1I(DataArrayIdType::New()),tmp4(DataArrayIdType::New()),tmp5(DataArrayIdType::New());
1144 MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc1,desc1I,tmp4,tmp5); tmp4=0; tmp5=0;
1146 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1147 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1148 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(),ret2=DataArrayIdType::New(); ret->alloc(0,1); ret2->alloc(0,1);
1150 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1151 const mcIdType *descPtr(desc1->begin()),*descIPtr(desc1I->begin()),*desc2Ptr(desc2->begin()),*desc2IPtr(desc2I->begin());
1152 DataArrayIdType *conn1D=0,*conn1DI=0,*conn2D=0,*conn2DI=0;
1153 std::set<INTERP_KERNEL::NormalizedCellType> types1D,types2D;
1154 DataArrayDouble *coordsTmp=0,*coordsTmp2=0;
1155 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=DataArrayIdType::New(); ret1D->alloc(0,1);
1156 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1157 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1158 MCAuto<DataArrayIdType> ret2D=m2D->convertLinearCellsToQuadratic2D1(conn2D,conn2DI,coordsTmp2,types2D); ret2D=DataArrayIdType::New(); ret2D->alloc(0,1);
1159 MCAuto<DataArrayDouble> coordsTmp2Safe(coordsTmp2);
1160 MCAuto<DataArrayIdType> conn2DSafe(conn2D),conn2DISafe(conn2DI);
1161 const mcIdType *c1DPtr=conn1D->begin(),*c1DIPtr=conn1DI->begin(),*c2DPtr=conn2D->begin(),*c2DIPtr=conn2DI->begin();
1162 mcIdType nbOfCells=getNumberOfCells();
1163 const mcIdType *cPtr=_nodal_connec->begin();
1164 const mcIdType *icPtr=_nodal_connec_index->begin();
1166 mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1167 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++,desc2IPtr++)
1169 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1170 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1171 if(!cm.isQuadratic())
1173 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1174 if(typ2==INTERP_KERNEL::NORM_ERROR)
1176 std::ostringstream oss; oss << "MEDCouplingUMesh::convertLinearCellsToQuadratic3D1 : On cell #" << i << " the linear cell type does not support advanced quadratization !";
1177 throw INTERP_KERNEL::Exception(oss.str());
1179 types.insert(typ2); newConn->pushBackSilent(typ2);
1180 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1181 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1182 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1183 for(const mcIdType *d=desc2Ptr+desc2IPtr[0];d!=desc2Ptr+desc2IPtr[1];d++)
1185 mcIdType nodeId2=c2DPtr[c2DIPtr[(*d)+1]-1];
1186 mcIdType tmpPos=newConn->getNumberOfTuples();
1187 newConn->pushBackSilent(nodeId2);
1188 ret2D->pushBackSilent(nodeId2); ret1D->pushBackSilent(tmpPos);
1190 newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1191 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+(desc2IPtr[1]-desc2IPtr[0])+1;
1192 newConnI->pushBackSilent(lastVal);
1193 ret->pushBackSilent(i);
1198 lastVal+=(icPtr[1]-icPtr[0]);
1199 newConnI->pushBackSilent(lastVal);
1200 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1203 MCAuto<DataArrayIdType> diffRet2D=ret2D->getDifferentValues();
1204 MCAuto<DataArrayIdType> o2nRet2D=diffRet2D->invertArrayN2O2O2N(coordsTmp2Safe->getNumberOfTuples());
1205 coordsTmp2Safe=coordsTmp2Safe->selectByTupleId(diffRet2D->begin(),diffRet2D->end());
1206 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1207 std::vector<const DataArrayDouble *> v(3); v[0]=coordsTmpSafe; v[1]=coordsTmp2Safe; v[2]=tmp;
1208 mcIdType *c=newConn->getPointer();
1209 const mcIdType *cI(newConnI->begin());
1210 for(const mcIdType *elt=ret1D->begin();elt!=ret1D->end();elt++)
1211 c[*elt]=o2nRet2D->getIJ(c[*elt],0)+offset;
1212 offset=coordsTmp2Safe->getNumberOfTuples();
1213 for(const mcIdType *elt=ret->begin();elt!=ret->end();elt++)
1214 c[cI[(*elt)+1]-1]+=offset;
1215 coords=DataArrayDouble::Aggregate(v); conn=newConn.retn(); connI=newConnI.retn();
1219 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1221 mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1222 const mcIdType *n2oPtr(n2o->begin());
1223 MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1224 skin->getReverseNodalConnectivity(revNodal,revNodalI);
1225 const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1226 const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1227 const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1228 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1229 mcIdType *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_POLYGON;
1230 if(nbOfNodesExpected<1)
1232 mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1233 *work++=n2oPtr[prevNode];
1234 for(mcIdType i=1;i<nbOfNodesExpected;i++)
1236 if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==3)
1238 std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1239 conn.erase(prevNode);
1242 mcIdType curNode(*(conn.begin()));
1243 *work++=n2oPtr[curNode];
1244 std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1245 shar.erase(prevCell);
1248 prevCell=*(shar.begin());
1252 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
1255 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
1258 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
1263 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1265 mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1266 mcIdType nbOfTurn(nbOfNodesExpected/2);
1267 const mcIdType *n2oPtr(n2o->begin());
1268 MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1269 skin->getReverseNodalConnectivity(revNodal,revNodalI);
1270 const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1271 const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1272 const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1273 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1274 mcIdType *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_QPOLYG;
1275 if(nbOfNodesExpected<1)
1277 mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1278 *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
1279 for(mcIdType i=1;i<nbOfTurn;i++)
1281 if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
1283 std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1284 conn.erase(prevNode);
1287 mcIdType curNode(*(conn.begin()));
1288 *work=n2oPtr[curNode];
1289 std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1290 shar.erase(prevCell);
1293 mcIdType curCell(*(shar.begin()));
1294 work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
1300 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
1303 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
1306 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
1311 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
1314 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
1315 std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
1316 int meshDim=(*it)->getMeshDimension();
1317 mcIdType nbOfCells=ToIdType((*it)->getNumberOfCells());
1318 mcIdType meshLgth=(*it++)->getNodalConnectivityArrayLen();
1319 for(;it!=a.end();it++)
1321 if(meshDim!=(*it)->getMeshDimension())
1322 throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
1323 nbOfCells+=ToIdType((*it)->getNumberOfCells());
1324 meshLgth+=(*it)->getNodalConnectivityArrayLen();
1326 std::vector<const MEDCouplingPointSet *> aps(a.size());
1327 std::copy(a.begin(),a.end(),aps.begin());
1328 MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1329 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
1330 ret->setCoords(pts);
1331 MCAuto<DataArrayIdType> c=DataArrayIdType::New();
1332 c->alloc(meshLgth,1);
1333 mcIdType *cPtr=c->getPointer();
1334 MCAuto<DataArrayIdType> cI=DataArrayIdType::New();
1335 cI->alloc(nbOfCells+1,1);
1336 mcIdType *cIPtr=cI->getPointer();
1340 for(it=a.begin();it!=a.end();it++)
1342 mcIdType curNbOfCell=ToIdType((*it)->getNumberOfCells());
1343 const mcIdType *curCI=(*it)->_nodal_connec_index->begin();
1344 const mcIdType *curC=(*it)->_nodal_connec->begin();
1345 cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind2nd(std::plus<mcIdType>(),offset));
1346 for(mcIdType j=0;j<curNbOfCell;j++)
1348 const mcIdType *src=curC+curCI[j];
1350 for(;src!=curC+curCI[j+1];src++,cPtr++)
1358 offset+=curCI[curNbOfCell];
1359 offset2+=(*it)->getNumberOfNodes();
1362 ret->setConnectivity(c,cI,true);
1368 * \param [in] pt the start pointer (included) of the coordinates of the point
1369 * \param [in] cellIdsBg the start pointer (included) of cellIds
1370 * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1371 * \param [in] nc nodal connectivity
1372 * \param [in] ncI nodal connectivity index
1373 * \param [in,out] ret0 the min distance between \a this and the external input point
1374 * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1375 * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1377 void MEDCouplingUMesh::DistanceToPoint3DSurfAlg(const double *pt, const mcIdType *cellIdsBg, const mcIdType *cellIdsEnd, const double *coords, const mcIdType *nc, const mcIdType *ncI, double& ret0, mcIdType& cellId)
1380 ret0=std::numeric_limits<double>::max();
1381 for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1383 switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1385 case INTERP_KERNEL::NORM_TRI3:
1387 double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
1389 { ret0=tmp; cellId=*zeCell; }
1392 case INTERP_KERNEL::NORM_QUAD4:
1393 case INTERP_KERNEL::NORM_POLYGON:
1395 double tmp=INTERP_KERNEL::DistanceFromPtToPolygonInSpaceDim3(pt,nc+ncI[*zeCell]+1,nc+ncI[*zeCell+1],coords);
1397 { ret0=tmp; cellId=*zeCell; }
1401 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
1407 * \param [in] pt the start pointer (included) of the coordinates of the point
1408 * \param [in] cellIdsBg the start pointer (included) of cellIds
1409 * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1410 * \param [in] nc nodal connectivity
1411 * \param [in] ncI nodal connectivity index
1412 * \param [in,out] ret0 the min distance between \a this and the external input point
1413 * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1414 * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1416 void MEDCouplingUMesh::DistanceToPoint2DCurveAlg(const double *pt, const mcIdType *cellIdsBg, const mcIdType *cellIdsEnd, const double *coords, const mcIdType *nc, const mcIdType *ncI, double& ret0, mcIdType& cellId)
1419 ret0=std::numeric_limits<double>::max();
1420 for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1422 switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1424 case INTERP_KERNEL::NORM_SEG2:
1426 std::size_t uselessEntry=0;
1427 double tmp=INTERP_KERNEL::SquareDistanceFromPtToSegInSpaceDim2(pt,coords+2*nc[ncI[*zeCell]+1],coords+2*nc[ncI[*zeCell]+2],uselessEntry);
1430 { ret0=tmp; cellId=*zeCell; }
1434 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
1438 DataArrayIdType *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg(std::vector<bool>& fetched, const mcIdType *seedBg, const mcIdType *seedEnd, const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, mcIdType nbOfDepthPeeling, mcIdType& nbOfDepthPeelingPerformed)
1440 nbOfDepthPeelingPerformed=0;
1441 if(!seedBg || !seedEnd || !arrIn || !arrIndxIn)
1442 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : some input pointer is NULL !");
1443 mcIdType nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
1444 std::vector<bool> fetched2(nbOfTuples,false);
1446 for(const mcIdType *seedElt=seedBg;seedElt!=seedEnd;seedElt++,i++)
1448 if(*seedElt>=0 && *seedElt<nbOfTuples)
1449 { fetched[*seedElt]=true; fetched2[*seedElt]=true; }
1451 { std::ostringstream oss; oss << "MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : At pos #" << i << " of seeds value is " << *seedElt << "! Should be in [0," << nbOfTuples << ") !"; throw INTERP_KERNEL::Exception(oss.str()); }
1453 const mcIdType *arrInPtr=arrIn->begin();
1454 const mcIdType *arrIndxPtr=arrIndxIn->begin();
1455 mcIdType targetNbOfDepthPeeling=nbOfDepthPeeling!=-1?nbOfDepthPeeling:std::numeric_limits<mcIdType>::max();
1456 std::vector<mcIdType> idsToFetch1(seedBg,seedEnd);
1457 std::vector<mcIdType> idsToFetch2;
1458 std::vector<mcIdType> *idsToFetch=&idsToFetch1;
1459 std::vector<mcIdType> *idsToFetchOther=&idsToFetch2;
1460 while(!idsToFetch->empty() && nbOfDepthPeelingPerformed<targetNbOfDepthPeeling)
1462 for(std::vector<mcIdType>::const_iterator it=idsToFetch->begin();it!=idsToFetch->end();it++)
1463 for(const mcIdType *it2=arrInPtr+arrIndxPtr[*it];it2!=arrInPtr+arrIndxPtr[*it+1];it2++)
1465 { fetched[*it2]=true; fetched2[*it2]=true; idsToFetchOther->push_back(*it2); }
1466 std::swap(idsToFetch,idsToFetchOther);
1467 idsToFetchOther->clear();
1468 nbOfDepthPeelingPerformed++;
1470 mcIdType lgth=ToIdType(std::count(fetched2.begin(),fetched2.end(),true));
1472 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(lgth,1);
1473 mcIdType *retPtr=ret->getPointer();
1474 for(std::vector<bool>::const_iterator it=fetched2.begin();it!=fetched2.end();it++,i++)
1481 * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
1482 * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
1484 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, mcIdType nbOfNodesInCell, const mcIdType *conn, const double *coo, double *zipFrmt)
1488 for(int i=0;i<nbOfNodesInCell;i++)
1489 w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
1490 else if(spaceDim==2)
1492 for(int i=0;i<nbOfNodesInCell;i++)
1494 w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
1499 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
1503 * This method takes in input a cell defined by its MEDcouplingUMesh connectivity [ \a connBg , \a connEnd ) and returns its extruded cell by inserting the result at the end of ret.
1504 * \param nbOfNodesPerLev in parameter that specifies the number of nodes of one slice of global dataset
1505 * \param isQuad specifies the policy of connectivity.
1506 * @ret in/out parameter in which the result will be append
1508 void MEDCouplingUMesh::AppendExtrudedCell(const mcIdType *connBg, const mcIdType *connEnd, mcIdType nbOfNodesPerLev, bool isQuad, std::vector<mcIdType>& ret)
1510 INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
1511 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
1512 ret.push_back(cm.getExtrudedType());
1513 mcIdType deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
1516 case INTERP_KERNEL::NORM_POINT1:
1518 ret.push_back(connBg[1]);
1519 ret.push_back(connBg[1]+nbOfNodesPerLev);
1522 case INTERP_KERNEL::NORM_SEG2:
1524 mcIdType conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
1525 ret.insert(ret.end(),conn,conn+4);
1528 case INTERP_KERNEL::NORM_SEG3:
1530 mcIdType conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
1531 ret.insert(ret.end(),conn,conn+8);
1534 case INTERP_KERNEL::NORM_QUAD4:
1536 mcIdType conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
1537 ret.insert(ret.end(),conn,conn+8);
1540 case INTERP_KERNEL::NORM_TRI3:
1542 mcIdType conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
1543 ret.insert(ret.end(),conn,conn+6);
1546 case INTERP_KERNEL::NORM_TRI6:
1548 mcIdType conn[15]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4],connBg[5],connBg[6],connBg[4]+deltaz,connBg[5]+deltaz,connBg[6]+deltaz,
1549 connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
1550 ret.insert(ret.end(),conn,conn+15);
1553 case INTERP_KERNEL::NORM_QUAD8:
1556 connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
1557 connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
1558 connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
1560 ret.insert(ret.end(),conn,conn+20);
1563 case INTERP_KERNEL::NORM_POLYGON:
1565 std::back_insert_iterator< std::vector<mcIdType> > ii(ret);
1566 std::copy(connBg+1,connEnd,ii);
1568 std::reverse_iterator<const mcIdType *> rConnBg(connEnd);
1569 std::reverse_iterator<const mcIdType *> rConnEnd(connBg+1);
1570 std::transform(rConnBg,rConnEnd,ii,std::bind2nd(std::plus<mcIdType>(),deltaz));
1571 std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
1572 for(std::size_t i=0;i<nbOfRadFaces;i++)
1575 mcIdType conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
1576 std::copy(conn,conn+4,ii);
1581 throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
1587 * This method is part of the Slice3D algorithm. It is the first step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
1588 * This method allows to compute given the status of 3D curve cells and the descending connectivity 3DSurf->3DCurve to deduce the intersection of each 3D surf cells
1589 * with a plane. The result will be put in 'cut3DSuf' out parameter.
1590 * \param [in] cut3DCurve input parameter that gives for each 3DCurve cell if it owns fully to the plane or partially.
1591 * \param [out] nodesOnPlane, returns all the nodes that are on the plane.
1592 * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
1593 * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
1594 * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
1595 * \param [in] nodal3DIndxCurve is the nodal connectivity index of 3D curve mesh.
1596 * \param [in] desc is the descending connectivity 3DSurf->3DCurve
1597 * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
1598 * \param [out] cut3DSuf input/output param.
1600 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
1601 const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
1602 const mcIdType *desc, const mcIdType *descIndx,
1603 std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf)
1605 std::set<mcIdType> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
1606 mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
1607 for(mcIdType i=0;i<nbOf3DSurfCell;i++)
1609 std::vector<mcIdType> res;
1610 mcIdType offset=descIndx[i];
1611 mcIdType nbOfSeg=descIndx[i+1]-offset;
1612 for(mcIdType j=0;j<nbOfSeg;j++)
1614 mcIdType edgeId=desc[offset+j];
1615 mcIdType status=cut3DCurve[edgeId];
1619 res.push_back(status);
1622 res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
1623 res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
1631 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1637 std::set<mcIdType> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
1638 std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<mcIdType> >(res));
1641 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1645 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
1650 {// case when plane is on a multi colinear edge of a polyhedron
1651 if(ToIdType(res.size())==2*nbOfSeg)
1653 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
1656 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
1664 * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
1665 * This method is part of the Slice3D algorithm. It is the second step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
1666 * This method allows to compute given the result of 3D surf cells with plane and the descending connectivity 3D->3DSurf to deduce the intersection of each 3D cells
1667 * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
1668 * \param cut3DSurf input parameter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
1669 * \param desc is the descending connectivity 3D->3DSurf
1670 * \param descIndx is the descending connectivity index 3D->3DSurf
1672 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1673 const mcIdType *desc, const mcIdType *descIndx,
1674 DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const
1676 checkFullyDefined();
1677 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1678 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1679 const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1680 mcIdType nbOfCells=getNumberOfCells();
1681 for(mcIdType i=0;i<nbOfCells;i++)
1683 std::map<mcIdType, std::set<mcIdType> > m;
1684 mcIdType offset=descIndx[i];
1685 mcIdType nbOfFaces=descIndx[i+1]-offset;
1688 for(int j=0;j<nbOfFaces;j++)
1690 const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1691 if(p.first!=-1 && p.second!=-1)
1695 start=p.first; end=p.second;
1696 m[p.first].insert(p.second);
1697 m[p.second].insert(p.first);
1701 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
1702 mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1703 INTERP_KERNEL::AutoPtr<mcIdType> tmp=new mcIdType[sz];
1704 INTERP_KERNEL::NormalizedCellType cmsId;
1705 unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
1706 start=tmp[0]; end=tmp[nbOfNodesSon-1];
1707 for(unsigned k=0;k<nbOfNodesSon;k++)
1709 m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
1710 m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
1717 std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
1721 std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
1722 const std::set<mcIdType>& s=(*it).second;
1723 std::set<mcIdType> s2; s2.insert(prev);
1724 std::set<mcIdType> s3;
1725 std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
1728 mcIdType val=*s3.begin();
1729 conn.push_back(start);
1736 conn.push_back(end);
1739 nodalRes->insertAtTheEnd(conn.begin(),conn.end());
1740 nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
1741 cellIds->pushBackSilent(i);
1747 void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayIdType *nodalConnec, const DataArrayIdType *nodalConnecIndex)
1749 if(nodalConnec && nodalConnecIndex)
1752 const mcIdType *conn(nodalConnec->begin()),*connIndex(nodalConnecIndex->begin());
1753 mcIdType nbOfElem=ToIdType(nodalConnecIndex->getNbOfElems())-1;
1755 for(const mcIdType *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
1756 types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
1761 * This method expects that \a this a quadratic 1D, 2D or 3D mesh.
1762 * This method will 'attract' middle points of seg3 (deduced from this by explosion if needed) of cells connected to nodes specified in [\a nodeIdsBg, \a nodeIdsEnd )
1763 * For those selected mid points, their coordinates will be modified by applying a dilation between node in input [\a nodeIdsBg, \a nodeIdsEnd ) and the corresponding mid points using \a ratio input value.
1764 * So this method is non const because coordinates are modified.
1765 * If there is a couple of 2 points in [\a nodeIdsBg, \a nodeIdsEnd ) that are boundaries of a seg3, the corresponding mid point will remain untouched.
1767 * \param [in] ratio - ratio of dilation
1768 * \param [in] nodeIdsBg - start (included) of input node list
1769 * \param [in] nodeIdsEnd - end (excluded) of input node list
1770 * \throw if there is a point in [\a nodeIdsBg, \a nodeIdsEnd ) that is a mid point of a seg3
1771 * \warning in case of throw the coordinates may be partially modified before the exception arises
1773 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodes(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1775 checkFullyDefined();
1776 int mdim(getMeshDimension());
1777 if(mdim==2 || mdim==3)
1779 MCAuto<MEDCouplingUMesh> edges;
1781 MCAuto<DataArrayIdType> a,b,c,d;
1782 edges=this->explodeIntoEdges(a,b,c,d);
1784 return edges->attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1787 return attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1788 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : not managed dimension ! Should be in [1,2,3] !");
1792 * \a this is expected to have meshdim==1.
1794 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodesUnderground(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1796 int spaceDim(getSpaceDimension());
1797 double *coords(getCoords()->getPointer());
1798 auto nbNodes(getNumberOfNodes());
1799 std::size_t nbCells(getNumberOfCells());
1800 std::vector<bool> fastFinder(nbNodes,false);
1801 for(auto work=nodeIdsBg;work!=nodeIdsEnd;work++)
1802 if(*work>=0 && *work<nbNodes)
1803 fastFinder[*work]=true;
1804 MCAuto<DataArrayIdType> cellsIds(getCellIdsLyingOnNodes(nodeIdsBg,nodeIdsEnd,false));
1805 const mcIdType *nc(_nodal_connec->begin()),*nci(_nodal_connec_index->begin());
1806 for(std::size_t cellId=0;cellId<nbCells;cellId++,nci++)
1808 const mcIdType *isSelected(std::find_if(nc+nci[0]+1,nc+nci[1],[&fastFinder](mcIdType v) { return fastFinder[v]; }));
1809 if(isSelected!=nc+nci[1])
1811 if((INTERP_KERNEL::NormalizedCellType)nc[nci[0]]==INTERP_KERNEL::NORM_SEG3 && nci[1]-nci[0]==4)
1813 bool aa(fastFinder[nc[*nci+1]]),bb(fastFinder[nc[*nci+2]]),cc(fastFinder[nc[*nci+3]]);
1818 auto ptToMove(nc[*nci+3]);
1819 auto attractor(aa?nc[*nci+1]:nc[*nci+2]),endPt(aa?nc[*nci+2]:nc[*nci+1]);
1820 std::transform(coords+spaceDim*attractor,coords+spaceDim*(attractor+1),coords+spaceDim*endPt,
1821 coords+spaceDim*ptToMove,[ratio](const double& stPt, const double& endPt) { return stPt+ratio*(endPt-stPt); });
1824 continue;//both 2 boundary nodes of current seg3 are un nodeIds input list -> skip it.
1828 std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " has a mid point " << nc[*nci+3] << " ! This node is in input list !";
1829 throw INTERP_KERNEL::Exception(oss.str());
1834 std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " sharing one of the input nodes list its geo type is NOT SEG3 !";
1835 throw INTERP_KERNEL::Exception(oss.str());