+/*!
+ * This method expects that \a this a 3D mesh (spaceDim=3 and meshDim=3) with all coordinates and connectivities set.
+ * All cells in \a this are expected to be linear 3D cells.
+ * This method will split **all** 3D cells in \a this into INTERP_KERNEL::NORM_TETRA4 cells and put them in the returned mesh.
+ * It leads to an increase to number of cells.
+ * This method contrary to MEDCouplingUMesh::simplexize can append coordinates in \a this to perform its work.
+ * The \a nbOfAdditionalPoints returned value informs about it. If > 0, the coordinates array in returned mesh will have \a nbOfAdditionalPoints
+ * more tuples (nodes) than in \a this. Anyway, all the nodes in \a this (with the same order) will be in the returned mesh.
+ *
+ * \param [in] policy - the policy of splitting that must be in (PLANAR_FACE_5, PLANAR_FACE_6, GENERAL_24, GENERAL_48). The policy will be used only for INTERP_KERNEL::NORM_HEXA8 cells.
+ * For all other cells, the splitting policy will be ignored.
+ * \param [out] nbOfAdditionalPoints - number of nodes added to \c this->_coords. If > 0 a new coordinates object will be constructed result of the aggregation of the old one and the new points added.
+ * \param [out] n2oCells - A new instance of DataArrayInt holding, for each new cell,
+ * an id of old cell producing it. The caller is to delete this array using
+ * decrRef() as it is no more needed.
+ * \return MEDCoupling1SGTUMesh * - the mesh containing only INTERP_KERNEL::NORM_TETRA4 cells.
+ *
+ * \warning This method operates on each cells in this independantly ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output
+ * the policy PLANAR_FACE_6 should be used on a mesh sorted with MEDCoupling1SGTUMesh::sortHexa8EachOther.
+ *
+ * \throw If \a this is not a 3D mesh (spaceDim==3 and meshDim==3).
+ * \throw If \a this is not fully constituted with linear 3D cells.
+ * \sa MEDCouplingUMesh::simplexize, MEDCoupling1SGTUMesh::sortHexa8EachOther
+ */
+MEDCoupling1SGTUMesh *MEDCouplingUMesh::tetrahedrize(int policy, DataArrayInt *& n2oCells, int& nbOfAdditionalPoints) const
+{
+ INTERP_KERNEL::SplittingPolicy pol((INTERP_KERNEL::SplittingPolicy)policy);
+ checkConnectivityFullyDefined();
+ if(getMeshDimension()!=3 || getSpaceDimension()!=3)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tetrahedrize : only available for mesh with meshdim == 3 and spacedim == 3 !");
+ int nbOfCells(getNumberOfCells()),nbNodes(getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret0(MEDCoupling1SGTUMesh::New(getName(),INTERP_KERNEL::NORM_TETRA4));
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfCells,1);
+ int *retPt(ret->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn(DataArrayInt::New()); newConn->alloc(0,1);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addPts(DataArrayDouble::New()); addPts->alloc(0,1);
+ const int *oldc(_nodal_connec->begin());
+ const int *oldci(_nodal_connec_index->begin());
+ const double *coords(_coords->begin());
+ for(int i=0;i<nbOfCells;i++,oldci++,retPt++)
+ {
+ std::vector<int> a; std::vector<double> b;
+ INTERP_KERNEL::SplitIntoTetras(pol,(INTERP_KERNEL::NormalizedCellType)oldc[oldci[0]],oldc+oldci[0]+1,oldc+oldci[1],coords,a,b);
+ std::size_t nbOfTet(a.size()/4); *retPt=(int)nbOfTet;
+ const int *aa(&a[0]);
+ if(!b.empty())
+ {
+ for(std::vector<int>::iterator it=a.begin();it!=a.end();it++)
+ if(*it<0)
+ *it=(-(*(it))-1+nbNodes);
+ addPts->insertAtTheEnd(b.begin(),b.end());
+ nbNodes+=(int)b.size()/3;
+ }
+ for(std::size_t j=0;j<nbOfTet;j++,aa+=4)
+ newConn->insertAtTheEnd(aa,aa+4);
+ }
+ if(!addPts->empty())
+ {
+ addPts->rearrange(3);
+ nbOfAdditionalPoints=addPts->getNumberOfTuples();
+ addPts=DataArrayDouble::Aggregate(getCoords(),addPts);
+ ret0->setCoords(addPts);
+ }
+ else
+ {
+ nbOfAdditionalPoints=0;
+ ret0->setCoords(getCoords());
+ }
+ ret0->setNodalConnectivity(newConn);
+ //
+ ret->computeOffsets2();
+ n2oCells=ret->buildExplicitArrOfSliceOnScaledArr(0,nbOfCells,1);
+ return ret0.retn();
+}
+
+/*!
+ * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
+ *
+ * \sa MEDCouplingUMesh::split2DCells
+ */
+void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
+{
+ checkConnectivityFullyDefined();
+ int ncells(getNumberOfCells()),lgthToReach(getMeshLength()+subNodesInSeg->getNumberOfTuples());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
+ const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
+ int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
+ int prevPosOfCi(ciPtr[0]);
+ for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
+ {
+ int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
+ *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
+ for(int j=0;j<sz;j++)
+ {
+ int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
+ for(int k=0;k<sz2;k++)
+ *cPtr++=subPtr[offset2+k];
+ if(j!=sz-1)
+ *cPtr++=oldConn[prevPosOfCi+j+2];
+ deltaSz+=sz2;
+ }
+ prevPosOfCi=ciPtr[1];
+ ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
+ }
+ if(c->end()!=cPtr)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
+ _nodal_connec->decrRef();
+ _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
+}
+
+int internalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
+{
+ if(id!=-1)
+ return id;
+ else
+ {
+ int ret(nodesCnter++);
+ double newPt[2];
+ e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
+ addCoo.insertAtTheEnd(newPt,newPt+2);
+ return ret;
+ }
+}
+
+/*!
+ * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
+ *
+ * \return int - the number of new nodes created.
+ * \sa MEDCouplingUMesh::split2DCells
+ */
+int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
+{
+ checkCoherency();
+ int ncells(getNumberOfCells()),lgthToReach(getMeshLength()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
+ const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
+ const int *midPtr(mid->begin()),*midIPtr(midI->begin());
+ const double *oldCoordsPtr(getCoords()->begin());
+ int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
+ int prevPosOfCi(ciPtr[0]);
+ for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
+ {
+ int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
+ for(int j=0;j<sz;j++)
+ { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
+ *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
+ for(int j=0;j<sz;j++)//loop over subedges of oldConn
+ {
+ int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
+ if(sz2==0)
+ {
+ if(j<sz-1)
+ cPtr[1]=oldConn[prevPosOfCi+2+j];
+ cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
+ continue;
+ }
+ std::vector<INTERP_KERNEL::Node *> ns(3);
+ ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
+ ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
+ ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
+ for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
+ {
+ cPtr[1]=subPtr[offset2+k];
+ cPtr[deltaSz]=internalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
+ }
+ int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
+ if(j!=sz-1)
+ { cPtr[1]=tmpEnd; }
+ cPtr[deltaSz]=internalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
+ }
+ prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
+ ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
+ }
+ if(c->end()!=cPtr)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
+ _nodal_connec->decrRef();
+ _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
+ addCoo->rearrange(2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
+ setCoords(coo);
+ return addCoo->getNumberOfTuples();
+}
+