]> OCCT Git - occt.git/commitdiff
0028089: Mesh - New algorithm for triangulation of 2d polygons
authoroan <oan@opencascade.com>
Wed, 10 Jul 2019 10:04:25 +0000 (13:04 +0300)
committerbugmaster <bugmaster@opencascade.com>
Sun, 27 Sep 2020 09:00:07 +0000 (12:00 +0300)
Added custom meshing core algorithm to generate base mesh using Delabella library,
which can be enabled via IMeshTools_Parameters::MeshAlgo option or CSF_MeshAlgo environment variable.

Do not fill cirles filter upon explicit initialization.
Call base postProcessMesh functionality after initialization of circles in BRepMesh_CustomDelaunayBaseMeshAlgo.

Added Vsprintf() wrapper for vsprintf() preserving C locale.

30 files changed:
dox/overview/overview.md
dox/user_guides/modeling_algos/modeling_algos.md
src/BRepMesh/BRepMesh_BaseMeshAlgo.hxx
src/BRepMesh/BRepMesh_ConstrainedBaseMeshAlgo.hxx
src/BRepMesh/BRepMesh_Context.cxx
src/BRepMesh/BRepMesh_Context.hxx
src/BRepMesh/BRepMesh_CustomBaseMeshAlgo.hxx
src/BRepMesh/BRepMesh_CustomDelaunayBaseMeshAlgo.hxx
src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx [new file with mode: 0644]
src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx [new file with mode: 0644]
src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx [new file with mode: 0644]
src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx [new file with mode: 0644]
src/BRepMesh/BRepMesh_Delaun.cxx
src/BRepMesh/BRepMesh_Delaun.hxx
src/BRepMesh/BRepMesh_DelaunayBaseMeshAlgo.hxx
src/BRepMesh/BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx
src/BRepMesh/BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx
src/BRepMesh/BRepMesh_IncrementalMesh.cxx
src/BRepMesh/BRepMesh_MeshAlgoFactory.hxx
src/BRepMesh/FILES
src/BRepMesh/delabella.cpp [new file with mode: 0644]
src/BRepMesh/delabella.pxx [new file with mode: 0644]
src/IMeshTools/FILES
src/IMeshTools/IMeshTools_MeshAlgoType.hxx [new file with mode: 0644]
src/IMeshTools/IMeshTools_Parameters.hxx
src/MeshTest/MeshTest.cxx
src/Standard/Standard_CString.cxx
src/Standard/Standard_CString.hxx
tests/bugs/mesh/bug28089_1 [new file with mode: 0644]
tests/bugs/mesh/bug28089_2 [new file with mode: 0644]

