+ ret->setCoords(getCoords());
+ ret->setNodalConnectivity(c);
+ return ret.retn();
+}
+
+/*!
+ * This method for each cell in \a this the triangle height for each edge in a newly allocated/created array instance.
+ *
+ * \return DataArrayDouble * - a newly allocated instance with this->getNumberOfCells() tuples and 3 components storing for each cell in \a this the corresponding height.
+ * \throw If \a this is not a mesh containing only NORM_TRI3 cells.
+ * \throw If \a this is not properly allocated.
+ * \throw If spaceDimension is not in 2 or 3.
+ */
+MCAuto<DataArrayDouble> MEDCoupling1SGTUMesh::computeTriangleHeight() const
+{
+ checkConsistencyLight();
+ const INTERP_KERNEL::CellModel& cm(getCellModel());
+ if(cm.getEnum()!=INTERP_KERNEL::NORM_TRI3)
+ THROW_IK_EXCEPTION("MEDCoupling1SGTUMesh::computeTriangleHeight : this method can be applied only on TRI3 mesh !");
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ mcIdType nbTri3( getNumberOfCells() );
+ const double *coordPtr( this->getCoords()->begin() );
+ const mcIdType *inConnPtr(getNodalConnectivity()->begin());
+ ret->alloc(nbTri3,3);
+ double *retPtr( ret->getPointer() );
+ switch( this->getSpaceDimension())
+ {
+ case 2:
+ {
+ constexpr unsigned SPACEDIM = 2;
+ for(mcIdType iCell = 0 ; iCell < nbTri3 ; ++iCell)
+ {
+ INTERP_KERNEL::ComputeTriangleHeight<SPACEDIM>(coordPtr + SPACEDIM*inConnPtr[3*iCell+0], coordPtr + SPACEDIM*inConnPtr[3*iCell+1], coordPtr + SPACEDIM*inConnPtr[3*iCell+2],retPtr+3*iCell);
+ }
+ break;
+ }
+ case 3:
+ {
+ constexpr unsigned SPACEDIM = 3;
+ for(mcIdType iCell = 0 ; iCell < nbTri3 ; ++iCell)
+ {
+ INTERP_KERNEL::ComputeTriangleHeight<SPACEDIM>(coordPtr + SPACEDIM*inConnPtr[3*iCell+0], coordPtr + SPACEDIM*inConnPtr[3*iCell+1], coordPtr + SPACEDIM*inConnPtr[3*iCell+2],retPtr+3*iCell);
+ }
+ break;
+ }
+ default:
+ THROW_IK_EXCEPTION("MEDCoupling1SGTUMesh::computeTriangleHeight : only spacedim in [2,3] supported !");
+ }
+ return ret;
+}
+
+/*!
+ * This method starts from an unstructured mesh that hides in reality a cartesian mesh.
+ * If it is not the case, an exception will be thrown.
+ * This method returns three objects : The cartesian mesh geometrically equivalent to \a this (within a precision of \a eps) and a permutation of cells
+ * and a permutation of nodes.
+ *
+ * - this[cellPerm[i]]=ret[i]
+ *
+ * \param [out] cellPerm the permutation array of size \c this->getNumberOfCells()
+ * \param [out] nodePerm the permutation array of size \c this->getNumberOfNodes()
+ * \return MEDCouplingCMesh * - a newly allocated mesh that is the result of the structurization of \a this.
+ */
+MEDCouplingCMesh *MEDCoupling1SGTUMesh::structurizeMe(DataArrayIdType *& cellPerm, DataArrayIdType *& nodePerm, double eps) const
+{
+ checkConsistencyLight();
+ int spaceDim(getSpaceDimension()),meshDim(getMeshDimension()); mcIdType nbNodes(getNumberOfNodes());
+ if(MEDCouplingStructuredMesh::GetGeoTypeGivenMeshDimension(meshDim)!=getCellModelEnum())
+ throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::structurizeMe : the unique geo type in this is not compatible with the geometric type regarding mesh dimension !");
+ MCAuto<MEDCouplingCMesh> cm(MEDCouplingCMesh::New());
+ for(int i=0;i<spaceDim;i++)
+ {
+ std::vector<std::size_t> tmp(1,i);
+ MCAuto<DataArrayDouble> elt(static_cast<DataArrayDouble*>(getCoords()->keepSelectedComponents(tmp)));
+ elt=elt->getDifferentValues(eps);
+ elt->sort(true);
+ cm->setCoordsAt(i,elt);
+ }
+ if(nbNodes!=cm->getNumberOfNodes())
+ throw INTERP_KERNEL::Exception("MEDCoupling1SGTUMesh::structurizeMe : considering the number of nodes after split per components in space this can't be a cartesian mesh ! Maybe your epsilon parameter is invalid ?");
+ try
+ { cm->copyTinyInfoFrom(this); }
+ catch(INTERP_KERNEL::Exception&) { }
+ MCAuto<MEDCouplingUMesh> um(cm->buildUnstructured()),self(buildUnstructured());
+ self->checkGeoEquivalWith(um,12,eps,cellPerm,nodePerm);
+ return cm.retn();
+}
+
+/// @cond INTERNAL
+
+bool UpdateHexa8Cell(int validAxis, mcIdType neighId, const mcIdType *validConnQuad4NeighSide, mcIdType *allFacesNodalConn, mcIdType *myNeighbours)
+{
+ static const int TAB[48]={
+ 0,1,2,3,4,5,6,7,//0
+ 4,7,6,5,0,3,2,1,//1
+ 0,3,7,4,1,2,6,5,//2
+ 4,0,3,7,5,1,2,6,//3
+ 5,1,0,4,6,2,3,7,//4
+ 3,7,4,0,2,6,5,1 //5
+ };
+ static const int TAB2[6]={0,0,3,3,3,3};
+ if(myNeighbours[validAxis]==neighId && allFacesNodalConn[4*validAxis+0]==validConnQuad4NeighSide[TAB2[validAxis]])
+ return true;
+ mcIdType oldAxis(ToIdType(std::distance(myNeighbours,std::find(myNeighbours,myNeighbours+6,neighId))));
+ std::size_t pos(std::distance(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,std::find(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS,MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS+6,oldAxis)));
+ std::size_t pos0(pos/2),pos1(pos%2);
+ int oldAxisOpp(MEDCoupling1SGTUMesh::HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2]);
+ mcIdType oldConn[8],myConn2[8]={-1,-1,-1,-1,-1,-1,-1,-1},myConn[8],edgeConn[2],allFacesTmp[24],neighTmp[6];
+ oldConn[0]=allFacesNodalConn[0]; oldConn[1]=allFacesNodalConn[1]; oldConn[2]=allFacesNodalConn[2]; oldConn[3]=allFacesNodalConn[3];
+ oldConn[4]=allFacesNodalConn[4]; oldConn[5]=allFacesNodalConn[7]; oldConn[6]=allFacesNodalConn[6]; oldConn[7]=allFacesNodalConn[5];
+ const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_HEXA8));
+ for(int i=0;i<4;i++)
+ myConn2[i]=validConnQuad4NeighSide[(4-i+TAB2[validAxis])%4];
+ for(int i=0;i<4;i++)
+ {
+ mcIdType nodeId(myConn2[i]);//the node id for which the opposite one will be found
+ bool found(false);
+ INTERP_KERNEL::NormalizedCellType typeOfSon;
+ for(int j=0;j<12 && !found;j++)
+ {
+ cm.fillSonEdgesNodalConnectivity3D(j,oldConn,-1,edgeConn,typeOfSon);
+ if(edgeConn[0]==nodeId || edgeConn[1]==nodeId)
+ {
+ if(std::find(allFacesNodalConn+4*oldAxisOpp,allFacesNodalConn+4*oldAxisOpp+4,edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0])!=allFacesNodalConn+4*oldAxisOpp+4)
+ {
+ myConn2[i+4]=edgeConn[0]==nodeId?edgeConn[1]:edgeConn[0];
+ found=true;
+ }
+ }
+ }
+ if(!found)
+ throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error !");
+ }
+ const int *myTab(TAB+8*validAxis);
+ for(int i=0;i<8;i++)
+ myConn[i]=myConn2[myTab[i]];
+ for(int i=0;i<6;i++)
+ {
+ cm.fillSonCellNodalConnectivity(i,myConn,allFacesTmp+4*i);
+ std::set<mcIdType> s(allFacesTmp+4*i,allFacesTmp+4*i+4);
+ bool found(false);
+ for(int j=0;j<6 && !found;j++)
+ {
+ std::set<mcIdType> s1(allFacesNodalConn+4*j,allFacesNodalConn+4*j+4);
+ if(s==s1)
+ {
+ neighTmp[i]=myNeighbours[j];
+ found=true;
+ }
+ }
+ if(!found)
+ throw INTERP_KERNEL::Exception("UpdateHexa8Cell : Internal Error #2 !");
+ }
+ std::copy(allFacesTmp,allFacesTmp+24,allFacesNodalConn);
+ std::copy(neighTmp,neighTmp+6,myNeighbours);
+ return false;
+}
+
+/// @endcond
+
+/*!
+ * This method expects the \a this contains NORM_HEXA8 cells only. This method will sort each cells in \a this so that their numbering was
+ * homogeneous. If it succeeds the result of MEDCouplingUMesh::tetrahedrize will return a conform mesh.
+ *
+ * \return DataArrayIdType * - a newly allocated array (to be managed by the caller) containing renumbered cell ids.
+ *
+ * \throw If \a this is not a mesh containing only NORM_HEXA8 cells.
+ * \throw If \a this is not properly allocated.
+ * \sa MEDCouplingUMesh::tetrahedrize, MEDCouplingUMesh::simplexize.
+ */
+DataArrayIdType *MEDCoupling1SGTUMesh::sortHexa8EachOther()
+{
+ MCAuto<MEDCoupling1SGTUMesh> quads(explodeEachHexa8To6Quad4());//checks that only hexa8
+ mcIdType nbHexa8=getNumberOfCells();
+ mcIdType *cQuads(quads->getNodalConnectivity()->getPointer());
+ MCAuto<DataArrayIdType> neighOfQuads(DataArrayIdType::New()); neighOfQuads->alloc(nbHexa8*6,1); neighOfQuads->fillWithValue(-1);
+ mcIdType *ptNeigh(neighOfQuads->getPointer());
+ {//neighOfQuads tells for each face of each Quad8 which cell (if!=-1) is connected to this face.
+ MCAuto<MEDCouplingUMesh> quadsTmp(quads->buildUnstructured());
+ MCAuto<DataArrayIdType> ccSafe,cciSafe;
+ DataArrayIdType *cc(0),*cci(0);
+ quadsTmp->findCommonCells(3,0,cc,cci);
+ ccSafe=cc; cciSafe=cci;
+ const mcIdType *ccPtr(ccSafe->begin());
+ mcIdType nbOfPair=cci->getNumberOfTuples()-1;
+ for(mcIdType i=0;i<nbOfPair;i++)
+ { ptNeigh[ccPtr[2*i+0]]=ccPtr[2*i+1]/6; ptNeigh[ccPtr[2*i+1]]=ccPtr[2*i+0]/6; }
+ }
+ MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+ std::vector<bool> fetched(nbHexa8,false);
+ std::vector<bool>::iterator it(std::find(fetched.begin(),fetched.end(),false));
+ while(it!=fetched.end())//it will turns as time as number of connected zones
+ {
+ mcIdType cellId(ToIdType(std::distance(fetched.begin(),it)));//it is the seed of the connected zone.
+ std::set<mcIdType> s; s.insert(cellId);//s contains already organized.
+ while(!s.empty())
+ {
+ std::set<mcIdType> sNext;
+ for(std::set<mcIdType>::const_iterator it0=s.begin();it0!=s.end();it0++)
+ {
+ fetched[*it0]=true;
+ mcIdType *myNeighb(ptNeigh+6*(*it0));
+ for(int i=0;i<6;i++)
+ {
+ if(myNeighb[i]!=-1 && !fetched[myNeighb[i]])
+ {
+ std::size_t pos(std::distance(HEXA8_FACE_PAIRS,std::find(HEXA8_FACE_PAIRS,HEXA8_FACE_PAIRS+6,i)));
+ std::size_t pos0(pos/2),pos1(pos%2);
+ if(!UpdateHexa8Cell(HEXA8_FACE_PAIRS[2*pos0+(pos1+1)%2],*it0,cQuads+6*4*(*it0)+4*i,cQuads+6*4*myNeighb[i],ptNeigh+6*myNeighb[i]))
+ ret->pushBackSilent(myNeighb[i]);
+ fetched[myNeighb[i]]=true;
+ sNext.insert(myNeighb[i]);
+ }
+ }
+ }
+ s=sNext;
+ }
+ it=std::find(fetched.begin(),fetched.end(),false);
+ }
+ if(!ret->empty())
+ {
+ mcIdType *conn(getNodalConnectivity()->getPointer());
+ for(const mcIdType *pt=ret->begin();pt!=ret->end();pt++)
+ {
+ mcIdType cellId(*pt);
+ conn[8*cellId+0]=cQuads[24*cellId+0]; conn[8*cellId+1]=cQuads[24*cellId+1]; conn[8*cellId+2]=cQuads[24*cellId+2]; conn[8*cellId+3]=cQuads[24*cellId+3];
+ conn[8*cellId+4]=cQuads[24*cellId+4]; conn[8*cellId+5]=cQuads[24*cellId+7]; conn[8*cellId+6]=cQuads[24*cellId+6]; conn[8*cellId+7]=cQuads[24*cellId+5];
+ }
+ declareAsNew();
+ }
+ return ret.retn();