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 if(getNodalConnectivity()->empty())
515 cellIdsKeptArr=cellIdsKept.retn();
518 mcIdType sz=getNodalConnectivity()->getMaxValue(tmp); sz=std::max(sz,ToIdType(0))+1;
519 std::vector<bool> fastFinder(sz,false);
520 for(const mcIdType *work=begin;work!=end;work++)
521 if(*work>=0 && *work<sz)
522 fastFinder[*work]=true;
523 mcIdType nbOfCells=getNumberOfCells();
524 const mcIdType *conn=getNodalConnectivity()->getConstPointer();
525 const mcIdType *connIndex=getNodalConnectivityIndex()->getConstPointer();
526 for(mcIdType i=0;i<nbOfCells;i++)
528 mcIdType ref=0,nbOfHit=0;
529 for(const mcIdType *work2=conn+connIndex[i]+1;work2!=conn+connIndex[i+1];work2++)
533 if(fastFinder[*work2])
536 if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
537 cellIdsKept->pushBackSilent(i);
539 cellIdsKeptArr=cellIdsKept.retn();
543 * This method works on a 3D curve linear mesh that is to say (meshDim==1 and spaceDim==3).
544 * If it is not the case an exception will be thrown.
545 * This method is non const because the coordinate of \a this can be appended with some new points issued from
546 * intersection of plane defined by ('origin','vec').
547 * This method has one in/out parameter : 'cut3DCurve'.
548 * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
549 * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
550 * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
551 * This method will throw an exception if \a this contains a non linear segment.
553 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<mcIdType>& cut3DCurve)
556 if(getMeshDimension()!=1 || getSpaceDimension()!=3)
557 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
558 mcIdType ncells=getNumberOfCells();
559 mcIdType nnodes=getNumberOfNodes();
560 double vec2[3],vec3[3],vec4[3];
561 double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
563 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
564 vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
565 const mcIdType *conn=_nodal_connec->getConstPointer();
566 const mcIdType *connI=_nodal_connec_index->getConstPointer();
567 const double *coo=_coords->getConstPointer();
568 std::vector<double> addCoo;
569 for(mcIdType i=0;i<ncells;i++)
571 if(conn[connI[i]]==ToIdType(INTERP_KERNEL::NORM_SEG2))
573 if(cut3DCurve[i]==-2)
575 mcIdType st=conn[connI[i]+1],endd=conn[connI[i]+2];
576 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];
577 double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
578 double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
579 if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
581 const double *st2=coo+3*st;
582 vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
583 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]));
584 if(pos>eps && pos<1-eps)
586 mcIdType nNode=ToIdType(addCoo.size())/3;
587 vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
588 addCoo.insert(addCoo.end(),vec4,vec4+3);
589 cut3DCurve[i]=nnodes+nNode;
595 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
599 mcIdType newNbOfNodes=nnodes+ToIdType(addCoo.size())/3;
600 MCAuto<DataArrayDouble> coo2=DataArrayDouble::New();
601 coo2->alloc(newNbOfNodes,3);
602 double *tmp=coo2->getPointer();
603 tmp=std::copy(_coords->begin(),_coords->end(),tmp);
604 std::copy(addCoo.begin(),addCoo.end(),tmp);
605 DataArrayDouble::SetArrayIn(coo2,_coords);
610 * This method incarnates the policy 0 for MEDCouplingUMesh::buildExtrudedMesh method.
611 * \param mesh1D is the input 1D mesh used for translation computation.
612 * \return newCoords new coords filled by this method.
614 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
616 mcIdType oldNbOfNodes=getNumberOfNodes();
617 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
618 std::size_t spaceDim=getSpaceDimension();
619 DataArrayDouble *ret=DataArrayDouble::New();
620 std::vector<bool> isQuads;
621 mcIdType nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
622 ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
623 double *retPtr=ret->getPointer();
624 const double *coords=getCoords()->getConstPointer();
625 double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
626 std::vector<mcIdType> v;
627 std::vector<double> c;
631 for(mcIdType i=0;i<nbOf1DCells;i++)
634 mesh1D->getNodeIdsOfCell(i,v);
636 mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
637 mesh1D->getCoordinatesOfNode(v[0],c);
638 std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
639 for(mcIdType j=0;j<oldNbOfNodes;j++)
640 work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
644 mesh1D->getCoordinatesOfNode(v[1],c);
645 mesh1D->getCoordinatesOfNode(v[0],c);
646 std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
647 for(mcIdType j=0;j<oldNbOfNodes;j++)
648 work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
651 ret->copyStringInfoFrom(*getCoords());
656 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
657 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
658 * \return newCoords new coords filled by this method.
660 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
662 if(mesh1D->getSpaceDimension()==2)
663 return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
664 if(mesh1D->getSpaceDimension()==3)
665 return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
666 throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
670 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
671 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
672 * \return newCoords new coords filled by this method.
674 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
677 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
678 mcIdType oldNbOfNodes=getNumberOfNodes();
679 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
681 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
682 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
683 mcIdType nbOfLevsInVec=nbOf1DCells+1;
684 ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
685 double *retPtr=ret->getPointer();
686 retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
687 MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
688 MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
689 tmp->setCoords(tmp2);
690 const double *coo1D=mesh1D->getCoords()->getConstPointer();
691 const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
692 const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
693 for(mcIdType i=1;i<nbOfLevsInVec;i++)
695 const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
696 const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
697 const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
698 const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
700 double tmp3[2],radius,alpha,alpha0;
701 const double *p0=i+1<nbOfLevsInVec?begin:third;
702 const double *p1=i+1<nbOfLevsInVec?end:begin;
703 const double *p2=i+1<nbOfLevsInVec?third:end;
704 INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
705 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]);
706 double angle=acos(cosangle/(radius*radius));
707 tmp->rotate(end,0,angle);
708 retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
714 * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
715 * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
716 * \return newCoords new coords filled by this method.
718 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
721 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
722 mcIdType oldNbOfNodes=getNumberOfNodes();
723 mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
725 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
726 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
727 mcIdType nbOfLevsInVec=nbOf1DCells+1;
728 ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
729 double *retPtr=ret->getPointer();
730 retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
731 MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
732 MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
733 tmp->setCoords(tmp2);
734 const double *coo1D=mesh1D->getCoords()->getConstPointer();
735 const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
736 const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
737 for(mcIdType i=1;i<nbOfLevsInVec;i++)
739 const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
740 const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
741 const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
742 const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
744 double tmp3[2],radius,alpha,alpha0;
745 const double *p0=i+1<nbOfLevsInVec?begin:third;
746 const double *p1=i+1<nbOfLevsInVec?end:begin;
747 const double *p2=i+1<nbOfLevsInVec?third:end;
749 (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
750 (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
751 (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
753 double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
756 vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
757 double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
758 double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
760 double c2=cos(asin(s2));
762 {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
763 {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
764 {-vec2[1]*s2, vec2[0]*s2, c2}
766 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]};
767 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]};
768 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]};
769 INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
770 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]);
771 double angle=acos(cosangle/(radius*radius));
772 tmp->rotate(end,vecPlane,angle);
774 retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
780 * This method is private because not easy to use for end user. This method is const contrary to
781 * MEDCouplingUMesh::buildExtrudedMesh method because this->_coords are expected to contain
782 * the coords sorted slice by slice.
783 * \param isQuad specifies presence of quadratic cells.
785 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(mcIdType nbOfNodesOf1Lev, bool isQuad) const
787 mcIdType nbOf1DCells(getNumberOfNodes()/nbOfNodesOf1Lev-1);
788 mcIdType nbOf2DCells=getNumberOfCells();
789 mcIdType nbOf3DCells(nbOf2DCells*nbOf1DCells);
790 MEDCouplingUMesh *ret(MEDCouplingUMesh::New("Extruded",getMeshDimension()+1));
791 const mcIdType *conn(_nodal_connec->begin()),*connI(_nodal_connec_index->begin());
792 MCAuto<DataArrayIdType> newConn(DataArrayIdType::New()),newConnI(DataArrayIdType::New());
793 newConnI->alloc(nbOf3DCells+1,1);
794 mcIdType *newConnIPtr(newConnI->getPointer());
796 std::vector<mcIdType> newc;
797 for(mcIdType j=0;j<nbOf2DCells;j++)
799 AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
800 *newConnIPtr++=ToIdType(newc.size());
802 newConn->alloc(newc.size()*nbOf1DCells,1);
803 mcIdType *newConnPtr(newConn->getPointer());
804 mcIdType deltaPerLev(isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev);
805 newConnIPtr=newConnI->getPointer();
806 for(mcIdType iz=0;iz<nbOf1DCells;iz++)
809 std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind2nd(std::plus<mcIdType>(),newConnIPtr[iz*nbOf2DCells]));
810 const mcIdType *posOfTypeOfCell(newConnIPtr);
811 for(std::vector<mcIdType>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
813 mcIdType icell(ToIdType(iter-newc.begin()));//std::distance unfortunately cannot been called here in C++98
814 if(icell!=*posOfTypeOfCell)
817 *newConnPtr=(*iter)+iz*deltaPerLev;
828 ret->setConnectivity(newConn,newConnI,true);
829 ret->setCoords(getCoords());
835 * This method find in candidate pool defined by 'candidates' the cells equal following the polycy 'compType'.
836 * If any true is returned and the results will be put at the end of 'result' output parameter. If not false is returned
837 * and result remains unchanged.
838 * The semantic of 'compType' is specified in MEDCouplingPointSet::zipConnectivityTraducer method.
839 * If in 'candidates' pool -1 value is considered as an empty value.
840 * WARNING this method returns only ONE set of result !
842 bool MEDCouplingUMesh::AreCellsEqualInPool(const std::vector<mcIdType>& candidates, int compType, const mcIdType *conn, const mcIdType *connI, DataArrayIdType *result)
844 if(candidates.size()<1)
847 std::vector<mcIdType>::const_iterator iter=candidates.begin();
848 mcIdType start=(*iter++);
849 for(;iter!=candidates.end();iter++)
851 int status=AreCellsEqual(conn,connI,start,*iter,compType);
856 result->pushBackSilent(start);
860 result->pushBackSilent(*iter);
862 result->pushBackSilent(status==2?(*iter+1):-(*iter+1));
869 * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelf.
870 * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
871 * The return newly allocated mesh will share the same coordinates as \a this.
873 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const mcIdType *begin, const mcIdType *end) const
875 checkConnectivityFullyDefined();
876 mcIdType ncell=getNumberOfCells();
877 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
878 ret->_mesh_dim=_mesh_dim;
879 ret->setCoords(_coords);
880 std::size_t nbOfElemsRet=std::distance(begin,end);
881 mcIdType *connIndexRet=(mcIdType *)malloc((nbOfElemsRet+1)*sizeof(mcIdType));
883 const mcIdType *conn=_nodal_connec->getConstPointer();
884 const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
885 mcIdType newNbring=0;
886 for(const mcIdType *work=begin;work!=end;work++,newNbring++)
888 if(*work>=0 && *work<ncell)
889 connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
893 std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
894 throw INTERP_KERNEL::Exception(oss.str());
897 mcIdType *connRet=(mcIdType *)malloc(connIndexRet[nbOfElemsRet]*sizeof(mcIdType));
898 mcIdType *connRetWork=connRet;
899 std::set<INTERP_KERNEL::NormalizedCellType> types;
900 for(const mcIdType *work=begin;work!=end;work++)
902 types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*work]]);
903 connRetWork=std::copy(conn+connIndex[*work],conn+connIndex[*work+1],connRetWork);
905 MCAuto<DataArrayIdType> connRetArr=DataArrayIdType::New();
906 connRetArr->useArray(connRet,true,DeallocType::C_DEALLOC,connIndexRet[nbOfElemsRet],1);
907 MCAuto<DataArrayIdType> connIndexRetArr=DataArrayIdType::New();
908 connIndexRetArr->useArray(connIndexRet,true,DeallocType::C_DEALLOC,nbOfElemsRet+1,1);
909 ret->setConnectivity(connRetArr,connIndexRetArr,false);
911 ret->copyTinyInfoFrom(this);
916 * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelfSlice.
917 * CellIds are given using range specified by a start an end and step.
919 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice(mcIdType start, mcIdType end, mcIdType step) const
922 mcIdType ncell=getNumberOfCells();
923 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
924 ret->_mesh_dim=_mesh_dim;
925 ret->setCoords(_coords);
926 mcIdType newNbOfCells=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : ");
927 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(newNbOfCells+1,1);
928 mcIdType *newConnIPtr=newConnI->getPointer(); *newConnIPtr=0;
930 const mcIdType *conn=_nodal_connec->getConstPointer();
931 const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
932 for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
934 if(work>=0 && work<ncell)
936 newConnIPtr[1]=newConnIPtr[0]+connIndex[work+1]-connIndex[work];
940 std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : On pos #" << i << " input cell id =" << work << " should be in [0," << ncell << ") !";
941 throw INTERP_KERNEL::Exception(oss.str());
944 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(newConnIPtr[0],1);
945 mcIdType *newConnPtr=newConn->getPointer();
946 std::set<INTERP_KERNEL::NormalizedCellType> types;
948 for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
950 types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[work]]);
951 newConnPtr=std::copy(conn+connIndex[work],conn+connIndex[work+1],newConnPtr);
953 ret->setConnectivity(newConn,newConnI,false);
955 ret->copyTinyInfoFrom(this);
960 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
965 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
971 if(cm.getOrientationStatus(nb,conn1,conn2))
980 * Implements \a conversionType 0 for meshes with meshDim = 1, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
981 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
982 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
984 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic1D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
986 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
987 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
988 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
989 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
990 mcIdType nbOfCells=getNumberOfCells();
991 mcIdType nbOfNodes=getNumberOfNodes();
992 const mcIdType *cPtr=_nodal_connec->begin();
993 const mcIdType *icPtr=_nodal_connec_index->begin();
994 mcIdType lastVal=0,offset=nbOfNodes;
995 for(mcIdType i=0;i<nbOfCells;i++,icPtr++)
997 INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
998 if(type==INTERP_KERNEL::NORM_SEG2)
1000 types.insert(INTERP_KERNEL::NORM_SEG3);
1001 newConn->pushBackSilent(ToIdType(INTERP_KERNEL::NORM_SEG3));
1002 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[0]+3);
1003 newConn->pushBackSilent(offset++);
1005 newConnI->pushBackSilent(lastVal);
1006 ret->pushBackSilent(i);
1011 lastVal+=(icPtr[1]-icPtr[0]);
1012 newConnI->pushBackSilent(lastVal);
1013 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1016 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1017 coords=DataArrayDouble::Aggregate(getCoords(),tmp); conn=newConn.retn(); connI=newConnI.retn();
1021 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
1023 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1024 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1025 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1027 const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1028 DataArrayIdType *conn1D=0,*conn1DI=0;
1029 std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1030 DataArrayDouble *coordsTmp=0;
1031 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1032 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1033 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1034 const mcIdType *c1DPtr=conn1D->begin();
1035 const mcIdType *c1DIPtr=conn1DI->begin();
1036 mcIdType nbOfCells=getNumberOfCells();
1037 const mcIdType *cPtr=_nodal_connec->begin();
1038 const mcIdType *icPtr=_nodal_connec_index->begin();
1040 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1042 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1043 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1044 if(!cm.isQuadratic())
1046 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType();
1047 types.insert(typ2); newConn->pushBackSilent(typ2);
1048 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1049 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1050 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1051 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0]);
1052 newConnI->pushBackSilent(lastVal);
1053 ret->pushBackSilent(i);
1058 lastVal+=(icPtr[1]-icPtr[0]);
1059 newConnI->pushBackSilent(lastVal);
1060 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1063 conn=newConn.retn(); connI=newConnI.retn(); coords=coordsTmpSafe.retn();
1068 * Implements \a conversionType 0 for meshes with meshDim = 2, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1069 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1070 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1072 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1074 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1075 MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1076 return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1079 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1081 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1082 MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1084 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1085 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1086 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1088 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1089 const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1090 DataArrayIdType *conn1D=0,*conn1DI=0;
1091 std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1092 DataArrayDouble *coordsTmp=0;
1093 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1094 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1095 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1096 const mcIdType *c1DPtr=conn1D->begin();
1097 const mcIdType *c1DIPtr=conn1DI->begin();
1098 mcIdType nbOfCells=getNumberOfCells();
1099 const mcIdType *cPtr=_nodal_connec->begin();
1100 const mcIdType *icPtr=_nodal_connec_index->begin();
1102 mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1103 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1105 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1106 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1107 if(!cm.isQuadratic())
1109 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1110 types.insert(typ2); newConn->pushBackSilent(typ2);
1111 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1112 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1113 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1114 newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1115 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+1;
1116 newConnI->pushBackSilent(lastVal);
1117 ret->pushBackSilent(i);
1122 lastVal+=(icPtr[1]-icPtr[0]);
1123 newConnI->pushBackSilent(lastVal);
1124 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1127 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1128 coords=DataArrayDouble::Aggregate(coordsTmpSafe,tmp); conn=newConn.retn(); connI=newConnI.retn();
1133 * Implements \a conversionType 0 for meshes with meshDim = 3, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1134 * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1135 * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1137 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1139 MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1140 MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1141 return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1144 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1146 MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),desc2I(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1147 MCAuto<MEDCouplingUMesh> m2D=buildDescendingConnectivityGen<MinusOneSonsGeneratorBiQuadratic>(desc2,desc2I,tmp2,tmp3,MEDCouplingFastNbrer); tmp2=0; tmp3=0;
1148 MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),desc1I(DataArrayIdType::New()),tmp4(DataArrayIdType::New()),tmp5(DataArrayIdType::New());
1149 MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc1,desc1I,tmp4,tmp5); tmp4=0; tmp5=0;
1151 MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1152 MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1153 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(),ret2=DataArrayIdType::New(); ret->alloc(0,1); ret2->alloc(0,1);
1155 MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1156 const mcIdType *descPtr(desc1->begin()),*descIPtr(desc1I->begin()),*desc2Ptr(desc2->begin()),*desc2IPtr(desc2I->begin());
1157 DataArrayIdType *conn1D=0,*conn1DI=0,*conn2D=0,*conn2DI=0;
1158 std::set<INTERP_KERNEL::NormalizedCellType> types1D,types2D;
1159 DataArrayDouble *coordsTmp=0,*coordsTmp2=0;
1160 MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=DataArrayIdType::New(); ret1D->alloc(0,1);
1161 MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1162 MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1163 MCAuto<DataArrayIdType> ret2D=m2D->convertLinearCellsToQuadratic2D1(conn2D,conn2DI,coordsTmp2,types2D); ret2D=DataArrayIdType::New(); ret2D->alloc(0,1);
1164 MCAuto<DataArrayDouble> coordsTmp2Safe(coordsTmp2);
1165 MCAuto<DataArrayIdType> conn2DSafe(conn2D),conn2DISafe(conn2DI);
1166 const mcIdType *c1DPtr=conn1D->begin(),*c1DIPtr=conn1DI->begin(),*c2DPtr=conn2D->begin(),*c2DIPtr=conn2DI->begin();
1167 mcIdType nbOfCells=getNumberOfCells();
1168 const mcIdType *cPtr=_nodal_connec->begin();
1169 const mcIdType *icPtr=_nodal_connec_index->begin();
1171 mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1172 for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++,desc2IPtr++)
1174 INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1175 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1176 if(!cm.isQuadratic())
1178 INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1179 if(typ2==INTERP_KERNEL::NORM_ERROR)
1181 std::ostringstream oss; oss << "MEDCouplingUMesh::convertLinearCellsToQuadratic3D1 : On cell #" << i << " the linear cell type does not support advanced quadratization !";
1182 throw INTERP_KERNEL::Exception(oss.str());
1184 types.insert(typ2); newConn->pushBackSilent(typ2);
1185 newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1186 for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1187 newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1188 for(const mcIdType *d=desc2Ptr+desc2IPtr[0];d!=desc2Ptr+desc2IPtr[1];d++)
1190 mcIdType nodeId2=c2DPtr[c2DIPtr[(*d)+1]-1];
1191 mcIdType tmpPos=newConn->getNumberOfTuples();
1192 newConn->pushBackSilent(nodeId2);
1193 ret2D->pushBackSilent(nodeId2); ret1D->pushBackSilent(tmpPos);
1195 newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1196 lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+(desc2IPtr[1]-desc2IPtr[0])+1;
1197 newConnI->pushBackSilent(lastVal);
1198 ret->pushBackSilent(i);
1203 lastVal+=(icPtr[1]-icPtr[0]);
1204 newConnI->pushBackSilent(lastVal);
1205 newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1208 MCAuto<DataArrayIdType> diffRet2D=ret2D->getDifferentValues();
1209 MCAuto<DataArrayIdType> o2nRet2D=diffRet2D->invertArrayN2O2O2N(coordsTmp2Safe->getNumberOfTuples());
1210 coordsTmp2Safe=coordsTmp2Safe->selectByTupleId(diffRet2D->begin(),diffRet2D->end());
1211 MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1212 std::vector<const DataArrayDouble *> v(3); v[0]=coordsTmpSafe; v[1]=coordsTmp2Safe; v[2]=tmp;
1213 mcIdType *c=newConn->getPointer();
1214 const mcIdType *cI(newConnI->begin());
1215 for(const mcIdType *elt=ret1D->begin();elt!=ret1D->end();elt++)
1216 c[*elt]=o2nRet2D->getIJ(c[*elt],0)+offset;
1217 offset=coordsTmp2Safe->getNumberOfTuples();
1218 for(const mcIdType *elt=ret->begin();elt!=ret->end();elt++)
1219 c[cI[(*elt)+1]-1]+=offset;
1220 coords=DataArrayDouble::Aggregate(v); conn=newConn.retn(); connI=newConnI.retn();
1224 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1226 mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1227 const mcIdType *n2oPtr(n2o->begin());
1228 MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1229 skin->getReverseNodalConnectivity(revNodal,revNodalI);
1230 const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1231 const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1232 const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1233 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1234 mcIdType *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_POLYGON;
1235 if(nbOfNodesExpected<1)
1237 mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1238 *work++=n2oPtr[prevNode];
1239 for(mcIdType i=1;i<nbOfNodesExpected;i++)
1241 if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==3)
1243 std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1244 conn.erase(prevNode);
1247 mcIdType curNode(*(conn.begin()));
1248 *work++=n2oPtr[curNode];
1249 std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1250 shar.erase(prevCell);
1253 prevCell=*(shar.begin());
1257 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
1260 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
1263 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
1268 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1270 mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1271 mcIdType nbOfTurn(nbOfNodesExpected/2);
1272 const mcIdType *n2oPtr(n2o->begin());
1273 MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1274 skin->getReverseNodalConnectivity(revNodal,revNodalI);
1275 const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1276 const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1277 const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1278 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1279 mcIdType *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_QPOLYG;
1280 if(nbOfNodesExpected<1)
1282 mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1283 *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
1284 for(mcIdType i=1;i<nbOfTurn;i++)
1286 if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
1288 std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1289 conn.erase(prevNode);
1292 mcIdType curNode(*(conn.begin()));
1293 *work=n2oPtr[curNode];
1294 std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1295 shar.erase(prevCell);
1298 mcIdType curCell(*(shar.begin()));
1299 work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
1305 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
1308 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
1311 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
1316 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
1319 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
1320 std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
1321 int meshDim=(*it)->getMeshDimension();
1322 mcIdType nbOfCells=ToIdType((*it)->getNumberOfCells());
1323 mcIdType meshLgth=(*it++)->getNodalConnectivityArrayLen();
1324 for(;it!=a.end();it++)
1326 if(meshDim!=(*it)->getMeshDimension())
1327 throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
1328 nbOfCells+=ToIdType((*it)->getNumberOfCells());
1329 meshLgth+=(*it)->getNodalConnectivityArrayLen();
1331 std::vector<const MEDCouplingPointSet *> aps(a.size());
1332 std::copy(a.begin(),a.end(),aps.begin());
1333 MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1334 MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
1335 ret->setCoords(pts);
1336 MCAuto<DataArrayIdType> c=DataArrayIdType::New();
1337 c->alloc(meshLgth,1);
1338 mcIdType *cPtr=c->getPointer();
1339 MCAuto<DataArrayIdType> cI=DataArrayIdType::New();
1340 cI->alloc(nbOfCells+1,1);
1341 mcIdType *cIPtr=cI->getPointer();
1345 for(it=a.begin();it!=a.end();it++)
1347 mcIdType curNbOfCell=ToIdType((*it)->getNumberOfCells());
1348 const mcIdType *curCI=(*it)->_nodal_connec_index->begin();
1349 const mcIdType *curC=(*it)->_nodal_connec->begin();
1350 cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind2nd(std::plus<mcIdType>(),offset));
1351 for(mcIdType j=0;j<curNbOfCell;j++)
1353 const mcIdType *src=curC+curCI[j];
1355 for(;src!=curC+curCI[j+1];src++,cPtr++)
1363 offset+=curCI[curNbOfCell];
1364 offset2+=(*it)->getNumberOfNodes();
1367 ret->setConnectivity(c,cI,true);
1373 * \param [in] pt the start pointer (included) of the coordinates of the point
1374 * \param [in] cellIdsBg the start pointer (included) of cellIds
1375 * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1376 * \param [in] nc nodal connectivity
1377 * \param [in] ncI nodal connectivity index
1378 * \param [in,out] ret0 the min distance between \a this and the external input point
1379 * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1380 * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1382 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)
1385 ret0=std::numeric_limits<double>::max();
1386 for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1388 switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1390 case INTERP_KERNEL::NORM_TRI3:
1392 double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
1394 { ret0=tmp; cellId=*zeCell; }
1397 case INTERP_KERNEL::NORM_QUAD4:
1398 case INTERP_KERNEL::NORM_POLYGON:
1400 double tmp=INTERP_KERNEL::DistanceFromPtToPolygonInSpaceDim3(pt,nc+ncI[*zeCell]+1,nc+ncI[*zeCell+1],coords);
1402 { ret0=tmp; cellId=*zeCell; }
1406 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
1412 * \param [in] pt the start pointer (included) of the coordinates of the point
1413 * \param [in] cellIdsBg the start pointer (included) of cellIds
1414 * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1415 * \param [in] nc nodal connectivity
1416 * \param [in] ncI nodal connectivity index
1417 * \param [in,out] ret0 the min distance between \a this and the external input point
1418 * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1419 * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1421 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)
1424 ret0=std::numeric_limits<double>::max();
1425 for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1427 switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1429 case INTERP_KERNEL::NORM_SEG2:
1431 std::size_t uselessEntry=0;
1432 double tmp=INTERP_KERNEL::SquareDistanceFromPtToSegInSpaceDim2(pt,coords+2*nc[ncI[*zeCell]+1],coords+2*nc[ncI[*zeCell]+2],uselessEntry);
1435 { ret0=tmp; cellId=*zeCell; }
1439 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
1443 DataArrayIdType *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg(std::vector<bool>& fetched, const mcIdType *seedBg, const mcIdType *seedEnd, const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, mcIdType nbOfDepthPeeling, mcIdType& nbOfDepthPeelingPerformed)
1445 nbOfDepthPeelingPerformed=0;
1446 if(!seedBg || !seedEnd || !arrIn || !arrIndxIn)
1447 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : some input pointer is NULL !");
1448 mcIdType nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
1449 std::vector<bool> fetched2(nbOfTuples,false);
1451 for(const mcIdType *seedElt=seedBg;seedElt!=seedEnd;seedElt++,i++)
1453 if(*seedElt>=0 && *seedElt<nbOfTuples)
1454 { fetched[*seedElt]=true; fetched2[*seedElt]=true; }
1456 { 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()); }
1458 const mcIdType *arrInPtr=arrIn->begin();
1459 const mcIdType *arrIndxPtr=arrIndxIn->begin();
1460 mcIdType targetNbOfDepthPeeling=nbOfDepthPeeling!=-1?nbOfDepthPeeling:std::numeric_limits<mcIdType>::max();
1461 std::vector<mcIdType> idsToFetch1(seedBg,seedEnd);
1462 std::vector<mcIdType> idsToFetch2;
1463 std::vector<mcIdType> *idsToFetch=&idsToFetch1;
1464 std::vector<mcIdType> *idsToFetchOther=&idsToFetch2;
1465 while(!idsToFetch->empty() && nbOfDepthPeelingPerformed<targetNbOfDepthPeeling)
1467 for(std::vector<mcIdType>::const_iterator it=idsToFetch->begin();it!=idsToFetch->end();it++)
1468 for(const mcIdType *it2=arrInPtr+arrIndxPtr[*it];it2!=arrInPtr+arrIndxPtr[*it+1];it2++)
1470 { fetched[*it2]=true; fetched2[*it2]=true; idsToFetchOther->push_back(*it2); }
1471 std::swap(idsToFetch,idsToFetchOther);
1472 idsToFetchOther->clear();
1473 nbOfDepthPeelingPerformed++;
1475 mcIdType lgth=ToIdType(std::count(fetched2.begin(),fetched2.end(),true));
1477 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(lgth,1);
1478 mcIdType *retPtr=ret->getPointer();
1479 for(std::vector<bool>::const_iterator it=fetched2.begin();it!=fetched2.end();it++,i++)
1486 * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
1487 * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
1489 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, mcIdType nbOfNodesInCell, const mcIdType *conn, const double *coo, double *zipFrmt)
1493 for(int i=0;i<nbOfNodesInCell;i++)
1494 w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
1495 else if(spaceDim==2)
1497 for(int i=0;i<nbOfNodesInCell;i++)
1499 w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
1504 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
1508 * 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.
1509 * \param nbOfNodesPerLev in parameter that specifies the number of nodes of one slice of global dataset
1510 * \param isQuad specifies the policy of connectivity.
1511 * @ret in/out parameter in which the result will be append
1513 void MEDCouplingUMesh::AppendExtrudedCell(const mcIdType *connBg, const mcIdType *connEnd, mcIdType nbOfNodesPerLev, bool isQuad, std::vector<mcIdType>& ret)
1515 INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
1516 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
1517 ret.push_back(cm.getExtrudedType());
1518 mcIdType deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
1521 case INTERP_KERNEL::NORM_POINT1:
1523 ret.push_back(connBg[1]);
1524 ret.push_back(connBg[1]+nbOfNodesPerLev);
1527 case INTERP_KERNEL::NORM_SEG2:
1529 mcIdType conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
1530 ret.insert(ret.end(),conn,conn+4);
1533 case INTERP_KERNEL::NORM_SEG3:
1535 mcIdType conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
1536 ret.insert(ret.end(),conn,conn+8);
1539 case INTERP_KERNEL::NORM_QUAD4:
1541 mcIdType conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
1542 ret.insert(ret.end(),conn,conn+8);
1545 case INTERP_KERNEL::NORM_TRI3:
1547 mcIdType conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
1548 ret.insert(ret.end(),conn,conn+6);
1551 case INTERP_KERNEL::NORM_TRI6:
1553 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,
1554 connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
1555 ret.insert(ret.end(),conn,conn+15);
1558 case INTERP_KERNEL::NORM_QUAD8:
1561 connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
1562 connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
1563 connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
1565 ret.insert(ret.end(),conn,conn+20);
1568 case INTERP_KERNEL::NORM_POLYGON:
1570 std::back_insert_iterator< std::vector<mcIdType> > ii(ret);
1571 std::copy(connBg+1,connEnd,ii);
1573 std::reverse_iterator<const mcIdType *> rConnBg(connEnd);
1574 std::reverse_iterator<const mcIdType *> rConnEnd(connBg+1);
1575 std::transform(rConnBg,rConnEnd,ii,std::bind2nd(std::plus<mcIdType>(),deltaz));
1576 std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
1577 for(std::size_t i=0;i<nbOfRadFaces;i++)
1580 mcIdType conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
1581 std::copy(conn,conn+4,ii);
1586 throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
1592 * 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).
1593 * 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
1594 * with a plane. The result will be put in 'cut3DSuf' out parameter.
1595 * \param [in] cut3DCurve input parameter that gives for each 3DCurve cell if it owns fully to the plane or partially.
1596 * \param [out] nodesOnPlane returns all the nodes that are on the plane.
1597 * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
1598 * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
1599 * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
1600 * \param [in] nodalIndx3DCurve is the nodal connectivity index of 3D curve mesh.
1601 * \param [in] desc is the descending connectivity 3DSurf->3DCurve
1602 * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
1603 * \param [out] cut3DSurf input/output param.
1605 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
1606 const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
1607 const mcIdType *desc, const mcIdType *descIndx,
1608 std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf)
1610 std::set<mcIdType> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
1611 mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
1612 for(mcIdType i=0;i<nbOf3DSurfCell;i++)
1614 std::vector<mcIdType> res;
1615 mcIdType offset=descIndx[i];
1616 mcIdType nbOfSeg=descIndx[i+1]-offset;
1617 for(mcIdType j=0;j<nbOfSeg;j++)
1619 mcIdType edgeId=desc[offset+j];
1620 mcIdType status=cut3DCurve[edgeId];
1624 res.push_back(status);
1627 res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
1628 res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
1636 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1642 std::set<mcIdType> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
1643 std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<mcIdType> >(res));
1646 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1650 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
1655 {// case when plane is on a multi colinear edge of a polyhedron
1656 if(ToIdType(res.size())==2*nbOfSeg)
1658 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
1661 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
1669 * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
1670 * 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).
1671 * 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
1672 * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
1673 * \param cut3DSurf input parameter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
1674 * \param desc is the descending connectivity 3D->3DSurf
1675 * \param descIndx is the descending connectivity index 3D->3DSurf
1677 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1678 const mcIdType *desc, const mcIdType *descIndx,
1679 DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const
1681 checkFullyDefined();
1682 if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1683 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1684 const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1685 mcIdType nbOfCells=getNumberOfCells();
1686 for(mcIdType i=0;i<nbOfCells;i++)
1688 std::map<mcIdType, std::set<mcIdType> > m;
1689 mcIdType offset=descIndx[i];
1690 mcIdType nbOfFaces=descIndx[i+1]-offset;
1693 for(int j=0;j<nbOfFaces;j++)
1695 const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1696 if(p.first!=-1 && p.second!=-1)
1700 start=p.first; end=p.second;
1701 m[p.first].insert(p.second);
1702 m[p.second].insert(p.first);
1706 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
1707 mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1708 INTERP_KERNEL::AutoPtr<mcIdType> tmp=new mcIdType[sz];
1709 INTERP_KERNEL::NormalizedCellType cmsId;
1710 unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
1711 start=tmp[0]; end=tmp[nbOfNodesSon-1];
1712 for(unsigned k=0;k<nbOfNodesSon;k++)
1714 m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
1715 m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
1722 std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
1726 std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
1727 const std::set<mcIdType>& s=(*it).second;
1728 std::set<mcIdType> s2; s2.insert(prev);
1729 std::set<mcIdType> s3;
1730 std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
1733 mcIdType val=*s3.begin();
1734 conn.push_back(start);
1741 conn.push_back(end);
1744 nodalRes->insertAtTheEnd(conn.begin(),conn.end());
1745 nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
1746 cellIds->pushBackSilent(i);
1752 void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayIdType *nodalConnec, const DataArrayIdType *nodalConnecIndex)
1754 if(nodalConnec && nodalConnecIndex)
1757 const mcIdType *conn(nodalConnec->begin()),*connIndex(nodalConnecIndex->begin());
1758 mcIdType nbOfElem=ToIdType(nodalConnecIndex->getNbOfElems())-1;
1760 for(const mcIdType *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
1761 types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
1766 * This method expects that \a this a quadratic 1D, 2D or 3D mesh.
1767 * 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 )
1768 * 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.
1769 * So this method is non const because coordinates are modified.
1770 * 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.
1772 * \param [in] ratio - ratio of dilation
1773 * \param [in] nodeIdsBg - start (included) of input node list
1774 * \param [in] nodeIdsEnd - end (excluded) of input node list
1775 * \throw if there is a point in [\a nodeIdsBg, \a nodeIdsEnd ) that is a mid point of a seg3
1776 * \warning in case of throw the coordinates may be partially modified before the exception arises
1778 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodes(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1780 checkFullyDefined();
1781 int mdim(getMeshDimension());
1782 if(mdim==2 || mdim==3)
1784 MCAuto<MEDCouplingUMesh> edges;
1786 MCAuto<DataArrayIdType> a,b,c,d;
1787 edges=this->explodeIntoEdges(a,b,c,d);
1789 return edges->attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1792 return attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1793 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : not managed dimension ! Should be in [1,2,3] !");
1797 * \a this is expected to have meshdim==1.
1799 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodesUnderground(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1801 int spaceDim(getSpaceDimension());
1802 double *coords(getCoords()->getPointer());
1803 auto nbNodes(getNumberOfNodes());
1804 std::size_t nbCells(getNumberOfCells());
1805 std::vector<bool> fastFinder(nbNodes,false);
1806 for(auto work=nodeIdsBg;work!=nodeIdsEnd;work++)
1807 if(*work>=0 && *work<nbNodes)
1808 fastFinder[*work]=true;
1809 MCAuto<DataArrayIdType> cellsIds(getCellIdsLyingOnNodes(nodeIdsBg,nodeIdsEnd,false));
1810 const mcIdType *nc(_nodal_connec->begin()),*nci(_nodal_connec_index->begin());
1811 for(std::size_t cellId=0;cellId<nbCells;cellId++,nci++)
1813 const mcIdType *isSelected(std::find_if(nc+nci[0]+1,nc+nci[1],[&fastFinder](mcIdType v) { return fastFinder[v]; }));
1814 if(isSelected!=nc+nci[1])
1816 if((INTERP_KERNEL::NormalizedCellType)nc[nci[0]]==INTERP_KERNEL::NORM_SEG3 && nci[1]-nci[0]==4)
1818 bool aa(fastFinder[nc[*nci+1]]),bb(fastFinder[nc[*nci+2]]),cc(fastFinder[nc[*nci+3]]);
1823 auto ptToMove(nc[*nci+3]);
1824 auto attractor(aa?nc[*nci+1]:nc[*nci+2]),endPt(aa?nc[*nci+2]:nc[*nci+1]);
1825 std::transform(coords+spaceDim*attractor,coords+spaceDim*(attractor+1),coords+spaceDim*endPt,
1826 coords+spaceDim*ptToMove,[ratio](const double& stPt2, const double& endPt2) { return stPt2+ratio*(endPt2-stPt2); });
1829 continue;//both 2 boundary nodes of current seg3 are un nodeIds input list -> skip it.
1833 std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " has a mid point " << nc[*nci+3] << " ! This node is in input list !";
1834 throw INTERP_KERNEL::Exception(oss.str());
1839 std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " sharing one of the input nodes list its geo type is NOT SEG3 !";
1840 throw INTERP_KERNEL::Exception(oss.str());