index 6eb6c1b2b9542eeae170bfa86f6c2239d951384e..24d666b6561b10b2a7fb474a165e254b9ad81971 100644 (file)
@@ -156,6 +156,10 @@ RapidJSON is optionally used by OCCT for reading glTF files (https://rapidjson.o
 **DejaVu** fonts are a font family based on the Vera Fonts under a permissive license (MIT-like, https://dejavu-fonts.github.io/License.html).
 DejaVu Sans (basic Latin sub-set) is used by OCCT as fallback font when no system font is available.
 
+**Delabella** is an open-source, cross-platform implementation of the Newton Apple Wrapper algorithm producing 2D Delaunay triangulation.
+Delabella is used by BRepMesh as one of alternative 2D triangulation algorithms.
+Delabella is licensed under the MIT license (https://github.com/msokalski/delabella).
+
 Adobe Systems, Inc. provides **Adobe Reader**, which can be used to view files in Portable Document Format (PDF). 
 
 @section OCCT_OVW_SECTION_3 Documentation
index df0cdcd104ab189ad535a45cdef95c43acb741a0..9977a562e8c4593f71294829287810b3c3cd795b 100644 (file)
@@ -3251,6 +3251,40 @@ In general, face meshing algorithms have the following structure:
   * Classes *BRepMesh_DelaunayNodeInsertionMeshAlgo* and *BRepMesh_SweepLineNodeInsertionMeshAlgo* implement algorithm-specific functionality related to addition of internal nodes supplementing functionality provided by *BRepMesh_NodeInsertionMeshAlgo*;
   * *BRepMesh_DelaunayDeflectionControlMeshAlgo* extends functionality of *BRepMesh_DelaunayNodeInsertionMeshAlgo* by additional procedure controlling deflection of generated triangles.
 
+BRepMesh provides user a way to switch default triangulation algorithm to a custom one, either implemented by user or available worldwide.
+There are three base classes that can be currently used to integrate 3rd-party algorithms:
+  * *BRepMesh_ConstrainedBaseMeshAlgo* base class for tools providing generation of triangulations with constraints requiring no common processing by BRepMesh;
+  * *BRepMesh_CustomBaseMeshAlgo* provides the entry point for generic algorithms without support of constraints and supposed for generation of base mesh only.
+    Constraint edges are processed using standard functionality provided by the component itself upon background mesh produced by 3rd-party solver;
+  * *BRepMesh_CustomDelaunayBaseMeshAlgo* contains initialization part for tools used by BRepMesh for checks or optimizations using results of 3rd-party algorithm.
+
+Meshing algorithms could be provided by implemeting IMeshTools_MeshAlgoFactory with related interfaces and passing it to BRepMesh_Context::SetFaceDiscret().
+OCCT comes with two base 2D meshing algorithms: BRepMesh_MeshAlgoFactory (used by default) and BRepMesh_DelabellaMeshAlgoFactory.
+
+The following example demonstrates how it could be done from Draw environment:
+~~~~~
+psphere s 10
+
+### Default Algo ###
+incmesh s 0.0001 -algo default
+
+### Delabella Algo ###
+incmesh s 0.0001 -algo delabella
+~~~~~
+
+The code snippet below shows passing a custom mesh factory to BRepMesh_IncrementalMesh:
+~~~~~
+IMeshTools_Parameters aMeshParams;
+Handle(IMeshTools_Context) aContext = new BRepMesh_Context();
+aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_DelabellaMeshAlgoFactory()));
+
+BRepMesh_IncrementalMesh aMesher;
+aMesher.SetShape (aShape);
+aMesher.ChangeParameters() = aMeshParams;
+
+aMesher.Perform (aContext);
+~~~~~
+
 #### Range splitter
 Range splitter tools provide functionality to generate internal surface nodes defined within the range computed using discrete model data. The base functionality is provided by *BRepMesh_DefaultRangeSplitter* which can be used without modifications in case of planar surface. The default splitter does not generate any internal node.
 
index 1ab210fe8e5376d77e2325acbc73d8c98597f2fc..ddf4419933797c105ccf1ec1850bf66e56ee9868 100644 (file)
@@ -25,7 +25,7 @@
 class BRepMesh_DataStructureOfDelaun;
 class BRepMesh_Delaun;
 
-//! Class provides base fuctionality for algorithms building face triangulation.
+//! Class provides base functionality for algorithms building face triangulation.
 //! Performs initialization of BRepMesh_DataStructureOfDelaun and nodes map structures.
 class BRepMesh_BaseMeshAlgo : public IMeshTools_MeshAlgo
 {
@@ -108,7 +108,7 @@ protected:
 
 private:
 
-  //! If the given edge has another pcurve for current face coinsiding with specified one,
+  //! If the given edge has another pcurve for current face coinciding with specified one,
   //! returns TopAbs_INTERNAL flag. Elsewhere returns orientation of specified pcurve.
   TopAbs_Orientation fixSeamEdgeOrientation(
     const IMeshData::IEdgeHandle&   theDEdge,
index 38a63aab73453388c77a878c323a9506f556b33e..59e983602778fbad0e8c6f57ad5b26a3aebbd286 100644 (file)
@@ -23,7 +23,7 @@
 class BRepMesh_DataStructureOfDelaun;
 class BRepMesh_Delaun;
 
-//! Class provides base fuctionality to build face triangulation using Dealunay approach.
+//! Class provides base functionality to build face triangulation using Dealunay approach.
 //! Performs generation of mesh using raw data from model.
 class BRepMesh_ConstrainedBaseMeshAlgo : public BRepMesh_BaseMeshAlgo
 {
@@ -49,7 +49,7 @@ protected:
     return std::pair<Standard_Integer, Standard_Integer> (-1, -1);
   }
 
-  //! Perfroms processing of generated mesh.
+  //! Performs processing of generated mesh.
   //! By default does nothing.
   //! Expected to be called from method generateMesh() in successor classes.
   virtual void postProcessMesh (BRepMesh_Delaun&              /*theMesher*/,
index e48d038467c3efea18e13a1fc2bd8f387f505dec..fe9b6520d70cdf8ffab7d867aa0cc279d30a46c2 100644 (file)
 #include <BRepMesh_FaceDiscret.hxx>
 #include <BRepMesh_ModelPreProcessor.hxx>
 #include <BRepMesh_ModelPostProcessor.hxx>
+
 #include <BRepMesh_MeshAlgoFactory.hxx>
+#include <BRepMesh_DelabellaMeshAlgoFactory.hxx>
+#include <Message.hxx>
+#include <OSD_Environment.hxx>
 
 //=======================================================================
 // Function: Constructor
 // Purpose : 
 //=======================================================================
-BRepMesh_Context::BRepMesh_Context ()
+BRepMesh_Context::BRepMesh_Context (IMeshTools_MeshAlgoType theMeshType)
 {
+  if (theMeshType == IMeshTools_MeshAlgoType_DEFAULT)
+  {
+    TCollection_AsciiString aValue = OSD_Environment ("CSF_MeshAlgo").Value();
+    aValue.LowerCase();
+    if (aValue == "watson"
+     || aValue == "0")
+    {
+      theMeshType = IMeshTools_MeshAlgoType_Watson;
+    }
+    else if (aValue == "delabella"
+          || aValue == "1")
+    {
+      theMeshType = IMeshTools_MeshAlgoType_Delabella;
+    }
+    else
+    {
+      if (!aValue.IsEmpty())
+      {
+        Message::SendWarning (TCollection_AsciiString("BRepMesh_Context, ignore unknown algorithm '") + aValue + "' specified in CSF_MeshAlgo variable");
+      }
+      theMeshType = IMeshTools_MeshAlgoType_Watson;
+    }
+  }
+
+  Handle (IMeshTools_MeshAlgoFactory) aAlgoFactory;
+  switch (theMeshType)
+  {
+    case IMeshTools_MeshAlgoType_DEFAULT:
+    case IMeshTools_MeshAlgoType_Watson:
+      aAlgoFactory = new BRepMesh_MeshAlgoFactory();
+      break;
+    case IMeshTools_MeshAlgoType_Delabella:
+      aAlgoFactory = new BRepMesh_DelabellaMeshAlgoFactory();
+      break;
+  }
+
   SetModelBuilder (new BRepMesh_ModelBuilder);
   SetEdgeDiscret  (new BRepMesh_EdgeDiscret);
   SetModelHealer  (new BRepMesh_ModelHealer);
   SetPreProcessor (new BRepMesh_ModelPreProcessor);
-  SetFaceDiscret  (new BRepMesh_FaceDiscret(new BRepMesh_MeshAlgoFactory));
+  SetFaceDiscret  (new BRepMesh_FaceDiscret (aAlgoFactory));
   SetPostProcessor(new BRepMesh_ModelPostProcessor);
 }
 
index a802a6b44bf05184561e903f8f801f9cabd277f3..e2b004d1d6521db77d1d2a9624838edd674e7eb3 100644 (file)
@@ -25,7 +25,7 @@ class BRepMesh_Context : public IMeshTools_Context
 public:
 
   //! Constructor.
-  Standard_EXPORT BRepMesh_Context ();
+  Standard_EXPORT BRepMesh_Context (IMeshTools_MeshAlgoType theMeshType = IMeshTools_MeshAlgoType_DEFAULT);
 
   //! Destructor.
   Standard_EXPORT virtual ~BRepMesh_Context ();
index e175091de90818829d2886a37ee01196e9942edb..74d6b0c81f4b12e7845d23684dfc76ca88f82df7 100644 (file)
@@ -25,7 +25,7 @@
 
 class BRepMesh_DataStructureOfDelaun;
 
-//! Class provides base fuctionality to build face triangulation using custom triangulation algorithm.
+//! Class provides base functionality to build face triangulation using custom triangulation algorithm.
 //! Performs generation of mesh using raw data from model.
 class BRepMesh_CustomBaseMeshAlgo : public BRepMesh_ConstrainedBaseMeshAlgo
 {
@@ -46,19 +46,42 @@ public:
 protected:
 
   //! Generates mesh for the contour stored in data structure.
-  Standard_EXPORT virtual void generateMesh () Standard_OVERRIDE
+  Standard_EXPORT virtual void generateMesh (const Message_ProgressRange& theRange) Standard_OVERRIDE
   {
     const Handle (BRepMesh_DataStructureOfDelaun)& aStructure = this->getStructure ();
+    const Standard_Integer aNodesNb = aStructure->NbNodes ();
+
     buildBaseTriangulation ();
 
     std::pair<Standard_Integer, Standard_Integer> aCellsCount = this->getCellsCount (aStructure->NbNodes ());
     BRepMesh_Delaun aMesher (aStructure, aCellsCount.first, aCellsCount.second, Standard_False);
+
+    const Standard_Integer aNewNodesNb = aStructure->NbNodes ();
+    const Standard_Boolean isRemoveAux = aNewNodesNb > aNodesNb;
+    if (isRemoveAux)
+    {
+      IMeshData::VectorOfInteger aAuxVertices (aNewNodesNb - aNodesNb);
+      for (Standard_Integer aExtNodesIt = aNodesNb + 1; aExtNodesIt <= aNewNodesNb; ++aExtNodesIt)
+      {
+        aAuxVertices.Append (aExtNodesIt);
+      }
+
+      // Set aux vertices if there are some to clean up mesh correctly.
+      aMesher.SetAuxVertices (aAuxVertices);
+    }
+
     aMesher.ProcessConstraints ();
 
+    // Destruction of triangles containing aux vertices added (possibly) during base mesh computation.
+    if (isRemoveAux)
+    {
+      aMesher.RemoveAuxElements ();
+    }
+
     BRepMesh_MeshTool aCleaner (aStructure);
     aCleaner.EraseFreeLinks ();
 
-    postProcessMesh (aMesher);
+    postProcessMesh (aMesher, theRange);
   }
 
 protected:
index 6bde7b90a8903c2c121d872068f752990f81eca9..2143db63591cf86e00e128964141124809f33610 100644 (file)
@@ -19,7 +19,7 @@
 class BRepMesh_DataStructureOfDelaun;
 class BRepMesh_Delaun;
 
-//! Class provides base fuctionality to build face triangulation using custom
+//! Class provides base functionality to build face triangulation using custom
 //! triangulation algorithm with possibility to modify final mesh.
 //! Performs generation of mesh using raw data from model.
 template<class BaseAlgo>
@@ -39,14 +39,15 @@ public:
 
 protected:
 
-  //! Perfroms processing of generated mesh.
-  virtual void postProcessMesh(BRepMesh_Delaun& theMesher)
+  //! Performs processing of generated mesh.
+  virtual void postProcessMesh (BRepMesh_Delaun& theMesher,
+                                const Message_ProgressRange& theRange)
   {
-    BaseAlgo::postProcessMesh (theMesher);
-
     const Handle(BRepMesh_DataStructureOfDelaun)& aStructure  = this->getStructure();
     std::pair<Standard_Integer, Standard_Integer> aCellsCount = this->getCellsCount (aStructure->NbNodes());
     theMesher.InitCirclesTool (aCellsCount.first, aCellsCount.second);
+
+    BaseAlgo::postProcessMesh (theMesher, theRange);
   }
 };
 
diff --git a/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.cxx
new file mode 100644 (file)
index 0000000..40bbb36
--- /dev/null
@@ -0,0 +1,193 @@
+// Created on: 2019-07-05
+// Copyright (c) 2019 OPEN CASCADE SAS
+// Created by: Oleg AGASHIN
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <BRepMesh_DelabellaBaseMeshAlgo.hxx>
+
+#include <BRepMesh_MeshTool.hxx>
+#include <BRepMesh_Delaun.hxx>
+#include <Message.hxx>
+
+#include <string.h>
+#include <stdarg.h>
+
+#include "delabella.pxx"
+
+IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_DelabellaBaseMeshAlgo, BRepMesh_CustomBaseMeshAlgo)
+
+namespace
+{
+  //! Redirect algorithm messages to OCCT messenger.
+  static int logDelabella2Occ (void* theStream, const char* theFormat, ...)
+  {
+    (void )theStream;
+    char aBuffer[1024]; // should be more than enough for Delabella messages
+
+    va_list anArgList;
+    va_start(anArgList, theFormat);
+    Vsprintf(aBuffer, theFormat, anArgList);
+    va_end(anArgList);
+
+    Message_Gravity aGravity = Message_Warning;
+    switch ((int )theFormat[1])
+    {
+      case int('E'): aGravity = Message_Fail;  break; // [ERR]
+      case int('W'): aGravity = Message_Trace; break; // [WRN]
+      case int('N'): aGravity = Message_Trace; break; // [NFO]
+    }
+    Message::Send (aBuffer, aGravity);
+    return 0;
+  }
+}
+
+//=======================================================================
+// Function: Constructor
+// Purpose :
+//=======================================================================
+BRepMesh_DelabellaBaseMeshAlgo::BRepMesh_DelabellaBaseMeshAlgo ()
+{
+}
+
+//=======================================================================
+// Function: Destructor
+// Purpose :
+//=======================================================================
+BRepMesh_DelabellaBaseMeshAlgo::~BRepMesh_DelabellaBaseMeshAlgo ()
+{
+}
+
+//=======================================================================
+//function : buildBaseTriangulation
+//purpose  :
+//=======================================================================
+void BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation()
+{
+  const Handle(BRepMesh_DataStructureOfDelaun)& aStructure = this->getStructure();
+
+  Bnd_B2d aBox;
+  const Standard_Integer aNodesNb = aStructure->NbNodes ();
+  std::vector<Standard_Real> aPoints (2 * (aNodesNb + 4));
+  for (Standard_Integer aNodeIt = 0; aNodeIt < aNodesNb; ++aNodeIt)
+  {
+    const BRepMesh_Vertex& aVertex = aStructure->GetNode (aNodeIt + 1);
+
+    const size_t aBaseIdx = 2 * static_cast<size_t> (aNodeIt);
+    aPoints[aBaseIdx + 0] = aVertex.Coord ().X ();
+    aPoints[aBaseIdx + 1] = aVertex.Coord ().Y ();
+
+    aBox.Add (gp_Pnt2d(aVertex.Coord ()));
+  }
+
+  aBox.Enlarge (0.1 * (aBox.CornerMax () - aBox.CornerMin ()).Modulus ());
+  const gp_XY aMin = aBox.CornerMin ();
+  const gp_XY aMax = aBox.CornerMax ();
+
+  aPoints[2 * aNodesNb + 0] = aMin.X ();
+  aPoints[2 * aNodesNb + 1] = aMin.Y ();
+  aStructure->AddNode (BRepMesh_Vertex (
+    aPoints[2 * aNodesNb + 0],
+    aPoints[2 * aNodesNb + 1], BRepMesh_Free));
+
+  aPoints[2 * aNodesNb + 2] = aMax.X ();
+  aPoints[2 * aNodesNb + 3] = aMin.Y ();
+  aStructure->AddNode (BRepMesh_Vertex (
+    aPoints[2 * aNodesNb + 2],
+    aPoints[2 * aNodesNb + 3], BRepMesh_Free));
+
+  aPoints[2 * aNodesNb + 4] = aMax.X ();
+  aPoints[2 * aNodesNb + 5] = aMax.Y ();
+  aStructure->AddNode (BRepMesh_Vertex (
+    aPoints[2 * aNodesNb + 4],
+    aPoints[2 * aNodesNb + 5], BRepMesh_Free));
+
+  aPoints[2 * aNodesNb + 6] = aMin.X ();
+  aPoints[2 * aNodesNb + 7] = aMax.Y ();
+  aStructure->AddNode (BRepMesh_Vertex (
+    aPoints[2 * aNodesNb + 6],
+    aPoints[2 * aNodesNb + 7], BRepMesh_Free));
+
+  const Standard_Real aDiffX = (aMax.X () - aMin.X ());
+  const Standard_Real aDiffY = (aMax.Y () - aMin.Y ());
+  for (size_t i = 0; i < aPoints.size(); i += 2)
+  {
+    aPoints[i + 0] = (aPoints[i + 0] - aMin.X ()) / aDiffX - 0.5;
+    aPoints[i + 1] = (aPoints[i + 1] - aMin.Y ()) / aDiffY - 0.5;
+  }
+
+  IDelaBella* aTriangulator = IDelaBella::Create();
+  if (aTriangulator == NULL) // should never happen
+  {
+    throw Standard_ProgramError ("BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation: unable creating a triangulation algorithm");
+  }
+
+  aTriangulator->SetErrLog (logDelabella2Occ, NULL);
+  try
+  {
+    const int aVerticesNb = aTriangulator->Triangulate (
+      static_cast<int>(aPoints.size () / 2),
+      &aPoints[0], &aPoints[1], 2 * sizeof (Standard_Real));
+
+    if (aVerticesNb > 0)
+    {
+      const DelaBella_Triangle* aTrianglePtr = aTriangulator->GetFirstDelaunayTriangle();
+      while (aTrianglePtr != NULL)
+      {
+        Standard_Integer aNodes[3] = {
+          aTrianglePtr->v[0]->i + 1,
+          aTrianglePtr->v[2]->i + 1,
+          aTrianglePtr->v[1]->i + 1
+        };
+
+        Standard_Integer aEdges       [3];
+        Standard_Boolean aOrientations[3];
+        for (Standard_Integer k = 0; k < 3; ++k)
+        {
+          const BRepMesh_Edge aLink (aNodes[k], aNodes[(k + 1) % 3], BRepMesh_Free);
+
+          const Standard_Integer aLinkInfo = aStructure->AddLink (aLink);
+          aEdges       [k] = Abs (aLinkInfo);
+          aOrientations[k] = aLinkInfo > 0;
+        }
+
+        const BRepMesh_Triangle aTriangle (aEdges, aOrientations, BRepMesh_Free);
+        aStructure->AddElement (aTriangle);
+
+        aTrianglePtr = aTrianglePtr->next;
+      }
+    }
+
+    aTriangulator->Destroy ();
+    aTriangulator = NULL;
+  }
+  catch (Standard_Failure const& theException)
+  {
+    if (aTriangulator != NULL)
+    {
+      aTriangulator->Destroy ();
+      aTriangulator = NULL;
+    }
+
+    throw Standard_Failure (theException);
+  }
+  catch (...)
+  {
+    if (aTriangulator != NULL)
+    {
+      aTriangulator->Destroy ();
+      aTriangulator = NULL;
+    }
+
+    throw Standard_Failure ("BRepMesh_DelabellaBaseMeshAlgo::buildBaseTriangulation: exception in triangulation algorithm");
+  }
+}
diff --git a/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx b/src/BRepMesh/BRepMesh_DelabellaBaseMeshAlgo.hxx
new file mode 100644 (file)
index 0000000..9d6f2e1
--- /dev/null
@@ -0,0 +1,46 @@
+// Created on: 2019-07-05
+// Copyright (c) 2019 OPEN CASCADE SAS
+// Created by: Oleg AGASHIN
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BRepMesh_DelabellaBaseMeshAlgo_HeaderFile
+#define _BRepMesh_DelabellaBaseMeshAlgo_HeaderFile
+
+#include <BRepMesh_CustomBaseMeshAlgo.hxx>
+#include <NCollection_Shared.hxx>
+#include <IMeshTools_Parameters.hxx>
+
+class BRepMesh_DataStructureOfDelaun;
+class BRepMesh_Delaun;
+
+//! Class provides base functionality to build face triangulation using Delabella project.
+//! Performs generation of mesh using raw data from model.
+class BRepMesh_DelabellaBaseMeshAlgo : public BRepMesh_CustomBaseMeshAlgo
+{
+public:
+
+  //! Constructor.
+  Standard_EXPORT BRepMesh_DelabellaBaseMeshAlgo ();
+
+  //! Destructor.
+  Standard_EXPORT virtual ~BRepMesh_DelabellaBaseMeshAlgo ();
+
+  DEFINE_STANDARD_RTTIEXT(BRepMesh_DelabellaBaseMeshAlgo, BRepMesh_CustomBaseMeshAlgo)
+
+protected:
+
+  //! Builds base triangulation using Delabella project.
+  Standard_EXPORT virtual void buildBaseTriangulation() Standard_OVERRIDE;
+};
+
+#endif
diff --git a/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.cxx
new file mode 100644 (file)
index 0000000..8b02c4a
--- /dev/null
@@ -0,0 +1,145 @@
+// Created on: 2019-07-05
+// Copyright (c) 2019 OPEN CASCADE SAS
+// Created by: Oleg AGASHIN
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <BRepMesh_DelabellaMeshAlgoFactory.hxx>
+#include <BRepMesh_DefaultRangeSplitter.hxx>
+#include <BRepMesh_NURBSRangeSplitter.hxx>
+#include <BRepMesh_SphereRangeSplitter.hxx>
+#include <BRepMesh_CylinderRangeSplitter.hxx>
+#include <BRepMesh_ConeRangeSplitter.hxx>
+#include <BRepMesh_TorusRangeSplitter.hxx>
+#include <BRepMesh_DelaunayBaseMeshAlgo.hxx>
+#include <BRepMesh_DelabellaBaseMeshAlgo.hxx>
+#include <BRepMesh_CustomDelaunayBaseMeshAlgo.hxx>
+#include <BRepMesh_DelaunayNodeInsertionMeshAlgo.hxx>
+#include <BRepMesh_DelaunayDeflectionControlMeshAlgo.hxx>
+#include <BRepMesh_BoundaryParamsRangeSplitter.hxx>
+
+namespace
+{
+  struct DefaultBaseMeshAlgo
+  {
+    typedef BRepMesh_DelaunayBaseMeshAlgo Type;
+  };
+
+  template<class RangeSplitter>
+  struct DefaultNodeInsertionMeshAlgo
+  {
+    typedef BRepMesh_DelaunayNodeInsertionMeshAlgo<RangeSplitter, BRepMesh_DelaunayBaseMeshAlgo> Type;
+  };
+
+  struct BaseMeshAlgo
+  {
+    typedef BRepMesh_DelabellaBaseMeshAlgo Type;
+  };
+
+  template<class RangeSplitter>
+  struct NodeInsertionMeshAlgo
+  {
+    typedef BRepMesh_DelaunayNodeInsertionMeshAlgo<RangeSplitter, BRepMesh_CustomDelaunayBaseMeshAlgo<BRepMesh_DelabellaBaseMeshAlgo> > Type;
+  };
+
+  template<class RangeSplitter>
+  struct DeflectionControlMeshAlgo
+  {
+    typedef BRepMesh_DelaunayDeflectionControlMeshAlgo<RangeSplitter, BRepMesh_CustomDelaunayBaseMeshAlgo<BRepMesh_DelabellaBaseMeshAlgo> > Type;
+  };
+}
+
+IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_DelabellaMeshAlgoFactory, IMeshTools_MeshAlgoFactory)
+
+//=======================================================================
+// Function: Constructor
+// Purpose :
+//=======================================================================
+BRepMesh_DelabellaMeshAlgoFactory::BRepMesh_DelabellaMeshAlgoFactory ()
+{
+}
+
+//=======================================================================
+// Function: Destructor
+// Purpose :
+//=======================================================================
+BRepMesh_DelabellaMeshAlgoFactory::~BRepMesh_DelabellaMeshAlgoFactory ()
+{
+}
+
+//=======================================================================
+// Function: GetAlgo
+// Purpose :
+//=======================================================================
+Handle(IMeshTools_MeshAlgo) BRepMesh_DelabellaMeshAlgoFactory::GetAlgo(
+  const GeomAbs_SurfaceType    theSurfaceType,
+  const IMeshTools_Parameters& theParameters) const
+{
+  switch (theSurfaceType)
+  {
+  case GeomAbs_Plane:
+    return theParameters.InternalVerticesMode ?
+      new NodeInsertionMeshAlgo<BRepMesh_DefaultRangeSplitter>::Type :
+      new BaseMeshAlgo::Type;
+    break;
+
+  case GeomAbs_Sphere:
+    {
+      NodeInsertionMeshAlgo<BRepMesh_SphereRangeSplitter>::Type* aMeshAlgo =
+        new NodeInsertionMeshAlgo<BRepMesh_SphereRangeSplitter>::Type;
+      aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True);
+      return aMeshAlgo;
+    }
+    break;
+
+  case GeomAbs_Cylinder:
+    return theParameters.InternalVerticesMode ?
+      new DefaultNodeInsertionMeshAlgo<BRepMesh_CylinderRangeSplitter>::Type :
+      new DefaultBaseMeshAlgo::Type;
+    break;
+
+  case GeomAbs_Cone:
+    {
+      NodeInsertionMeshAlgo<BRepMesh_ConeRangeSplitter>::Type* aMeshAlgo =
+        new NodeInsertionMeshAlgo<BRepMesh_ConeRangeSplitter>::Type;
+      aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True);
+      return aMeshAlgo;
+    }
+    break;
+
+  case GeomAbs_Torus:
+    {
+      NodeInsertionMeshAlgo<BRepMesh_TorusRangeSplitter>::Type* aMeshAlgo =
+        new NodeInsertionMeshAlgo<BRepMesh_TorusRangeSplitter>::Type;
+      aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True);
+      return aMeshAlgo;
+    }
+    break;
+
+  case GeomAbs_SurfaceOfRevolution:
+    {
+      DeflectionControlMeshAlgo<BRepMesh_BoundaryParamsRangeSplitter>::Type* aMeshAlgo =
+        new DeflectionControlMeshAlgo<BRepMesh_BoundaryParamsRangeSplitter>::Type;
+      aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True);
+      return aMeshAlgo;
+    }
+    break;
+
+  default:
+    {
+      DeflectionControlMeshAlgo<BRepMesh_NURBSRangeSplitter>::Type* aMeshAlgo =
+        new DeflectionControlMeshAlgo<BRepMesh_NURBSRangeSplitter>::Type;
+      aMeshAlgo->SetPreProcessSurfaceNodes (Standard_True);
+      return aMeshAlgo;
+    }
+  }
+}
diff --git a/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx b/src/BRepMesh/BRepMesh_DelabellaMeshAlgoFactory.hxx
new file mode 100644 (file)
index 0000000..486e0d5
--- /dev/null
@@ -0,0 +1,44 @@
+// Created on: 2019-07-05
+// Copyright (c) 2019 OPEN CASCADE SAS
+// Created by: Oleg AGASHIN
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BRepMesh_DelabellaMeshAlgoFactory_HeaderFile
+#define _BRepMesh_DelabellaMeshAlgoFactory_HeaderFile
+
+#include <Standard_Transient.hxx>
+#include <Standard_Type.hxx>
+#include <GeomAbs_SurfaceType.hxx>
+#include <IMeshTools_MeshAlgoFactory.hxx>
+
+//! Implementation of IMeshTools_MeshAlgoFactory providing Delabella-based
+//! algorithms of different complexity depending on type of target surface.
+class BRepMesh_DelabellaMeshAlgoFactory : public IMeshTools_MeshAlgoFactory
+{
+public:
+
+  //! Constructor.
+  Standard_EXPORT BRepMesh_DelabellaMeshAlgoFactory ();
+
+  //! Destructor.
+  Standard_EXPORT virtual ~BRepMesh_DelabellaMeshAlgoFactory ();
+
+  //! Creates instance of meshing algorithm for the given type of surface.
+  Standard_EXPORT virtual Handle(IMeshTools_MeshAlgo) GetAlgo(
+    const GeomAbs_SurfaceType    theSurfaceType,
+    const IMeshTools_Parameters& theParameters) const Standard_OVERRIDE;
+
+  DEFINE_STANDARD_RTTIEXT(BRepMesh_DelabellaMeshAlgoFactory, IMeshTools_MeshAlgoFactory)
+};
+
+#endif
index 5a8d30ffe50a10f75da82f214318a7d9848cc43d..6ca5cee049fa27eacbaf3d0f370e6e0fabb3303f 100644 (file)
@@ -90,9 +90,10 @@ BRepMesh_Delaun::BRepMesh_Delaun (
   const Standard_Boolean                        isFillCircles)
 : myMeshData ( theOldMesh ),
   myCircles (new NCollection_IncAllocator(
-             IMeshData::MEMORY_BLOCK_SIZE_HUGE))
+             IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
+  mySupVert (3),
+  myInitCircles (Standard_False)
 {
-  memset (mySupVert, 0, sizeof (mySupVert));
   if (isFillCircles)
   {
     InitCirclesTool (theCellsCountU, theCellsCountV);
@@ -105,9 +106,10 @@ BRepMesh_Delaun::BRepMesh_Delaun (
 //=======================================================================
 BRepMesh_Delaun::BRepMesh_Delaun(IMeshData::Array1OfVertexOfDelaun& theVertices)
 : myCircles (theVertices.Length(), new NCollection_IncAllocator(
-             IMeshData::MEMORY_BLOCK_SIZE_HUGE))
+             IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
+  mySupVert (3),
+  myInitCircles (Standard_False)
 {
-  memset (mySupVert, 0, sizeof (mySupVert));
   if ( theVertices.Length() > 2 )
   {
     myMeshData = new BRepMesh_DataStructureOfDelaun(
@@ -126,9 +128,10 @@ BRepMesh_Delaun::BRepMesh_Delaun(
   IMeshData::Array1OfVertexOfDelaun&            theVertices)
 : myMeshData( theOldMesh ),
   myCircles ( theVertices.Length(), new NCollection_IncAllocator(
-             IMeshData::MEMORY_BLOCK_SIZE_HUGE))
+             IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
+  mySupVert (3),
+  myInitCircles (Standard_False)
 {
-  memset (mySupVert, 0, sizeof (mySupVert));
   if ( theVertices.Length() > 2 )
   {
     Init( theVertices );
@@ -144,9 +147,10 @@ BRepMesh_Delaun::BRepMesh_Delaun(
   IMeshData::VectorOfInteger&                   theVertexIndices)
 : myMeshData( theOldMesh ),
   myCircles ( theVertexIndices.Length(), new NCollection_IncAllocator(
-             IMeshData::MEMORY_BLOCK_SIZE_HUGE))
+             IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
+  mySupVert (3),
+  myInitCircles (Standard_False)
 {
-  memset (mySupVert, 0, sizeof (mySupVert));
   perform(theVertexIndices);
 }
 
@@ -160,9 +164,10 @@ BRepMesh_Delaun::BRepMesh_Delaun (const Handle (BRepMesh_DataStructureOfDelaun)&
                                   const Standard_Integer                         theCellsCountV)
 : myMeshData (theOldMesh),
   myCircles (theVertexIndices.Length (), new NCollection_IncAllocator(
-             IMeshData::MEMORY_BLOCK_SIZE_HUGE))
+             IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
+  mySupVert (3),
+  myInitCircles (Standard_False)
 {
-  memset (mySupVert, 0, sizeof (mySupVert));
   perform (theVertexIndices, theCellsCountU, theCellsCountV);
 }
 
@@ -240,6 +245,8 @@ void BRepMesh_Delaun::initCirclesTool (const Bnd_Box2d&       theBox,
   myCircles.SetMinMaxSize( gp_XY( aMinX, aMinY ), gp_XY( aMaxX, aMaxY ) );
   myCircles.SetCellSize  ( aDeltaX / Max (theCellsCountU, aScaler),
                            aDeltaY / Max (theCellsCountV, aScaler));
+
+  myInitCircles = Standard_True;
 }
 
 //=======================================================================
@@ -290,14 +297,14 @@ void BRepMesh_Delaun::superMesh(const Bnd_Box2d& theBox)
   Standard_Real aDeltaMax = Max( aDeltaX, aDeltaY );
   Standard_Real aDelta    = aDeltaX + aDeltaY;
 
-  mySupVert[0] = myMeshData->AddNode(
-    BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) );
+  mySupVert.Append (myMeshData->AddNode(
+    BRepMesh_Vertex( ( aMinX + aMaxX ) / 2, aMaxY + aDeltaMax, BRepMesh_Free ) ) );
 
-  mySupVert[1] = myMeshData->AddNode(
-    BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
+  mySupVert.Append (myMeshData->AddNode(
+    BRepMesh_Vertex( aMinX - aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ) );
 
-  mySupVert[2] = myMeshData->AddNode(
-    BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) );
+  mySupVert.Append (myMeshData->AddNode(
+    BRepMesh_Vertex( aMaxX + aDelta, aMinY - aDeltaMin, BRepMesh_Free ) ) );
 
   Standard_Integer e[3];
   Standard_Boolean o[3];
@@ -324,7 +331,7 @@ void BRepMesh_Delaun::superMesh(const Bnd_Box2d& theBox)
 void BRepMesh_Delaun::deleteTriangle(const Standard_Integer          theIndex, 
                                      IMeshData::MapOfIntegerInteger& theLoopEdges )
 {
-  if (!myCircles.IsEmpty())
+  if (myInitCircles)
   {
     myCircles.Delete (theIndex);
   }
@@ -373,12 +380,25 @@ void BRepMesh_Delaun::compute(IMeshData::VectorOfInteger& theVertexIndexes)
     createTrianglesOnNewVertices (theVertexIndexes, Message_ProgressRange());
   }
 
+  RemoveAuxElements ();
+}
+
+//=======================================================================
+//function : RemoveAuxElements
+//purpose  :
+//=======================================================================
+void BRepMesh_Delaun::RemoveAuxElements ()
+{
+  Handle (NCollection_IncAllocator) aAllocator = new NCollection_IncAllocator (
+    IMeshData::MEMORY_BLOCK_SIZE_HUGE);
+
+  IMeshData::MapOfIntegerInteger aLoopEdges (10, aAllocator);
+
   // Destruction of triangles containing a top of the super triangle
-  BRepMesh_SelectorOfDataStructureOfDelaun aSelector( myMeshData );
-  for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
+  BRepMesh_SelectorOfDataStructureOfDelaun aSelector (myMeshData);
+  for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size(); ++aSupVertId)
     aSelector.NeighboursOfNode( mySupVert[aSupVertId] );
-  
-  aLoopEdges.Clear();
+
   IMeshData::IteratorOfMapOfInteger aFreeTriangles( aSelector.Elements() );
   for ( ; aFreeTriangles.More(); aFreeTriangles.Next() )
     deleteTriangle( aFreeTriangles.Key(), aLoopEdges );
@@ -393,7 +413,7 @@ void BRepMesh_Delaun::compute(IMeshData::VectorOfInteger& theVertexIndexes)
   }
 
   // The tops of the super triangle are destroyed
-  for (Standard_Integer aSupVertId = 0; aSupVertId < 3; ++aSupVertId)
+  for (Standard_Integer aSupVertId = 0; aSupVertId < mySupVert.Size (); ++aSupVertId)
     myMeshData->RemoveNode( mySupVert[aSupVertId] );
 }
 
@@ -786,9 +806,7 @@ void BRepMesh_Delaun::cleanupMesh()
               myMeshData->ElementNodes (aCurTriangle, v);
               for (int aNodeIdx = 0; aNodeIdx < 3 && isCanNotBeRemoved; ++aNodeIdx)
               {
-                if (v[aNodeIdx] == mySupVert[0] ||
-                    v[aNodeIdx] == mySupVert[1] ||
-                    v[aNodeIdx] == mySupVert[2])
+                if (isSupVertex (v[aNodeIdx]))
                 {
                   isCanNotBeRemoved = Standard_False;
                 }
@@ -1246,11 +1264,15 @@ void BRepMesh_Delaun::addTriangle( const Standard_Integer (&theEdgesId)[3],
     myMeshData->AddElement(BRepMesh_Triangle(theEdgesId, 
       theEdgesOri, BRepMesh_Free));
 
-  Standard_Boolean isAdded = myCircles.Bind( 
-    aNewTriangleId,
-    GetVertex( theNodesId[0] ).Coord(), 
-    GetVertex( theNodesId[1] ).Coord(),
-    GetVertex( theNodesId[2] ).Coord() );
+  Standard_Boolean isAdded = Standard_True;
+  if (myInitCircles)
+  {
+    isAdded = myCircles.Bind( 
+      aNewTriangleId,
+      GetVertex( theNodesId[0] ).Coord(), 
+      GetVertex( theNodesId[1] ).Coord(),
+      GetVertex( theNodesId[2] ).Coord() );
+  }
     
   if ( !isAdded )
     myMeshData->RemoveElement( aNewTriangleId );
index d60b76dc1834789c06712cd26d83b410f326e1fb..f0812d48d19ae68107209e27da704e2203585697 100755 (executable)
@@ -149,6 +149,17 @@ public:
                                              const Standard_Real    theSqTolerance,
                                              Standard_Integer&      theEdgeOn) const;
 
+  //! Explicitly sets ids of auxiliary vertices used to build mesh and used by 3rd-party algorithms.
+  inline void SetAuxVertices (const IMeshData::VectorOfInteger& theSupVert)
+  {
+    mySupVert = theSupVert;
+  }
+
+  //! Destruction of auxiliary triangles containing the given vertices.
+  //! Removes auxiliary vertices also.
+  //! @param theAuxVertices auxiliary vertices to be cleaned up.
+  Standard_EXPORT void RemoveAuxElements ();
+
 private:
 
   enum ReplaceFlag
@@ -353,13 +364,27 @@ private:
   //! Performs insertion of internal edges into mesh.
   void insertInternalEdges();
 
+  //! Checks whether the given vertex id relates to super contour.
+  Standard_Boolean isSupVertex (const Standard_Integer theVertexIdx) const
+  {
+    for (IMeshData::VectorOfInteger::Iterator aIt (mySupVert); aIt.More (); aIt.Next ())
+    {
+      if (theVertexIdx == aIt.Value ())
+      {
+        return Standard_True;
+      }
+    }
+
+    return Standard_False;
+  }
+
 private:
 
   Handle(BRepMesh_DataStructureOfDelaun) myMeshData;
   BRepMesh_CircleTool                    myCircles;
-  Standard_Integer                       mySupVert[3];
+  IMeshData::VectorOfInteger             mySupVert;
+  Standard_Boolean                       myInitCircles;
   BRepMesh_Triangle                      mySupTrian;
-
 };
 
 #endif
index 60793d7f4e7e2e76131404d744684b3e487c9d7c..627d489cf023be4835e89107aa3e976560f4a1d4 100644 (file)
@@ -23,7 +23,7 @@
 class BRepMesh_DataStructureOfDelaun;
 class BRepMesh_Delaun;
 
-//! Class provides base fuctionality to build face triangulation using Dealunay approach.
+//! Class provides base functionality to build face triangulation using Dealunay approach.
 //! Performs generation of mesh using raw data from model.
 class BRepMesh_DelaunayBaseMeshAlgo : public BRepMesh_ConstrainedBaseMeshAlgo
 {
index 128901f558ce60a5f907e39d5e37125efca7acd4..ea355e797fa088453a169c74203e6cd4cce88627 100644 (file)
@@ -46,7 +46,7 @@ public:
 
 protected:
 
-  //! Perfroms processing of generated mesh. Generates surface nodes and inserts them into structure.
+  //! Performs processing of generated mesh. Generates surface nodes and inserts them into structure.
   virtual void postProcessMesh (BRepMesh_Delaun& theMesher,
                                 const Message_ProgressRange& theRange) Standard_OVERRIDE
   {
index d9f46e9c46f6c354d31cb5f7fe2468a149d7b836..a2e6d5037ef5fd2d03d91566d1c97b7eb220ef1c 100644 (file)
@@ -84,7 +84,7 @@ protected:
                                           &this->getRangeSplitter());
   }
 
-  //! Perfroms processing of generated mesh. Generates surface nodes and inserts them into structure.
+  //! Performs processing of generated mesh. Generates surface nodes and inserts them into structure.
   virtual void postProcessMesh (BRepMesh_Delaun& theMesher,
                                 const Message_ProgressRange& theRange) Standard_OVERRIDE
   {
index cfae699e78c047ecafbf972e1369d354a9b084a4..eb00e303ebaa2efade713f2e1ed0aae56dc6c06a 100644 (file)
@@ -88,7 +88,7 @@ BRepMesh_IncrementalMesh::~BRepMesh_IncrementalMesh()
 //=======================================================================
 void BRepMesh_IncrementalMesh::Perform(const Message_ProgressRange& theRange)
 {
-  Handle(BRepMesh_Context) aContext = new BRepMesh_Context;
+  Handle(BRepMesh_Context) aContext = new BRepMesh_Context (myParameters.MeshAlgo);
   Perform (aContext, theRange);
 }
 
index 72a27c55877d25bb6f4d758d734375149ffb46e5..ed6b9b1db5640c7e9c98ca7b4311ad99d3d764ea 100644 (file)
@@ -22,7 +22,7 @@
 #include <IMeshTools_MeshAlgoFactory.hxx>
 
 //! Default implementation of IMeshTools_MeshAlgoFactory providing algorithms 
-//! of different compexity depending on type of target surface.
+//! of different complexity depending on type of target surface.
 class BRepMesh_MeshAlgoFactory : public IMeshTools_MeshAlgoFactory
 {
 public:
index a37a4504e43404e644f1bca6f0935dbc2f5f8c42..44333fcb2f4b47881340629476f5b94a2dcacbc4 100755 (executable)
@@ -86,3 +86,9 @@ BRepMesh_VertexTool.cxx
 BRepMesh_VertexTool.hxx
 BRepMesh_CustomBaseMeshAlgo.hxx
 BRepMesh_CustomDelaunayBaseMeshAlgo.hxx
+delabella.pxx
+delabella.cpp
+BRepMesh_DelabellaBaseMeshAlgo.hxx
+BRepMesh_DelabellaBaseMeshAlgo.cxx
+BRepMesh_DelabellaMeshAlgoFactory.hxx
+BRepMesh_DelabellaMeshAlgoFactory.cxx
diff --git a/src/BRepMesh/delabella.cpp b/src/BRepMesh/delabella.cpp
new file mode 100644 (file)
index 0000000..e88ccf2
--- /dev/null
@@ -0,0 +1,1043 @@
+/*
+
+MIT License
+
+DELABELLA - Delaunay triangulation library
+Copyright (C) 2018 GUMIX - Marcin Sokalski
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <search.h>
+
+#if (defined(__APPLE__))
+#include <malloc/malloc.h>
+#else
+#include <malloc.h>
+#endif
+
+#include <algorithm>
+#include "delabella.pxx" // we just need LOG() macro
+
+// assuming BITS is max(X_BITS,Y_BITS)
+
+typedef double Signed14;               // BITS                 xy coords
+typedef double Signed15;               // BITS + 1             vect::xy
+typedef long double Unsigned28; // 2xBITS              z coord
+typedef long double Signed29;   // 2xBITS + 1  vect::z
+typedef long double Signed31;   // 2xBITS + 3  norm::z
+typedef long double Signed45;   // 3xBITS + 3  norm::xy
+typedef long double Signed62;   // 4xBITS + 6  dot(vect,norm)
+
+/*
+// above typedefs can be used to configure delabella arithmetic types
+// in example, EXACT SOLVER (with xy coords limited to 14bits integers in range: +/-8192) 
+// could be achieved with following configuration:
+
+typedef int16_t  Signed14;             // BITS                 xy coords
+typedef int16_t  Signed15;             // BITS + 1             vect::xy
+typedef uint32_t Unsigned28;   // 2xBITS               z coord
+typedef int32_t  Signed29;             // 2xBITS + 1   vect::z
+typedef int32_t  Signed31;             // 2xBITS + 3   norm::z
+typedef int64_t  Signed45;             // 3xBITS + 3   norm::xy
+typedef int64_t  Signed62;             // 4xBITS + 6   dot(vect,norm)
+
+// alternatively, one could use some BigInt implementation
+// in order to expand xy range
+*/
+
+
+static Unsigned28 s14sqr(const Signed14& s)
+{
+       return (Unsigned28)((Signed29)s*s);
+}
+
+struct Norm
+{
+       Signed45 x;
+       Signed45 y;
+       Signed31 z;
+};
+
+struct Vect
+{
+       Signed15 x, y;
+       Signed29 z;
+
+       Norm cross (const Vect& v) const // cross prod
+       {
+               Norm n;
+               n.x = (Signed45)y*v.z - (Signed45)z*v.y;
+               n.y = (Signed45)z*v.x - (Signed45)x*v.z;
+               n.z = (Signed29)x*v.y - (Signed29)y*v.x;
+               return n;
+       }
+};
+
+struct CDelaBella : IDelaBella
+{
+       struct Face;
+
+       struct Vert : DelaBella_Vertex
+       {
+               Unsigned28 z;
+               Face* sew;
+
+               Vect operator - (const Vert& v) const // diff
+               {
+                       Vect d;
+                       d.x = (Signed15)x - (Signed15)v.x;
+                       d.y = (Signed15)y - (Signed15)v.y;
+                       d.z = (Signed29)z - (Signed29)v.z;
+                       return d;
+               }
+
+               static bool overlap(const Vert* v1, const Vert* v2)
+               {
+                       return v1->x == v2->x && v1->y == v2->y;
+               }
+
+               bool operator < (const Vert& v) const
+               {
+                       return u28cmp(this, &v) < 0;
+               }
+
+               static int u28cmp(const void* a, const void* b)
+               {
+                       Vert* va = (Vert*)a;
+                       Vert* vb = (Vert*)b;
+                       if (va->z < vb->z)
+                               return -1;
+                       if (va->z > vb->z)
+                               return 1;
+                       if (va->y < vb->y)
+                               return -1;
+                       if (va->y > vb->y)
+                               return 1;
+                       if (va->x < vb->x)
+                               return -1;
+                       if (va->x > vb->x)
+                               return 1;
+                       if (va->i < vb->i)
+                               return -1;
+                       if (va->i > vb->i)
+                               return 1;
+                       return 0;
+               }
+       };
+
+       struct Face : DelaBella_Triangle
+       {
+               Norm n;
+
+               static Face* Alloc(Face** from)
+               {
+                       Face* f = *from;
+                       *from = (Face*)f->next;
+                       f->next = 0;
+                       return f;
+               }
+
+               void Free(Face** to)
+               {
+                       next = *to;
+                       *to = this;
+               }
+
+               Face* Next(const Vert* p) const
+               {
+                       // return next face in same direction as face vertices are (cw/ccw)
+
+                       if (v[0] == p)
+                               return (Face*)f[1];
+                       if (v[1] == p)
+                               return (Face*)f[2];
+                       if (v[2] == p)
+                               return (Face*)f[0];
+                       return 0;
+               }
+
+               Signed62 dot(const Vert& p) const // dot
+               {
+                       Vect d = p - *(Vert*)v[0];
+                       return (Signed62)n.x * d.x + (Signed62)n.y * d.y + (Signed62)n.z * d.z;
+               }
+
+               Norm cross() const // cross of diffs
+               {
+                       return (*(Vert*)v[1] - *(Vert*)v[0]).cross(*(Vert*)v[2] - *(Vert*)v[0]);
+               }
+       };
+
+       Vert* vert_alloc;
+       Face* face_alloc;
+       int max_verts;
+       int max_faces;
+
+       Face* first_dela_face;
+       Face* first_hull_face;
+       Vert* first_hull_vert;
+
+       int inp_verts;
+       int out_verts;
+
+       int(*errlog_proc)(void* file, const char* fmt, ...);
+       void* errlog_file;
+
+       virtual ~CDelaBella ()
+       {
+       }
+
+       int Triangulate()
+       {
+               int points = inp_verts;
+               std::sort(vert_alloc, vert_alloc + points);
+
+               // rmove dups
+               {
+                       int w = 0, r = 1; // skip initial no-dups block
+                       while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + w))
+                       {
+                               w++;
+                               r++;
+                       }
+
+                       w++;
+
+                       while (r < points)
+                       {
+                               r++;
+
+                               // skip dups
+                               while (r < points && Vert::overlap(vert_alloc + r, vert_alloc + r - 1))
+                                       r++;
+
+                               // copy next no-dups block
+                               while (r < points && !Vert::overlap(vert_alloc + r, vert_alloc + r - 1))
+                                       vert_alloc[w++] = vert_alloc[r++];
+                       }
+
+                       if (points - w)
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[WRN] detected %d dups in xy array!\n", points - w);
+                               points = w;
+                       }
+               }
+
+               if (points < 3)
+               {
+                       if (points == 2)
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[WRN] all input points are colinear, returning single segment!\n");
+                               first_hull_vert = vert_alloc + 0;
+                               vert_alloc[0].next = (DelaBella_Vertex*)vert_alloc + 1;
+                               vert_alloc[1].next = 0;
+                       }
+                       else
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[WRN] all input points are identical, returning signle point!\n");
+                               first_hull_vert = vert_alloc + 0;
+                               vert_alloc[0].next = 0;
+                       }
+
+                       return -points;
+               }
+
+               int hull_faces = 2 * points - 4;
+
+               if (max_faces < hull_faces)
+               {
+                       if (max_faces)
+                               free(face_alloc);
+                       max_faces = 0;
+                       face_alloc = (Face*)malloc(sizeof(Face) * hull_faces);
+                       if (face_alloc)
+                               max_faces = hull_faces;
+                       else
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
+                               return 0;
+                       }
+               }
+
+               for (int i = 1; i < hull_faces; i++)
+                       face_alloc[i - 1].next = face_alloc + i;
+               face_alloc[hull_faces - 1].next = 0;
+
+               Face* cache = face_alloc;
+               Face* hull = 0;
+
+               Face f; // tmp
+               f.v[0] = vert_alloc + 0;
+               f.v[1] = vert_alloc + 1;
+               f.v[2] = vert_alloc + 2;
+               f.n = f.cross();
+
+               bool colinear = f.n.z == 0;
+               int i = 3;
+
+               /////////////////////////////////////////////////////////////////////////
+               // UNTIL INPUT IS COPLANAR, GROW IT IN FORM OF A 2D CONTOUR
+               /*
+               . |                |         after adding     . |        ________* L
+               . \ Last points to / Head     next point      . \ ______/        /
+               .  *____          |             ----->        .H *____          |
+               .  |\_  \_____    |                           .  |\_  \_____    |
+               .   \ \_      \__* - Tail points to Last      .   \ \_      \__* T
+               .    \  \_      /                             .    \  \_      /
+               .     \__ \_ __/                              .     \__ \_ __/
+               .        \__* - Head points to Tail           .        \__/
+               */
+
+               Vert* head = (Vert*)f.v[0];
+               Vert* tail = (Vert*)f.v[1];
+               Vert* last = (Vert*)f.v[2];
+
+               head->next = tail;
+               tail->next = last;
+               last->next = head;
+
+               while (i < points && f.dot(vert_alloc[i]) == 0)
+               {
+                       Vert* v = vert_alloc + i;
+
+                       // it is enough to test just 1 non-zero coord
+                       // but we want also to test stability (assert) 
+                       // so we calc all signs...
+
+                       // why not testing sign of dot prod of 2 normals?
+                       // that way we'd fall into precission problems
+
+                       Norm LvH = (*v - *last).cross(*head - *last);
+                       bool lvh =
+                               (f.n.x > 0 && LvH.x > 0) || (f.n.x < 0 && LvH.x < 0) ||
+                               (f.n.y > 0 && LvH.y > 0) || (f.n.y < 0 && LvH.y < 0) ||
+                               (f.n.z > 0 && LvH.z > 0) || (f.n.z < 0 && LvH.z < 0);
+
+                       Norm TvL = (*v - *tail).cross(*last - *tail);
+                       bool tvl =
+                               (f.n.x > 0 && TvL.x > 0) || (f.n.x < 0 && TvL.x < 0) ||
+                               (f.n.y > 0 && TvL.y > 0) || (f.n.y < 0 && TvL.y < 0) ||
+                               (f.n.z > 0 && TvL.z > 0) || (f.n.z < 0 && TvL.z < 0);
+
+                       if (lvh && !tvl) // insert new f on top of e(2,0) = (last,head)
+                       {
+                               // f.v[0] = head;
+                               f.v[1] = last;
+                               f.v[2] = v;
+
+                               last->next = v;
+                               v->next = head;
+                               tail = last;
+                       }
+                       else
+                               if (tvl && !lvh) // insert new f on top of e(1,2) = (tail,last)
+                               {
+                                       f.v[0] = last;
+                                       //f.v[1] = tail;
+                                       f.v[2] = v;
+
+                                       tail->next = v;
+                                       v->next = last;
+                                       head = last;
+                               }
+                               else
+                               {
+                                       // wtf? dilithium crystals are fucked.
+                                       assert(0);
+                               }
+
+                       last = v;
+                       i++;
+               }
+
+               if (i == points)
+               {
+                       if (colinear)
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[WRN] all input points are colinear, returning segment list!\n");
+                               first_hull_vert = head;
+                               last->next = 0; // break contour, make it a list 
+                               return -points;
+                       }
+                       else
+                       {
+                               if (points > 3)
+                               {
+                                       if (errlog_proc)
+                                               errlog_proc(errlog_file, "[NFO] all input points are cocircular.\n");
+                               }
+                               else
+                               {
+                                       if (errlog_proc)
+                                               errlog_proc(errlog_file, "[NFO] trivial case of 3 points, thank you.\n");
+
+                                       first_dela_face = Face::Alloc(&cache);
+                                       first_dela_face->next = 0;
+                                       first_hull_face = Face::Alloc(&cache);
+                                       first_hull_face->next = 0;
+
+                                       first_dela_face->f[0] = first_dela_face->f[1] = first_dela_face->f[2] = first_hull_face;
+                                       first_hull_face->f[0] = first_hull_face->f[1] = first_hull_face->f[2] = first_dela_face;
+
+                                       head->sew = tail->sew = last->sew = first_hull_face;
+
+                                       if (f.n.z < 0)
+                                       {
+                                               first_dela_face->v[0] = head;
+                                               first_dela_face->v[1] = tail;
+                                               first_dela_face->v[2] = last;
+                                               first_hull_face->v[0] = last;
+                                               first_hull_face->v[1] = tail;
+                                               first_hull_face->v[2] = head;
+
+                                               // reverse silhouette
+                                               head->next = last;
+                                               last->next = tail;
+                                               tail->next = head;
+
+                                               first_hull_vert = last;
+                                       }
+                                       else
+                                       {
+                                               first_dela_face->v[0] = last;
+                                               first_dela_face->v[1] = tail;
+                                               first_dela_face->v[2] = head;
+                                               first_hull_face->v[0] = head;
+                                               first_hull_face->v[1] = tail;
+                                               first_hull_face->v[2] = last;
+
+                                               first_hull_vert = head;
+                                       }
+
+                                       first_dela_face->n = first_dela_face->cross();
+                                       first_hull_face->n = first_hull_face->cross();
+
+                                       return 3;
+                               }
+
+                               // retract last point it will be added as a cone's top later
+                               last = head;
+                               head = (Vert*)head->next;
+                               i--;
+                       }
+               }
+
+               /////////////////////////////////////////////////////////////////////////
+               // CREATE CONE HULL WITH TOP AT cloud[i] AND BASE MADE OF CONTOUR LIST
+               // in 2 ways :( - depending on at which side of the contour a top vertex appears
+
+               Vert* q = vert_alloc + i;
+
+               if (f.dot(*q) > 0)
+               {
+                       Vert* p = last;
+                       Vert* n = (Vert*)p->next;
+
+                       Face* first_side = Face::Alloc(&cache);
+                       first_side->v[0] = p;
+                       first_side->v[1] = n;
+                       first_side->v[2] = q;
+                       first_side->n = first_side->cross();
+                       hull = first_side;
+
+                       p = n;
+                       n = (Vert*)n->next;
+
+                       Face* prev_side = first_side;
+                       Face* prev_base = 0;
+                       Face* first_base = 0;
+
+                       do
+                       {
+                               Face* base = Face::Alloc(&cache);
+                               base->v[0] = n;
+                               base->v[1] = p;
+                               base->v[2] = last;
+                               base->n = base->cross();
+
+                               Face* side = Face::Alloc(&cache);
+                               side->v[0] = p;
+                               side->v[1] = n;
+                               side->v[2] = q;
+                               side->n = side->cross();
+
+                               side->f[2] = base;
+                               base->f[2] = side;
+
+                               side->f[1] = prev_side;
+                               prev_side->f[0] = side;
+
+                               base->f[0] = prev_base;
+                               if (prev_base)
+                                       prev_base->f[1] = base;
+                               else
+                                       first_base = base;
+
+                               prev_base = base;
+                               prev_side = side;
+
+                               p = n;
+                               n = (Vert*)n->next;
+                       } while (n != last);
+
+                       Face* last_side = Face::Alloc(&cache);
+                       last_side->v[0] = p;
+                       last_side->v[1] = n;
+                       last_side->v[2] = q;
+                       last_side->n = last_side->cross();
+
+                       last_side->f[1] = prev_side;
+                       prev_side->f[0] = last_side;
+
+                       last_side->f[0] = first_side;
+                       first_side->f[1] = last_side;
+
+                       first_base->f[0] = first_side;
+                       first_side->f[2] = first_base;
+
+                       last_side->f[2] = prev_base;
+                       prev_base->f[1] = last_side;
+
+                       i++;
+               }
+               else
+               {
+                       Vert* p = last;
+                       Vert* n = (Vert*)p->next;
+
+                       Face* first_side = Face::Alloc(&cache);
+                       first_side->v[0] = n;
+                       first_side->v[1] = p;
+                       first_side->v[2] = q;
+                       first_side->n = first_side->cross();
+                       hull = first_side;
+
+                       p = n;
+                       n = (Vert*)n->next;
+
+                       Face* prev_side = first_side;
+                       Face* prev_base = 0;
+                       Face* first_base = 0;
+
+                       do
+                       {
+                               Face* base = Face::Alloc(&cache);
+                               base->v[0] = p;
+                               base->v[1] = n;
+                               base->v[2] = last;
+                               base->n = base->cross();
+
+                               Face* side = Face::Alloc(&cache);
+                               side->v[0] = n;
+                               side->v[1] = p;
+                               side->v[2] = q;
+                               side->n = side->cross();
+
+                               side->f[2] = base;
+                               base->f[2] = side;
+
+                               side->f[0] = prev_side;
+                               prev_side->f[1] = side;
+
+                               base->f[1] = prev_base;
+                               if (prev_base)
+                                       prev_base->f[0] = base;
+                               else
+                                       first_base = base;
+
+                               prev_base = base;
+                               prev_side = side;
+
+                               p = n;
+                               n = (Vert*)n->next;
+                       } while (n != last);
+
+                       Face* last_side = Face::Alloc(&cache);
+                       last_side->v[0] = n;
+                       last_side->v[1] = p;
+                       last_side->v[2] = q;
+                       last_side->n = last_side->cross();
+
+                       last_side->f[0] = prev_side;
+                       prev_side->f[1] = last_side;
+
+                       last_side->f[1] = first_side;
+                       first_side->f[0] = last_side;
+
+                       first_base->f[1] = first_side;
+                       first_side->f[2] = first_base;
+
+                       last_side->f[2] = prev_base;
+                       prev_base->f[0] = last_side;
+
+                       i++;
+               }
+
+               /////////////////////////////////////////////////////////////////////////
+               // ACTUAL ALGORITHM
+
+               for (; i < points; i++)
+               {
+                       //ValidateHull(alloc, 2 * i - 4);
+
+                       Vert* _q = vert_alloc + i;
+                       Vert* _p = vert_alloc + i - 1;
+                       Face* _f = hull;
+
+                       // 1. FIND FIRST VISIBLE FACE 
+                       //    simply iterate around last vertex using last added triange adjecency info
+                       while (_f->dot(*_q) <= 0)
+                       {
+                               _f = _f->Next(_p);
+                               if (_f == hull)
+                               {
+                                       // if no visible face can be located at last vertex,
+                                       // let's run through all faces (approximately last to first), 
+                                       // yes this is emergency fallback and should not ever happen.
+                                       _f = face_alloc + 2 * i - 4 - 1;
+                                       while (_f->dot(*_q) <= 0)
+                                       {
+                                               assert(_f != face_alloc); // no face is visible? you must be kidding!
+                                               _f--;
+                                       }
+                               }
+                       }
+
+                       // 2. DELETE VISIBLE FACES & ADD NEW ONES
+                       //    (we also build silhouette (vertex loop) between visible & invisible faces)
+
+                       int del = 0;
+                       int add = 0;
+
+                       // push first visible face onto stack (of visible faces)
+                       Face* stack = _f;
+                       _f->next = _f; // old trick to use list pointers as 'on-stack' markers
+                       while (stack)
+                       {
+                               // pop, take care of last item ptr (it's not null!)
+                               _f = stack;
+                               stack = (Face*)_f->next;
+                               if (stack == _f)
+                                       stack = 0;
+                               _f->next = 0;
+
+                               // copy parts of old face that we still need after removal
+                               Vert* fv[3] = { (Vert*)_f->v[0],(Vert*)_f->v[1],(Vert*)_f->v[2] };
+                               Face* ff[3] = { (Face*)_f->f[0],(Face*)_f->f[1],(Face*)_f->f[2] };
+
+                               // delete visible face
+                               _f->Free(&cache);
+                               del++;
+
+                               // check all 3 neighbors
+                               for (int e = 0; e < 3; e++)
+                               {
+                                       Face* n = ff[e];
+                                       if (n && !n->next) // ensure neighbor is not processed yet & isn't on stack 
+                                       {
+                                               if (n->dot(*_q) <= 0) // if neighbor is not visible we have slihouette edge
+                                               {
+                                                       // build face
+                                                       add++;
+
+                                                       // ab: given face adjacency [index][], 
+                                                       // it provides [][2] vertex indices on shared edge (CCW order)
+                                                       const static int ab[3][2] = { { 1,2 },{ 2,0 },{ 0,1 } };
+
+                                                       Vert* a = fv[ab[e][0]];
+                                                       Vert* b = fv[ab[e][1]];
+
+                                                       Face* s = Face::Alloc(&cache);
+                                                       s->v[0] = a;
+                                                       s->v[1] = b;
+                                                       s->v[2] = _q;
+
+                                                       s->n = s->cross();
+                                                       s->f[2] = n;
+
+                                                       // change neighbour's adjacency from old visible face to cone side
+                                                       if (n->f[0] == _f)
+                                                               n->f[0] = s;
+                                                       else
+                                                               if (n->f[1] == _f)
+                                                                       n->f[1] = s;
+                                                               else
+                                                                       if (n->f[2] == _f)
+                                                                               n->f[2] = s;
+                                                                       else
+                                                                               assert(0);
+
+                                                       // build silhouette needed for sewing sides in the second pass
+                                                       a->sew = s;
+                                                       a->next = b;
+                                               }
+                                               else
+                                               {
+                                                       // disjoin visible faces
+                                                       // so they won't be processed more than once
+
+                                                       if (n->f[0] == _f)
+                                                               n->f[0] = 0;
+                                                       else
+                                                               if (n->f[1] == _f)
+                                                                       n->f[1] = 0;
+                                                               else
+                                                                       if (n->f[2] == _f)
+                                                                               n->f[2] = 0;
+                                                                       else
+                                                                               assert(0);
+
+                                                       // push neighbor face, it's visible and requires processing
+                                                       n->next = stack ? stack : n;
+                                                       stack = n;
+                                               }
+                                       }
+                               }
+                       }
+
+                       // if add<del+2 hungry hull has consumed some point
+                       // that means we can't do delaunay for some under precission reasons
+                       // althought convex hull would be fine with it
+                       assert(add == del + 2);
+
+                       // 3. SEW SIDES OF CONE BUILT ON SLIHOUTTE SEGMENTS
+
+                       hull = face_alloc + 2 * i - 4 + 1; // last added face
+
+                                                                                 // last face must contain part of the silhouette
+                                                                                 // (edge between its v[0] and v[1])
+                       Vert* entry = (Vert*)hull->v[0];
+
+                       Vert* pr = entry;
+                       do
+                       {
+                               // sew pr<->nx
+                               Vert* nx = (Vert*)pr->next;
+                               pr->sew->f[0] = nx->sew;
+                               nx->sew->f[1] = pr->sew;
+                               pr = nx;
+                       } while (pr != entry);
+               }
+
+               assert(2 * i - 4 == hull_faces);
+               //ValidateHull(alloc, hull_faces);
+
+               // needed?
+               for (i = 0; i < points; i++)
+               {
+                       vert_alloc[i].next = 0;
+                       vert_alloc[i].sew = 0;
+               }
+
+               i = 0;
+               Face** prev_dela = &first_dela_face;
+               Face** prev_hull = &first_hull_face;
+               for (int j = 0; j < hull_faces; j++)
+               {
+                       Face* _f = face_alloc + j;
+                       if (_f->n.z < 0)
+                       {
+                               *prev_dela = _f;
+                               prev_dela = (Face**)&_f->next;
+                               i++;
+                       }
+                       else
+                       {
+                               *prev_hull = _f;
+                               prev_hull = (Face**)&_f->next;
+                               if (((Face*)_f->f[0])->n.z < 0)
+                               {
+                                       _f->v[1]->next = _f->v[2];
+                                       ((Vert*)_f->v[1])->sew = _f;
+                               }
+                               if (((Face*)_f->f[1])->n.z < 0)
+                               {
+                                       _f->v[2]->next = _f->v[0];
+                                       ((Vert*)_f->v[2])->sew = _f;
+                               }
+                               if (((Face*)_f->f[2])->n.z < 0)
+                               {
+                                       _f->v[0]->next = _f->v[1];
+                                       ((Vert*)_f->v[0])->sew = _f;
+                               }
+                       }
+               }
+
+               *prev_dela = 0;
+               *prev_hull = 0;
+
+               first_hull_vert = (Vert*)first_hull_face->v[0];
+
+               // todo link slihouette verts into contour
+               // and sew its edges with hull faces
+
+               return 3*i;
+       }
+
+       bool ReallocVerts(int points)
+       {
+               inp_verts = points;
+               out_verts = 0;
+
+               first_dela_face = 0;
+               first_hull_face = 0;
+               first_hull_vert = 0;
+
+               if (max_verts < points)
+               {
+                       if (max_verts)
+                       {
+                               free(vert_alloc);
+                               vert_alloc = 0;
+                               max_verts = 0;
+                       }
+
+                       vert_alloc = (Vert*)malloc(sizeof(Vert)*points);
+                       if (vert_alloc)
+                               max_verts = points;
+                       else
+                       {
+                               if (errlog_proc)
+                                       errlog_proc(errlog_file, "[ERR] Not enough memory, shop for some more RAM. See you!\n");
+                               return false;
+                       }
+               }
+
+               return true;
+       }
+
+       virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0)
+       {
+               if (!x)
+                       return 0;
+
+               if (!y)
+                       y = x + 1;
+
+               if (advance_bytes < static_cast<int>(sizeof(float) * 2))
+                       advance_bytes = static_cast<int>(sizeof(float) * 2);
+
+               if (!ReallocVerts(points))
+                       return 0;
+
+               for (int i = 0; i < points; i++)
+               {
+                       Vert* v = vert_alloc + i;
+                       v->i = i;
+                       v->x = (Signed14)*(const float*)((const char*)x + i*advance_bytes);
+                       v->y = (Signed14)*(const float*)((const char*)y + i*advance_bytes);
+                       v->z = s14sqr(v->x) + s14sqr(v->y);
+               }
+
+               out_verts = Triangulate();
+               return out_verts;
+       }
+
+       virtual int Triangulate(int points, const double* x, const double* y, int advance_bytes)
+       {
+               if (!x)
+                       return 0;
+               
+               if (!y)
+                       y = x + 1;
+               
+               if (advance_bytes < static_cast<int>(sizeof(double) * 2))
+                       advance_bytes = static_cast<int>(sizeof(double) * 2);
+
+               if (!ReallocVerts(points))
+                       return 0;
+
+               for (int i = 0; i < points; i++)
+               {
+                       Vert* v = vert_alloc + i;
+                       v->i = i;
+                       v->x = (Signed14)*(const double*)((const char*)x + i*advance_bytes);
+                       v->y = (Signed14)*(const double*)((const char*)y + i*advance_bytes);
+                       v->z = s14sqr (v->x) + s14sqr (v->y);
+               }
+
+               out_verts = Triangulate();
+               return out_verts;
+       }
+
+       virtual void Destroy()
+       {
+               if (face_alloc)
+                       free(face_alloc);
+               if (vert_alloc)
+                       free(vert_alloc);
+               delete this;
+       }
+
+       // num of points passed to last call to Triangulate()
+       virtual int GetNumInputPoints() const
+       {
+               return inp_verts;
+       }
+
+       // num of verts returned from last call to Triangulate()
+       virtual int GetNumOutputVerts() const
+       {
+               return out_verts;
+       }
+
+       virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const
+       {
+               return first_dela_face;
+       }
+
+       virtual const DelaBella_Triangle* GetFirstHullTriangle() const
+       {
+               return first_hull_face;
+       }
+
+       virtual const DelaBella_Vertex* GetFirstHullVertex() const
+       {
+               return first_hull_vert;
+       }
+
+       virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream)
+       {
+               errlog_proc = proc;
+               errlog_file = stream;
+       }
+};
+
+IDelaBella* IDelaBella::Create()
+{
+       CDelaBella* db = new CDelaBella;
+       if (!db)
+               return 0;
+
+       db->vert_alloc = 0;
+       db->face_alloc = 0;
+       db->max_verts = 0;
+       db->max_faces = 0;
+
+       db->first_dela_face = 0;
+       db->first_hull_face = 0;
+       db->first_hull_vert = 0;
+
+       db->inp_verts = 0;
+       db->out_verts = 0;
+
+       db->errlog_proc = 0;
+       db->errlog_file = 0;
+
+       return db;
+}
+
+void* DelaBella_Create()
+{
+       return IDelaBella::Create();
+}
+
+void  DelaBella_Destroy(void* db)
+{
+       ((IDelaBella*)db)->Destroy();
+}
+
+void  DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream)
+{
+       ((IDelaBella*)db)->SetErrLog(proc, stream);
+}
+
+int   DelaBella_TriangulateFloat(void* db, int points, float* x, float* y, int advance_bytes)
+{
+       return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
+}
+
+int   DelaBella_TriangulateDouble(void* db, int points, double* x, double* y, int advance_bytes)
+{
+       return ((IDelaBella*)db)->Triangulate(points, x, y, advance_bytes);
+}
+
+int   DelaBella_GetNumInputPoints(void* db)
+{
+       return ((IDelaBella*)db)->GetNumInputPoints();
+}
+
+int   DelaBella_GetNumOutputVerts(void* db)
+{
+       return ((IDelaBella*)db)->GetNumOutputVerts();
+}
+
+const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db)
+{
+       return ((IDelaBella*)db)->GetFirstDelaunayTriangle();
+}
+
+const DelaBella_Triangle* GetFirstHullTriangle(void* db)
+{
+       return ((IDelaBella*)db)->GetFirstHullTriangle();
+}
+
+const DelaBella_Vertex*   GetFirstHullVertex(void* db)
+{
+       return ((IDelaBella*)db)->GetFirstHullVertex();
+}
+
+// depreciated!
+int DelaBella(int points, const double* xy, int* abc, int(*errlog)(const char* fmt, ...))
+{
+       if (errlog)
+               errlog("[WRN] Depreciated interface! errlog disabled.\n");
+
+       if (!xy || points <= 0)
+               return 0;
+
+       IDelaBella* db = IDelaBella::Create();
+       int verts = db->Triangulate(points, xy, 0, 0);
+
+       if (!abc)
+               return verts;
+
+       if (verts > 0)
+       {
+               int tris = verts / 3;
+               const DelaBella_Triangle* dela = db->GetFirstDelaunayTriangle();
+               for (int i = 0; i < tris; i++)
+               {
+                       for (int j=0; j<3; j++)
+                               abc[3 * i + j] = dela->v[j]->i;
+                       dela = dela->next;
+               }
+       }
+       else
+       {
+               int pnts = -verts;
+               const DelaBella_Vertex* line = db->GetFirstHullVertex();
+               for (int i = 0; i < pnts; i++)
+               {
+                       abc[i] = line->i;
+                       line = line->next;
+               }
+       }
+
+       return verts;
+}
diff --git a/src/BRepMesh/delabella.pxx b/src/BRepMesh/delabella.pxx
new file mode 100644 (file)
index 0000000..a7ff00b
--- /dev/null
@@ -0,0 +1,90 @@
+/*
+
+MIT License
+
+DELABELLA - Delaunay triangulation library
+Copyright (C) 2018 GUMIX - Marcin Sokalski
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#ifndef DELABELLA_H
+#define DELABELLA_H
+
+// returns: positive value: number of triangle indices, negative: number of line segment indices (degenerated input)
+//          triangle indices in abc array are always returned in clockwise order
+// DEPRECIATED. move to new API either extern "C" or IDelaBella (C++)
+int DelaBella(int points, const double* xy/*[points][2]*/, int* abc/*[2*points-5][3]*/, int (*errlog)(const char* fmt,...) = printf);
+
+struct DelaBella_Vertex
+{
+       int i; // index of original point
+       double x, y; // coordinates (input copy)
+       DelaBella_Vertex* next; // next silhouette vertex
+};
+
+struct DelaBella_Triangle
+{
+       DelaBella_Vertex* v[3]; // 3 vertices spanning this triangle
+       DelaBella_Triangle* f[3]; // 3 adjacent faces, f[i] is at the edge opposite to vertex v[i]
+       DelaBella_Triangle* next; // next triangle (of delaunay set or hull set)
+};
+
+#ifdef __cplusplus
+struct IDelaBella
+{
+       static IDelaBella* Create();
+
+       virtual void Destroy() = 0;
+
+       virtual void SetErrLog(int(*proc)(void* stream, const char* fmt, ...), void* stream) = 0;
+
+       // return 0: no output 
+       // negative: all points are colinear, output hull vertices form colinear segment list, no triangles on output
+       // positive: output hull vertices form counter-clockwise ordered segment contour, delaunay and hull triangles are available
+       // if 'y' pointer is null, y coords are treated to be located immediately after every x
+       // if advance_bytes is less than 2*sizeof coordinate type, it is treated as 2*sizeof coordinate type  
+       virtual int Triangulate(int points, const float* x, const float* y = 0, int advance_bytes = 0) = 0;
+       virtual int Triangulate(int points, const double* x, const double* y = 0, int advance_bytes = 0) = 0;
+
+       // num of points passed to last call to Triangulate()
+       virtual int GetNumInputPoints() const = 0;
+
+       // num of verts returned from last call to Triangulate()
+       virtual int GetNumOutputVerts() const = 0;
+
+       virtual const DelaBella_Triangle* GetFirstDelaunayTriangle() const = 0; // valid only if Triangulate() > 0
+       virtual const DelaBella_Triangle* GetFirstHullTriangle() const = 0; // valid only if Triangulate() > 0
+       virtual const DelaBella_Vertex*   GetFirstHullVertex() const = 0; // if Triangulate() < 0 it is list, otherwise closed contour! 
+};
+#else
+void* DelaBella_Create();
+void  DelaBella_Destroy(void* db);
+void  DelaBella_SetErrLog(void* db, int(*proc)(void* stream, const char* fmt, ...), void* stream);
+int   DelaBella_TriangulateFloat(void* db, int points, float* x, float* y = 0, int advance_bytes = 0);
+int   DelaBella_TriangulateDouble(void* db, int points, double* x, double* y = 0, int advance_bytes = 0);
+int   DelaBella_GetNumInputPoints(void* db);
+int   DelaBella_GetNumOutputVerts(void* db);
+const DelaBella_Triangle* GetFirstDelaunayTriangle(void* db);
+const DelaBella_Triangle* GetFirstHullTriangle(void* db);
+const DelaBella_Vertex*   GetFirstHullVertex(void* db);
+#endif
+
+#endif
index e582bc1d1b462e41b1f5252a5d58d4e7af538bfc..5c4537202c98da7574d53338bd72acb36c1de9e1 100644 (file)
@@ -2,6 +2,7 @@ IMeshTools_Context.hxx
 IMeshTools_CurveTessellator.hxx
 IMeshTools_MeshAlgo.hxx
 IMeshTools_MeshAlgoFactory.hxx
