From: Anthony Geay Date: Thu, 6 Apr 2017 12:20:32 +0000 (+0200) Subject: On the road X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=a10e758196a1039f418c9867adb056f2d1543b21;p=tools%2Fmedcoupling.git On the road --- diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 2757c52ae..87c0118f0 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -6391,8 +6391,15 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const for(int i=0;i=3) + if(nodalI[1]-nodalI[0]>=4) { + double aa[3]={coor[nodal[nodalI[0]+1+1]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0], + coor[nodal[nodalI[0]+1+1]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1], + coor[nodal[nodalI[0]+1+1]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]} + ,bb[3]={coor[nodal[nodalI[0]+1+2]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0], + coor[nodal[nodalI[0]+1+2]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1], + coor[nodal[nodalI[0]+1+2]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]}; + double cc[3]={aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]}; for(int j=0;j<3;j++) { int nodeId(nodal[nodalI[0]+1+j]); @@ -6404,14 +6411,34 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const throw INTERP_KERNEL::Exception(oss.str()); } } + if(sqrt(cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2])>1e-7) + { + INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); + retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; + } + else + { + if(nodalI[1]-nodalI[0]==4) + { + std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : cell" << i << " : Presence of The 3 colinear points !"; + throw INTERP_KERNEL::Exception(oss.str()); + } + // + double dd[3]={0.,0.,0.}; + for(int offset=nodalI[0]+1;offset()); + int nbOfNodesInCell(nodalI[1]-nodalI[0]-1); + std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies(),1./(double)nbOfNodesInCell)); + std::copy(dd,dd+3,matrix+4*2); + INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); + retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; + } } else { std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !"; throw INTERP_KERNEL::Exception(oss.str()); } - INTERP_KERNEL::inverseMatrix(matrix,4,matrix2); - retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15]; } return ret.retn(); } diff --git a/src/MEDCoupling/MEDCouplingVoronoi.cxx b/src/MEDCoupling/MEDCouplingVoronoi.cxx index b7412321f..40d6b1215 100644 --- a/src/MEDCoupling/MEDCouplingVoronoi.cxx +++ b/src/MEDCoupling/MEDCouplingVoronoi.cxx @@ -168,9 +168,11 @@ MCAuto MergeVorCells3D(const std::vector< MCAuto dummy3(p->mergeNodes(eps,dummy,dummy2)); } static int iii=0; - std::ostringstream oss; - oss << "merger_" << iii++ << ".vtu"; - p->writeVTK(oss.str()); + { + std::ostringstream oss; + oss << "merger_" << iii << ".vtu"; + p->writeVTK(oss.str()); + } MCAuto edgeToKeep; MCAuto p0; { @@ -218,6 +220,12 @@ MCAuto MergeVorCells3D(const std::vector< MCAutoinsertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]); + + { + std::ostringstream oss; + oss << "merger_resu_" << iii++ << ".vtu"; + ret->writeVTK(oss.str()); + } return ret; } @@ -480,15 +488,9 @@ MCAuto MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh * vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn); if(polygsToIterOn.size()<1) throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !"); - std::set elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone; - std::size_t ii(0); std::vector< MCAuto > newVorCells; - MCAuto d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New()); - MCAuto faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI)); - // - while(!elemsToDo.empty()) + for(int poly=0;polygetNumberOfCells();poly++) { - int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly); const double *seed(pts+3*poly); MCAuto tile(l0[poly]); tile->zipCoords(); @@ -506,49 +508,8 @@ MCAuto MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh * newVorCell->zipCoords(); MCAuto modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true)); modifiedCell->zipCoords(); - l0[poly]=modifiedCell; - if(std::find(polygsToIterOn.begin(),polygsToIterOn.end(),poly)!=polygsToIterOn.end())// we iterate on a polyhedron containg the point to add pt -> add cells sharing faces with just computed newVorCell - { -#if 0 - MCAuto faces2; - { - MCAuto d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New()); - faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I); - } - MCAuto faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face - MCAuto facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true)); - // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates - MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol); - MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3); - INTERP_KERNEL::Interpolation3DSurf interpolation; - interpolation.setMinDotBtwPlane3DSurfIntersect(eps2); - interpolation.setMaxDistance3DSurfIntersect(eps); - interpolation.setPrecision(1e-12); - std::vector > matrix; - interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0"); - std::set zeCandidates; - for(std::vector >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++) - for(std::map::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++) - { - int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0)); - for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++) - { - if(*it4!=poly) - zeCandidates.insert(*it4); - } - } -#endif - std::set zeCandidates; - for(int j=0;jgetNumberOfCells();j++) - zeCandidates.insert(j); - std::set tmp,newElementsToDo; - std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin())); - std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin())); - elemsToDo=newElementsToDo; - } - // newVorCells.push_back(newVorCell); - ii++; + l0[poly]=modifiedCell; } l0.push_back(MergeVorCells3D(newVorCells,eps)); } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py index dd9a2e8b9..d226c7acf 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest5.py @@ -4491,6 +4491,15 @@ class MEDCouplingBasicsTest5(unittest.TestCase): self.assertEqual(f2.getName(),"field") self.assertEqual(f2.getTime(),[1.,2,3]) pass + + def testBugInComputationOfEqOfPlane1(self): + coo=DataArrayDouble([-1.0, 1.0, -0.3872983455657959, -1.0, 1.0, 0.3872983455657959, -1.0, 1.0, 0.693649172782898, 1.0, 1.0, 0.693649172782898, 1.0, 1.0, 0.3872983455657959, 1.0, 1.0, -0.3872983455657959],6,3) + m=MEDCouplingUMesh("",2) + m.setCoords(coo) + m.allocateCells() + 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 pass