Salome HOME
Little refactory abn/fix_intersec
authorgeay <anthony.geay@cea.fr>
Thu, 20 Feb 2014 11:38:38 +0000 (12:38 +0100)
committergeay <anthony.geay@cea.fr>
Thu, 20 Feb 2014 11:38:38 +0000 (12:38 +0100)
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx
src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling/MEDCouplingUMesh.hxx
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 431ee8ad93e89862b115081ee6ddf34abf658809..9bbe6ce5527972076ed35979623475828e6196d2 100644 (file)
@@ -70,41 +70,63 @@ QuadraticPolygon::~QuadraticPolygon()
 {
 }
 
-QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes, bool isClosed)
+QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
 {
-  QuadraticPolygon *ret=new QuadraticPolygon;
+  QuadraticPolygon *ret(new QuadraticPolygon);
   std::size_t size=nodes.size();
-  for(std::size_t i=0;i< (size - (isClosed?0:1));i++)
+  for(std::size_t i=0;i<size;i++)
     {
       ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
       nodes[i]->decrRef();
     }
-  if(!isClosed)
-    nodes[size-1]->decrRef();
   return ret;
 }
 
-QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes, bool isClosed)
+QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
 {
-  QuadraticPolygon *ret=new QuadraticPolygon;
+  QuadraticPolygon *ret(new QuadraticPolygon);
   std::size_t size=nodes.size();
-  std::size_t quad_offset = isClosed ? (size/2) : (size/2+1);
-  for(std::size_t i = 0; i < size/2; i++)
+  for(std::size_t i=0;i<size/2;i++)
+
     {
       EdgeLin *e1,*e2;
-      e1=new EdgeLin(nodes[i],nodes[i+quad_offset]);
-      e2=new EdgeLin(nodes[i+quad_offset],nodes[(i+1)%quad_offset]);
+      e1=new EdgeLin(nodes[i],nodes[i+size/2]);
+      e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]);
       SegSegIntersector inters(*e1,*e2);
       bool colinearity=inters.areColinears();
       delete e1; delete e2;
       if(colinearity)
-        ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%quad_offset]));
+        ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)]));
       else
-        ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+quad_offset],nodes[(i+1)%quad_offset]));
-      nodes[i]->decrRef(); nodes[i+quad_offset]->decrRef();
+        ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)]));
+      nodes[i]->decrRef(); nodes[i+size/2]->decrRef();
     }
-  if(!isClosed)
-    nodes[quad_offset-1]->decrRef();
+  return ret;
+}
+
+Edge *QuadraticPolygon::BuildLinearEdge(std::vector<Node *>& nodes)
+{
+  if(nodes.size()!=2)
+    throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildLinearEdge : input vector is expected to be of size 2 !");
+  Edge *ret(new EdgeLin(nodes[0],nodes[1]));
+  nodes[0]->decrRef(); nodes[1]->decrRef();
+  return ret;
+}
+
+Edge *QuadraticPolygon::BuildArcCircleEdge(std::vector<Node *>& nodes)
+{
+  if(nodes.size()!=3)
+    throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildArcCircleEdge : input vector is expected to be of size 3 !");
+  EdgeLin *e1(new EdgeLin(nodes[0],nodes[2])),*e2(new EdgeLin(nodes[2],nodes[1]));
+  SegSegIntersector inters(*e1,*e2);
+  bool colinearity=inters.areColinears();
+  delete e1; delete e2;
+  Edge *ret(0);
+  if(colinearity)
+    ret=new EdgeLin(nodes[0],nodes[1]);
+  else
+    ret=new EdgeArcCircle(nodes[0],nodes[2],nodes[1]);
+  nodes[0]->decrRef(); nodes[1]->decrRef(); nodes[2]->decrRef();
   return ret;
 }
 
index a8097a2566829fba2c80ce8dcc1a1b620afd76b5..c6d235e1cd6ed8bb10c82d6d7a1c8641ab33282f 100644 (file)
@@ -46,8 +46,10 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT QuadraticPolygon() { }
     INTERPKERNEL_EXPORT QuadraticPolygon(const QuadraticPolygon& other):ComposedEdge(other) { }
     INTERPKERNEL_EXPORT QuadraticPolygon(const char *fileName);
-    INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector<Node *>& nodes, bool isClosed=true);
-    INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector<Node *>& nodes, bool isClosed=true);
+    INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector<Node *>& nodes);
+    INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector<Node *>& nodes);
+    INTERPKERNEL_EXPORT static Edge *BuildLinearEdge(std::vector<Node *>& nodes);
+    INTERPKERNEL_EXPORT static Edge *BuildArcCircleEdge(std::vector<Node *>& nodes);
     INTERPKERNEL_EXPORT static void BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName);
     INTERPKERNEL_EXPORT ~QuadraticPolygon();
     INTERPKERNEL_EXPORT void closeMe() const;