+IMeshTools_MeshAlgoType.hxx
 IMeshTools_MeshBuilder.hxx
 IMeshTools_MeshBuilder.cxx
 IMeshTools_ModelAlgo.hxx
diff --git a/src/IMeshTools/IMeshTools_MeshAlgoType.hxx b/src/IMeshTools/IMeshTools_MeshAlgoType.hxx
new file mode 100644 (file)
index 0000000..ca71c6c
--- /dev/null
@@ -0,0 +1,25 @@
+// Copyright (c) 2020 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _IMeshTools_MeshAlgoType_HeaderFile
+#define _IMeshTools_MeshAlgoType_HeaderFile
+
+//! Enumerates built-in meshing algorithms factories implementing IMeshTools_MeshAlgoFactory interface.
+enum IMeshTools_MeshAlgoType
+{
+  IMeshTools_MeshAlgoType_DEFAULT = -1, //!< use global default (IMeshTools_MeshAlgoType_Watson or CSF_MeshAlgo)
+  IMeshTools_MeshAlgoType_Watson  = 0,  //!< generate 2D Delaunay triangulation based on Watson algorithm (BRepMesh_MeshAlgoFactory)
+  IMeshTools_MeshAlgoType_Delabella,    //!< generate 2D Delaunay triangulation based on Delabella algorithm (BRepMesh_DelabellaMeshAlgoFactory)
+};
+
+#endif
index 17c35911132ce71c8c72360be30b5328ed423b3c..4bbff6a902f6b8b8ee988c6ac6d7c501ca855e83 100644 (file)
@@ -16,6 +16,7 @@
 #ifndef _IMeshTools_Parameters_HeaderFile
 #define _IMeshTools_Parameters_HeaderFile
 
