Intersect2DMeshWith1DLine: starting work to output refined tool mesh

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index e163c9a4070e4e5ef165704e44e33baafc29ef23..453a1a1b4699fa64e73db0d8370af809246fd9af 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -8743,7 +8743,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1
  *  \throw If m2 is not a (piecewise) line (i.e. if a point has more than 2 adjacent segments)
  */
 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *mesh2,
-                                                              double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2)
+                                                              double eps, MEDCouplingUMesh *& mesh2Refined,DataArrayInt *&cellNb1,
+                                                              DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2)
 {
   m1->checkFullyDefined();
   mesh2->checkFullyDefined();
@@ -8776,7 +8777,11 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingU
   BuildIntersectEdges(m1Desc,m2,addCoo,subDiv2,intersectEdge2);
   subDiv2.clear(); dd5=0;
 
-  // Step 3: this is where we have a significant difference with Intersect2DMeshes
+  // Step 3: build a finer tool mesh with all extra intersection points added. It will be
+  // part of the result
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m2fine = BuildRefinedMeshFromCut(m2, addCoo, intersectEdge2, renumb);
+
+  // Step 4: this is where we have a significant difference with Intersect2DMeshes
   std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
   std::vector<int> cNb1,cNb2,cNbI2; //no DataArrayInt because interface with Geometric2D
   cNbI2.push_back(0);
@@ -8801,6 +8806,9 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingU
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI2=DataArrayInt::New(); cI2->alloc((int)cNbI2.size(),1); std::copy(cNbI2.begin(),cNbI2.end(),cI2->getPointer());
 
+  // Renumber connectivity of refined tool:
+  m2fine->renumberCells(renumb->getConstPointer(), true);
+
   const int * ptr = renumb->getConstPointer();
   c2->transformWithIndArr(ptr, ptr+nc2);
   ret->setConnectivity(conn,connI,true);
@@ -8809,6 +8817,68 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingU
   return ret.retn();
 }
 