index 7f339ddeae33d9ecbdad116983e21ab47abe290f..05bd3c7a7dc6a96a65c57ebeb2d9f27bf4d9b0b0 100644 (file)
@@ -6272,9 +6272,9 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const
 DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) const
 {
   int mDim(getMeshDimension()),sDim(getSpaceDimension());
-  if(sDim!=2)  // Compute refined boundary box for quadratic elements only in 2D.
+  if((mDim==3 && sDim==3) || (mDim==2 && sDim==3) || (mDim==1 && sDim==1) || ( mDim==1 && sDim==3))  // Compute refined boundary box for quadratic elements only in 2D.
     return getBoundingBoxForBBTreeFast();
-  else
+  if((mDim==2 && sDim==2) || (mDim==1 && sDim==2))
     {
       bool presenceOfQuadratic(false);
       for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=_types.begin();it!=_types.end();it++)
@@ -6283,11 +6283,14 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) con
           if(cm.isQuadratic())
             presenceOfQuadratic=true;
         }
-      if(presenceOfQuadratic)
+      if(!presenceOfQuadratic)
+        return getBoundingBoxForBBTreeFast();
+      if(mDim==2 && sDim==2)
         return getBoundingBoxForBBTree2DQuadratic(arcDetEps);
       else
-        return getBoundingBoxForBBTreeFast();
+        return getBoundingBoxForBBTree1DQuadratic(arcDetEps);
     }
+  throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree : Managed dimensions are (mDim=1,sDim=1), (mDim=1,sDim=2), (mDim=1,sDim=3), (mDim=2,sDim=2), (mDim=2,sDim=3) and (mDim=3,sDim=3) !");
 }
 
 /*!
@@ -6349,12 +6352,13 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const
  * \throw If \a this is not fully defined.
  * \throw If \a this is not a mesh with meshDimension equal to 2.
  * \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic
  */
 DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
 {
   checkFullyDefined();
   int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
-  if(spaceDim!=2)
+  if(spaceDim!=2 || spaceDim!=2)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!");
   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
   double *bbox(ret->getPointer());
@@ -6364,7 +6368,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc
     {
       const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
       int sz(connI[1]-connI[0]-1);
-      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1e-12;
+      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
       std::vector<INTERP_KERNEL::Node *> nodes(sz);
       INTERP_KERNEL::QuadraticPolygon *pol(0);
       for(int j=0;j<sz;j++)
@@ -6373,15 +6377,60 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc
           nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]);
         }
       if(!cm.isQuadratic())
-        pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes, mDim==2);
+        pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
       else
-        pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes, mDim==2);
+        pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
       INTERP_KERNEL::Bounds b; pol->fillBounds(b); delete pol;
       bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); 
     }
   return ret.retn();
 }
 
+/*!
+ * This method aggregates the bbox of each 1D cell in \a this considering the whole shape. This method is particularly
+ * useful for 2D meshes having quadratic cells
+ * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers
+ * the two extremities of the arc of circle).
+ * 
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
+ * \throw If \a this is not fully defined.
+ * \throw If \a this is not a mesh with meshDimension equal to 1.
+ * \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic
+ */
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arcDetEps) const
+{
+  checkFullyDefined();
+  int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
+  if(spaceDim!=2 || spaceDim!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic : This method should be applied on mesh with mesh dimension equal to 1 and space dimension also equal to 2!");
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
+  double *bbox(ret->getPointer());
+  const double *coords(_coords->getConstPointer());
+  const int *conn(_nodal_connec->getConstPointer()),*connI(_nodal_connec_index->getConstPointer());
+  for(int i=0;i<nbOfCells;i++,bbox+=4,connI++)
+    {
+      const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
+      int sz(connI[1]-connI[0]-1);
+      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
+      std::vector<INTERP_KERNEL::Node *> nodes(sz);
+      INTERP_KERNEL::Edge *edge(0);
+      for(int j=0;j<sz;j++)
+        {
+          int nodeId(conn[*connI+1+j]);
+          nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]);
+        }
+      if(!cm.isQuadratic())
+        edge=INTERP_KERNEL::QuadraticPolygon::BuildLinearEdge(nodes);
+      else
+        edge=INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(nodes);
+      const INTERP_KERNEL::Bounds& b(edge->getBounds());
+      bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); edge->decrRef();
+    }
+  return ret.retn();
+}
+
 /// @cond INTERNAL
 
 namespace ParaMEDMEMImpl
index eaacda6841524fc7ea2a053683f54b8accd10736..232784dd14d90377537873607a00e0f22f8c8b19 100644 (file)
@@ -167,6 +167,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
     MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTreeFast() const;
     MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const;
+    MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const;
     MEDCOUPLING_EXPORT MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy);
     MEDCOUPLING_EXPORT bool isFullyQuadratic() const;
     MEDCOUPLING_EXPORT bool isPresenceOfQuadratic() const;
index 8723563df575f94a4290808acd7527851e8293c7..d3dce4efd53f1532339fbffda3b47d61302649d8 100644 (file)
@@ -270,6 +270,7 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf3DMesh;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTreeFast;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__;
 %newobject ParaMEDMEM::MEDCoupling1GTUMesh::New;
@@ -1540,6 +1541,7 @@ namespace ParaMEDMEM
     DataArrayInt *buildUnionOf3DMesh() const throw(INTERP_KERNEL::Exception);
     DataArrayDouble *getBoundingBoxForBBTreeFast() const throw(INTERP_KERNEL::Exception);
     DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
+    DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);