+#include <IMeshTools_MeshAlgoType.hxx>
 #include <Precision.hxx>
 
 //! Structure storing meshing parameters
@@ -24,6 +25,7 @@ struct IMeshTools_Parameters {
   //! Default constructor
   IMeshTools_Parameters ()
     :
+    MeshAlgo (IMeshTools_MeshAlgoType_DEFAULT),
     Angle(0.5),
     Deflection(0.001),
     AngleInterior(-1.0),
@@ -47,6 +49,9 @@ struct IMeshTools_Parameters {
     return 0.1;
   }
 
+  //! 2D Delaunay triangulation algorithm factory to use
+  IMeshTools_MeshAlgoType                          MeshAlgo;
+
   //! Angular deflection used to tessellate the boundary edges
   Standard_Real                                    Angle;
 
index 5d337dee444562b76ccad3ee14fec136bcb15882..d02cfcd917d8d444eeeb024e5f50cfe18cf7bca7 100644 (file)
 #include <TopExp_Explorer.hxx>
 #include <TopTools_MapIteratorOfMapOfShape.hxx>
 
+#include <BRepMesh_Context.hxx>
+#include <BRepMesh_FaceDiscret.hxx>
+#include <BRepMesh_MeshAlgoFactory.hxx>
+#include <BRepMesh_DelabellaMeshAlgoFactory.hxx>
+
 //epa Memory leaks test
 //OAN: for triepoints
 #ifdef _WIN32
@@ -94,7 +99,8 @@ options:\n\
         -adjust_min     enables local adjustment of min size depending on edge size (switched off by default)\n\
         -force_face_def disables usage of shape tolerances for computing face deflection (switched off by default)\n\
         -decrease       enforces the meshing of the shape even if current mesh satisfies the new criteria\
-                        (switched off by default).\n";
+                        (switched off by default).\n\
+        -algo {watson|delabella} changes core triangulation algorithm to one with specified id (watson is used by default)\n";
     return 0;
   }
 