+MEDCouplingUMesh * MEDCouplingUMesh::BuildRefinedMeshFromCut(const MEDCouplingUMesh *m2, std::vector<double>& addCoo,
+                                                             std::vector<std::vector<int> >& intesctEdges2)
+{
+  const int *conn2 = m2->getNodalConnectivity()->getConstPointer();
+  const int *connI2 = m2->getNodalConnectivityIndex()->getConstPointer();
+  const double *coo2 = m2->getCoords()->getConstPointer();
+  std::vector<double> addCoordsQuadratic;
+
+  std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s;
+  std::map<INTERP_KERNEL::Node *,int> mapp;
+  std::map<int,INTERP_KERNEL::Node *> mappRev;
+  int cntEdges = 0;
+  for (int i = 0; i < m2->getNumberOfCells();i++)
+    {
+      INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[i]];
+      const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
+      for(int j = 0; j < intesctEdges2[i].size()/2; j++,cntEdges++)
+        {
+          pol2s.push_back(INTERP_KERNEL::QuadraticPolygon());
+          //INTERP_KERNEL::appendEdgeFromCrudeDataArray<ALL_C_MODE>(i,0,2,mapp,cm2.isQuadratic(),nodalBg,coords,intersectEdges);
+//      pol2s.back().buildFromCrudeDataArrayOneSeg(mappRev,cm2.isQuadratic(),conn2+connI2[i]+1,coo2,i,intesctEdges2,
+//          pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
+        }
+    }
+
+  // Explode QP into c,cI format:
+  DataArrayInt *conn = DataArrayInt::New(), *connI = DataArrayInt::New();
+  for (std::vector<INTERP_KERNEL::QuadraticPolygon>::const_iterator it = pol2s.begin(); it != pol2s.end(); it++)
+    (*it)->appendCrudeData(mapp,0,addCoordsQuadratic,conn,connI);
+
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m2fine = MEDCouplingUMesh::New("m2_refined", 2);
+  m2fine->setConnectivity(conn, connI, true);
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> qDa=DataArrayDouble::New();
+  qDa->alloc((int)(addCoordsQuadratic.size())/2,2);
+  std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),qDa->getPointer());
+  std::vector<const DataArrayDouble *> coordss(2);
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa=DataArrayDouble::New();
+  addCooDa->alloc((int)(addCoo.size())/2,2);
+  std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
+  coordss[0] = m2->getCoords();
+  coordss[1] = addCooDa;
+  coordss[2] = qDa;
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo = DataArrayDouble::Aggregate(coordss);
+  m2fine->setCoords(coo);
+
+  // Now rewrite intesctEdges2:
+  std::vector< std::vector <int> > intesctSave(intesctEdges2);
+  intesctEdges2.resize(cntEdges);
+
+  for (int i = 0, l = 0; i < m2->getNumberOfCells(); i++)
+    {
+      for (int j = 0; j < intesctSave[i].size()/2; j++, l++)
+        {
+          intesctEdges2[l].resize(2);
+          intesctEdges2[l][0] = intesctSave[i][j*2];
+          intesctEdges2[l][1] = intesctSave[i][j*2+1];
+        }
+    }
+
+  return m2fine.retn();
+}
+
 /**
  * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
  * (newly created) nodes corresponding to the edge intersections.
@@ -8902,7 +8972,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo
 }
 
 void MEDCouplingUMesh::Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
-                                                  const MEDCouplingUMesh *m2, const std::vector<std::vector<int> >& intesctEdges2,
+                                                  const MEDCouplingUMesh *m2split, const std::vector<std::vector<int> >& intesctEdges2,
                                                   const std::vector<double>& addCoords,
                                                   std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI,
                                                   std::vector<int>& cNb1, std::vector<int>& cNb2, std::vector<int>& cNbI2)
@@ -8912,15 +8982,15 @@ void MEDCouplingUMesh::Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh
   const int *conn1=m1->getNodalConnectivity()->getConstPointer();
   const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer();
   int offset1=m1->getNumberOfNodes();
-  const double *coo2=m2->getCoords()->getConstPointer();
-  const int *conn2=m2->getNodalConnectivity()->getConstPointer();
-  const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer();
-  int offset2=offset1+m2->getNumberOfNodes();
+  const double *coo2=m2split->getCoords()->getConstPointer();
+  const int *conn2=m2split->getNodalConnectivity()->getConstPointer();
+  const int *connI2=m2split->getNodalConnectivityIndex()->getConstPointer();
+  int offset2=offset1+m2split->getNumberOfNodes();
   int offset3=offset2+((int)addCoords.size())/2;
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2split->getBoundingBoxForBBTree());
   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
   // Here a BBTree on 1D-segments from the tool mesh, not on segments:
-  BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
+  BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2split->getNumberOfCells(),eps);
   int ncell1=m1->getNumberOfCells();
   crI.push_back(0);
   // for all 2D cells in base mesh, find intersecting segments

diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 6cc91f35744680d8c641e24bc72ee3b00bf89bf1..669fe090dcede7bfa9fe15c28ee8774850a38a01 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -241,7 +241,8 @@ namespace ParaMEDMEM
     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);
     MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
-                                                          double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2);
+                                                          double eps, MEDCouplingUMesh *& mesh2Refined,
+                                                          DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2);
     MEDCOUPLING_EXPORT static bool BuildConvexEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, DataArrayInt *nodalConnecOut);
     MEDCOUPLING_EXPORT static bool RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, const int *idsToRemoveEnd, DataArrayInt *arr, DataArrayInt *arrIndx, int offsetForRemoval=0);
     MEDCOUPLING_EXPORT static void ExtractFromIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
@@ -328,6 +329,8 @@ namespace ParaMEDMEM
                                                       const std::vector<double>& addCoords,
                                                       std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI,
                                                       std::vector<int>& cNb1, std::vector<int>& cNb2, std::vector<int>& cNbI2);
+    static MEDCouplingUMesh * BuildRefinedMeshFromCut(const MEDCouplingUMesh *m2, std::vector<double>& addCoo,
+                                                      std::vector<std::vector<int> >& intesctEdges2);
     static void AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
                                               const int *nodal3DCurve, const int *nodalIndx3DCurve,
                                               const int *desc, const int *descIndx, std::vector< std::pair<int,int> >& cut3DSurf);

diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i
index b9c814669d9c1a99651769a28ac1ca8dc3b74ed0..42a79aa74b02c91dfc91ee9f8dc4918fc6b39c56 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingCommon.i
+++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i
@@ -2386,11 +2386,13 @@ namespace ParaMEDMEM
       static PyObject *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps) throw(INTERP_KERNEL::Exception)
       {
         DataArrayInt *cellNb1=0,*cellNb2=0;
-        MEDCouplingUMesh *mret=MEDCouplingUMesh::Intersect2DMeshes(m1,m2,eps,cellNb1,cellNb2);
+        MEDCouplingUMesh *m2fine=0;
+        MEDCouplingUMesh *mret=MEDCouplingUMesh::Intersect2DMeshes(m1,m2,eps,m2fine,cellNb1,cellNb2);
         PyObject *ret=PyTuple_New(3);
         PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(mret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 ));
-        PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb1),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
-        PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+        PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(mret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 ));
+        PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb1),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+        PyTuple_SetItem(ret,3,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
         return ret;
       }