From: abn Date: Wed, 4 Oct 2017 09:15:08 +0000 (+0200) Subject: Bug fix in simplifyPolyhedronCells(). Thanks Antoine G. ! X-Git-Tag: V8_4_0rc1~1 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=e6b37f140ed2fa71a16d52adab59ae1bc3d24244;p=tools%2Fmedcoupling.git Bug fix in simplifyPolyhedronCells(). Thanks Antoine G. ! --- diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 697b876ab..87871d83b 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -1308,12 +1308,14 @@ void MEDCouplingUMesh::simplifyPolyhedra(double eps) connINew->alloc(nbOfCells+1,1); int *connINewPtr=connINew->getPointer(); *connINewPtr++=0; MCAuto connNew=DataArrayInt::New(); connNew->alloc(0,1); + MCAuto E_Fi(DataArrayInt::New()), E_F(DataArrayInt::New()), F_Ei(DataArrayInt::New()), F_E(DataArrayInt::New()); + MCAuto m_faces(buildDescendingConnectivity(E_F, E_Fi, F_E, F_Ei)); bool changed=false; for(int i=0;igetIJ(index + 1, 0) - E_Fi->getIJ(index, 0); MCAuto v=DataArrayDouble::New(); v->alloc(nbFaces,3); double *vPtr=v->getPointer(); - MCAuto p=DataArrayDouble::New(); p->alloc(nbFaces,1); + MCAuto p=DataArrayDouble::New(); p->alloc(nbFaces,2); double *pPtr=p->getPointer(); - const int *stFaceConn=begin+1; + int *e_fi = E_Fi->getPointer(), *e_f = E_F->getPointer(), *f_ei = F_Ei->getPointer(), *f_e = F_E->getPointer(); + const int *f_idx = faces->getNodalConnectivityIndex()->getPointer(), *f_cnn = faces->getNodalConnectivity()->getPointer(); for(int i=0;ibegin(),stFaceConn,endFaceConn,vPtr,pPtr); - stFaceConn=endFaceConn+1; + int face = e_f[e_fi[index] + i]; + ComputeVecAndPtOfFace(eps, coords->begin(), f_cnn + f_idx[face] + 1, f_cnn + f_idx[face + 1], vPtr, pPtr); + // to differentiate faces going to different cells: + pPtr++, *pPtr = 0; + for (int j = f_ei[face]; j < f_ei[face + 1]; j++) + *pPtr += f_e[j]; } pPtr=p->getPointer(); vPtr=v->getPointer(); DataArrayInt *comm1=0,*commI1=0; v->findCommonTuples(eps,-1,comm1,commI1); + for (int i = 0; i < nbFaces; i++) + if (comm1->findIdFirstEqual(i) < 0) + { + comm1->pushBackSilent(i); + commI1->pushBackSilent(comm1->getNumberOfTuples()); + } MCAuto comm1Auto(comm1),commI1Auto(commI1); const int *comm1Ptr=comm1->begin(); const int *commI1Ptr=commI1->begin(); int nbOfGrps1=commI1Auto->getNumberOfTuples()-1; res->pushBackSilent((int)INTERP_KERNEL::NORM_POLYHED); // - MCAuto mm=MEDCouplingUMesh::New("",3); - mm->setCoords(const_cast(coords)); mm->allocateCells(1); mm->insertNextCell(INTERP_KERNEL::NORM_POLYHED,(int)std::distance(begin+1,end),begin+1); - mm->finishInsertingCells(); - // for(int i=0;i tmpgrp2=p->selectByTupleId(comm1Ptr+commI1Ptr[i],comm1Ptr+commI1Ptr[i+1]); DataArrayInt *comm2=0,*commI2=0; tmpgrp2->findCommonTuples(eps,-1,comm2,commI2); + for (int j = 0; j < commI1Ptr[i+1] - commI1Ptr[i]; j++) + if (comm2->findIdFirstEqual(j) < 0) + { + comm2->pushBackSilent(j); + commI2->pushBackSilent(comm2->getNumberOfTuples()); + } MCAuto comm2Auto(comm2),commI2Auto(commI2); const int *comm2Ptr=comm2->begin(); const int *commI2Ptr=commI2->begin(); int nbOfGrps2=commI2Auto->getNumberOfTuples()-1; for(int j=0;jinsertAtTheEnd(begin,end); + int face = e_f[e_fi[index] + comm1Ptr[commI1Ptr[i] + comm2Ptr[commI2Ptr[j]]]]; //hmmm + res->insertAtTheEnd(f_cnn + f_idx[face] + 1, f_cnn + f_idx[face + 1]); res->pushBackSilent(-1); } else @@ -7069,13 +7085,12 @@ void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble int pointId=comm1Ptr[commI1Ptr[i]+comm2Ptr[commI2Ptr[j]]]; MCAuto ids2=comm2->selectByTupleIdSafeSlice(commI2Ptr[j],commI2Ptr[j+1],1); ids2->transformWithIndArr(comm1Ptr+commI1Ptr[i],comm1Ptr+commI1Ptr[i+1]); - DataArrayInt *tmp0=DataArrayInt::New(),*tmp1=DataArrayInt::New(),*tmp2=DataArrayInt::New(),*tmp3=DataArrayInt::New(); - MCAuto mm2=mm->buildDescendingConnectivity(tmp0,tmp1,tmp2,tmp3); tmp0->decrRef(); tmp1->decrRef(); tmp2->decrRef(); tmp3->decrRef(); - MCAuto mm3=static_cast(mm2->buildPartOfMySelf(ids2->begin(),ids2->end(),true)); + ids2->transformWithIndArr(e_f + e_fi[index], e_f + e_fi[index + 1]); + MCAuto mm3=static_cast(faces->buildPartOfMySelf(ids2->begin(),ids2->end(),true)); MCAuto idsNodeTmp=mm3->zipCoordsTraducer(); MCAuto idsNode=idsNodeTmp->invertArrayO2N2N2O(mm3->getNumberOfNodes()); const int *idsNodePtr=idsNode->begin(); - double center[3]; center[0]=pPtr[pointId]*vPtr[3*vecId]; center[1]=pPtr[pointId]*vPtr[3*vecId+1]; center[2]=pPtr[pointId]*vPtr[3*vecId+2]; + double center[3]; center[0]=pPtr[2*pointId]*vPtr[3*vecId]; center[1]=pPtr[2*pointId]*vPtr[3*vecId+1]; center[2]=pPtr[2*pointId]*vPtr[3*vecId+2]; double vec[3]; vec[0]=vPtr[3*vecId+1]; vec[1]=-vPtr[3*vecId]; vec[2]=0.; double norm=vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]; if(std::abs(norm)>eps) diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 07cec54b9..fec5af450 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -251,7 +251,8 @@ namespace MEDCoupling MEDCOUPLING_EXPORT static void CorrectExtrudedStaticCell(int *begin, int *end); MEDCOUPLING_EXPORT static bool IsTetra4WellOriented(const int *begin, const int *end, const double *coords); MEDCOUPLING_EXPORT static bool IsPyra5WellOriented(const int *begin, const int *end, const double *coords); - MEDCOUPLING_EXPORT static void SimplifyPolyhedronCell(double eps, const DataArrayDouble *coords, const int *begin, const int *end, DataArrayInt *res); + MEDCOUPLING_EXPORT static void SimplifyPolyhedronCell(double eps, const DataArrayDouble *coords, int index, DataArrayInt *res, MEDCouplingUMesh *faces, + DataArrayInt *E_Fi, DataArrayInt *E_F, DataArrayInt *F_Ei, DataArrayInt *F_E); MEDCOUPLING_EXPORT static void ComputeVecAndPtOfFace(double eps, const double *coords, const int *begin, const int *end, double *v, double *p); MEDCOUPLING_EXPORT static void TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords); MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index 6035e19fb..c3a254a15 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -17,7 +17,6 @@ # # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com # - from MEDCoupling import * import unittest from math import pi,e,sqrt,cos,sin @@ -4304,7 +4303,7 @@ class MEDCouplingBasicsTest5(unittest.TestCase): mea=fieldOnCell.getMesh().getMeasureField(True).getArray(); mea.rearrange(48); mea2 = mea.sumPerTuple() self.assertTrue(mea2.isEqual(meaRef2,1e-9)) pass - + def testVoronoi3DSurf_1(self): tmp=MEDCouplingCMesh("mesh") arr=DataArrayDouble(5) ; arr.iota() @@ -4474,7 +4473,7 @@ class MEDCouplingBasicsTest5(unittest.TestCase): self.assertTrue(f3.getMesh().getMeasureField(False).getArray().isEqual(ref,1e-12)) self.assertTrue(f3.getArray().isEqual(DataArrayDouble([0,1,2,3]),1e-12)) pass - + def testVoronoi3D_4(self): """Idem testVoronoi3D_3 except that here quadratic cells are considered""" coo=DataArrayDouble([0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,0.5,0.0,0.0,0.0,0.5,0.0,0.5,0.5,0.5,0.5,0.0,0.5,0.0,0.0,0.5,0.0,0.5],10,3) @@ -4767,7 +4766,22 @@ class MEDCouplingBasicsTest5(unittest.TestCase): m.insertNextCell(NORM_POLYGON,[0,1,2,3,4,5]) self.assertTrue(m.computePlaneEquationOf3DFaces().isEqual(DataArrayDouble([0,1,0,-1],1,4),1e-12)) pass - + + def testSimplifyPolyhedra(self): + mesh = MEDCouplingUMesh('mesh', 3) + coo = DataArrayDouble([(-0.01225,-0.0212176,0.02),(-0.00634107,-0.0236652,0.02),(1.50019e-18,-0.0245,0.02),(0.00634107,-0.0236652,0.02),(0.01225,-0.0212176,0.02),(-0.0153864,-0.02665,0),(-0.00714085,-0.02665,0),(1.63184e-18,-0.02665,0),(0.00714085,-0.02665,0),(0.0153864,-0.02665,0),(-0.00714085,-0.02665,0.0101475),(1.63184e-18,-0.02665,0.013145),(0.00714085,-0.02665,0.0101475),(-0.013,-0.0225167,0.02),(-0.0067293,-0.0251141,0.02),(1.59204e-18,-0.026,0.02),(0.0067293,-0.0251141,0.02),(0.013,-0.0225167,0.02),(-0.0161658,-0.028,0),(-0.00750258,-0.028,0),(1.71451e-18,-0.028,0),(0.00750258,-0.028,0),(0.0161658,-0.028,0),(-0.00750258,-0.028,0.0105625),(1.71451e-18,-0.028,0.0136825),(0.00750258,-0.028,0.0105625)]) + mesh.setCoords(coo) + c = DataArrayInt([31, 13, 14, 15, 16, 17, 4, 3, 2, 1, 0, -1, 18, 5, 6, 7, 8, 9, 22, 21, 20, 19, -1, 19, 23, 18, -1, 23, 14, 13, 18, -1, 20, 24, 23, 19, -1, 24, 15, 14, 23, -1, 21, 25, 24, 20, -1, 25, 16, 15, 24, -1, 22, 25, 21, -1, 22, 17, 16, 25, -1, 9, 4, 17, 22, -1, 8, 12, 9, -1, 12, 3, 4, 9, -1, 7, 11, 12, 8, -1, 11, 2, 3, 12, -1, 6, 10, 11, 7, -1, 10, 1, 2, 11, -1, 5, 10, 6, -1, 5, 0, 1, 10, -1, 18, 13, 0, 5]) + cI = DataArrayInt([0, 108]) + mesh.setConnectivity(c, cI) + mesh.simplifyPolyhedra(1.0e-8) + c, cI = mesh.getNodalConnectivity(), mesh.getNodalConnectivityIndex() + tgt_c = DataArrayInt([31, 23, 18, 19, 20, 21, 22, 25, 24, -1, 12, 9, 8, 7, 6, 5, 10, 11, -1, 13, 14, 15, 16, 17, 4, 3, 2, 1, 0, -1, 18, 5, 6, 7, 8, 9, 22, 21, 20, 19, -1, 23, 14, 13, 18, -1, 24, 15, 14, 23, -1, 25, 16, 15, 24, -1, 22, 17, 16, 25, -1, 9, 4, 17, 22, -1, 12, 3, 4, 9, -1, 11, 2, 3, 12, -1, 10, 1, 2, 11, -1, 5, 0, 1, 10, -1, 18, 13, 0, 5]) + tgt_cI = DataArrayInt([0, 90]) + self.assertEqual(c.getValues(), tgt_c.getValues()) + self.assertEqual(cI.getValues(), tgt_cI.getValues()) + pass + pass if __name__ == '__main__':