@@ -109,6 +115,7 @@ options:\n\
   aMeshParams.Deflection = aMeshParams.DeflectionInterior = 
     Max(Draw::Atof(argv[2]), Precision::Confusion());
 
+  Handle (IMeshTools_Context) aContext = new BRepMesh_Context;
   if (nbarg > 3)
   {
     Standard_Integer i = 3;
@@ -135,23 +142,54 @@ options:\n\
         aMeshParams.AllowQualityDecrease = Standard_True;
       else if (i < nbarg)
       {
-        Standard_Real aVal = Draw::Atof(argv[i++]);
-        if (aOpt == "-a")
-        {
-          aMeshParams.Angle = aVal * M_PI / 180.;
-        }
-        else if (aOpt == "-ai")
+        if (aOpt == "-algo")
         {
-          aMeshParams.AngleInterior = aVal * M_PI / 180.;
+          TCollection_AsciiString anAlgoStr (argv[i++]);
+          anAlgoStr.LowerCase();
+          if (anAlgoStr == "watson"
+           || anAlgoStr == "0")
+          {
+            aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Watson;
+            aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_MeshAlgoFactory));
+          }
+          else if (anAlgoStr == "delabella"
+                || anAlgoStr == "1")
+          {
+            aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Delabella;
+            aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_DelabellaMeshAlgoFactory));
+          }
+          else if (anAlgoStr == "-1"
+                || anAlgoStr == "default")
+          {
+            // already handled by BRepMesh_Context constructor
+            //aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_DEFAULT;
+          }
+          else
+          {
+            di << "Syntax error at " << anAlgoStr;
+            return 1;
+          }
         }
