On the road

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 2757c52aec249dd9cab5e70cea4c85b534480186..87c0118f02c07dc438dc139596d2e00a4fe12faf 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -6391,8 +6391,15 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
   for(int i=0;i<nbOfCells;i++,nodalI++,retPtr+=4)
     {
       double matrix[16]={0,0,0,1,0,0,0,1,0,0,0,1,1,1,1,0},matrix2[16];
-      if(nodalI[1]-nodalI[0]>=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<nodalI[1];offset++)
+                std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
+              int nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
+              std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),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 b7412321fef96cf535a89de14f19cd6568c1332f..40d6b1215c7f400ab7246599ee94c7fdad5aa23e 100644 (file)
--- a/src/MEDCoupling/MEDCouplingVoronoi.cxx
+++ b/src/MEDCoupling/MEDCouplingVoronoi.cxx
@@ -168,9 +168,11 @@ MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUM
     MCAuto<DataArrayInt> 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<DataArrayInt> edgeToKeep;
   MCAuto<MEDCouplingUMesh> p0;
   {
@@ -218,6 +220,12 @@ MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUM
       }
   }
   ret->insertNextCell(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<MEDCouplingUMesh> 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<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
-      std::size_t ii(0);
       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
-      MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
-      MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
-      //
-      while(!elemsToDo.empty())
+      for(int poly=0;poly<vorTess->getNumberOfCells();poly++)
         {
-          int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
           const double *seed(pts+3*poly);
           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
           tile->zipCoords();
@@ -506,49 +508,8 @@ MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *
           newVorCell->zipCoords();
           MCAuto<MEDCouplingUMesh> 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<MEDCouplingUMesh> faces2;
-              {
-                MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
-                faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
-              }
-              MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
-              MCAuto<MEDCouplingUMesh> 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<std::map<int,double> > matrix;
-              interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
-              std::set<int> zeCandidates;
-              for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
-                for(std::map<int,double>::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<int> zeCandidates;
-              for(int j=0;j<vorTess->getNumberOfCells();j++)
-                zeCandidates.insert(j);
-              std::set<int> 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 dd9a2e8b9d486ef6f65d9ffdd973aed3f9b1bb7a..d226c7acf68f7bda4dbcc38b50a488cfff56d508 100644 (file)
--- 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