-        else if (aOpt == "-min")
-          aMeshParams.MinSize = aVal;
-        else if (aOpt == "-di")
+        else
         {
-          aMeshParams.DeflectionInterior = aVal;
+          Standard_Real aVal = Draw::Atof(argv[i++]);
+          if (aOpt == "-a")
+          {
+            aMeshParams.Angle = aVal * M_PI / 180.;
+          }
+          else if (aOpt == "-ai")
+          {
+            aMeshParams.AngleInterior = aVal * M_PI / 180.;
+          }
+          else if (aOpt == "-min")
+            aMeshParams.MinSize = aVal;
+          else if (aOpt == "-di")
+          {
+            aMeshParams.DeflectionInterior = aVal;
+          }
+          else
+            --i;
         }
-        else
-          --i;
       }
     }
   }
@@ -160,7 +198,11 @@ options:\n\
      << (aMeshParams.InParallel ? "ON" : "OFF") << "\n";
 
   Handle(Draw_ProgressIndicator) aProgress = new Draw_ProgressIndicator(di, 1);
-  BRepMesh_IncrementalMesh aMesher (aShape, aMeshParams, aProgress->Start());
+  BRepMesh_IncrementalMesh aMesher;
+  aMesher.SetShape (aShape);
+  aMesher.ChangeParameters() = aMeshParams;
+
+  aMesher.Perform (aContext, aProgress->Start());
 
   di << "Meshing statuses: ";
   const Standard_Integer aStatus = aMesher.GetStatusFlags();
index 9e70d381663e4a96dc3b6ab680937c0474b67b91..b9a40d69d2ef9316f7bd75384e27239ee85dc27e 100755 (executable)
@@ -133,3 +133,9 @@ int Sprintf (char* theBuffer, const char* theFormat, ...)
   va_end(argp);
   return result;
 }
+
+int Vsprintf (char* theBuffer, const char* theFormat, va_list theArgList)
+{
+  SAVE_TL();
+  return vsprintf_l(theBuffer, Standard_CLocaleSentry::GetCLocale(), theFormat, theArgList);
+}
index 319a63d22e9ad13550513b9a88fe7cb92a44de55..451fc2f5348621ad89c9bd9db7be02f132cb780f 100644 (file)
@@ -85,6 +85,15 @@ Standard_EXPORT int Fprintf (FILE* theFile, const char* theFormat, ...);
 //! Equivalent of standard C function sprintf() that always uses C locale
 Standard_EXPORT int Sprintf (char* theBuffer, const char* theFormat, ...);
 
+//! Equivalent of standard C function vsprintf() that always uses C locale.
+//! Note that this function does not check buffer bounds and should be used with precaution measures
+//! (only with format fitting into the buffer of known size).
+//! @param theBuffer  [in] [out] string buffer to fill
+//! @param theFormat  [in] format to apply
+//! @param theArgList [in] argument list for specified format
+//! @return the total number of characters written, or a negative number on error
+Standard_EXPORT int Vsprintf (char* theBuffer, const char* theFormat, va_list theArgList);
+
 #ifdef __cplusplus
 }
 #endif /* __cplusplus */
diff --git a/tests/bugs/mesh/bug28089_1 b/tests/bugs/mesh/bug28089_1
new file mode 100644 (file)
index 0000000..7e731c0
--- /dev/null
@@ -0,0 +1,24 @@
+puts "========"
+puts "0028089: Mesh - New algorithm for triangulation of 2d polygons"
+puts "========"
+puts ""
+
+pload MODELING VISUALIZATION
+restore [locate_data_file terrain.brep] t1
+tclean t1
+#tscale t1 0 0 0 0.01
+tcopy t1 t2
+ttranslate t1 0 0 0
+ttranslate t2 8000 0 0
+incmesh t1 1000 -algo 0 -a 50
+trinfo  t1
+incmesh t2 1000 -algo 1 -a 50
+trinfo  t2
+vclear
+vinit View1
+vtop
+vdefaults -autotriang 0
+vdisplay -dispMode 1 t1 t2
+vfit
+vaspects t1 -drawEdges 1
+vaspects t2 -drawEdges 1
diff --git a/tests/bugs/mesh/bug28089_2 b/tests/bugs/mesh/bug28089_2
new file mode 100644 (file)
index 0000000..02e3d13
--- /dev/null
@@ -0,0 +1,27 @@
+puts "========"
+puts "0028089: Mesh - New algorithm for triangulation of 2d polygons"
+puts "========"
+puts ""
+
+pload MODELING VISUALIZATION
+psphere s1 1
+psphere s2 1
+box b1 0 2 0 1 2 2
+box b2 2 2 0 1 2 2
+ttranslate s2 2 0 0
+incmesh s1 0.5 -algo 0 -a 50
+incmesh s2 0.5 -algo 1 -a 50
+incmesh b1 0.5 -algo 0
+incmesh b2 0.5 -algo 1
+
+vclear
+vinit View1
+vdefaults -autotriang 0
+vdisplay -dispMode 1 s1 s2 b1 b2
+vfit
+vaspects s1 -drawEdges 1
+vaspects s2 -drawEdges 1
+vaspects b1 -drawEdges 1
+vaspects b2 -drawEdges 1
+
+vdump ${imagedir}/${casename}.png