Script to create Menger sponge

src/Mod/Mesh/App/Core/Degeneration.cpp

@@ -380,6 +380,30 @@ bool MeshFixDuplicateFacets::Fixup()
 
 // ----------------------------------------------------------------------
 
+bool MeshEvalInternalFacets::Evaluate()
+{
+    _indices.clear();
+    unsigned long uIndex=0;
+    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
+
+    // get all facets
+    std::set<FaceIterator, MeshFacet_Less > aFaceSet;
+    MeshFacetArray::_TConstIterator first = rFaces.begin();
+    for (MeshFacetArray::_TConstIterator it = rFaces.begin(); it != rFaces.end(); ++it, uIndex++) {
+        std::pair<std::set<FaceIterator, MeshFacet_Less>::iterator, bool>
+        pI = aFaceSet.insert(it);
+        if (!pI.second) {
+            // collect both elements
+            _indices.push_back(*pI.first - first);
+            _indices.push_back(uIndex);
+        }
+    }
+
+    return _indices.empty();
+}
+
+// ----------------------------------------------------------------------
+
 bool MeshEvalDegeneratedFacets::Evaluate()
 {
   MeshFacetIterator it(_rclMesh);

src/Mod/Mesh/App/Core/Degeneration.h

@@ -188,6 +188,35 @@ public:
   bool Fixup ();
 };
 
+/**
+ * The MeshEvalInternalFacets class identifies internal facets of a volume mesh.
+ * @author Werner Mayer
+ */
+class MeshExport MeshEvalInternalFacets : public MeshEvaluation
+{
+public:
+  /**
+   * Construction.
+   */
+  MeshEvalInternalFacets (const MeshKernel &rclM) : MeshEvaluation( rclM ) { }
+  /** 
+   * Destruction.
+   */
+  ~MeshEvalInternalFacets () { }
+  /** 
+   * Identifiy internal facets.
+   */
+  bool Evaluate ();
+  /**
+   * Return the indices.
+   */
+  const std::vector<unsigned long>& GetIndices() const
+  { return _indices; }
+
+private:
+    std::vector<unsigned long> _indices;
+};
+
 /**
  * The MeshEvalDegeneratedFacets class searches for degenerated facets. A facet is degenerated either if its points
  * are collinear, i.e. they lie on a line or two points are coincident. In the latter case these points are duplicated.

src/Mod/Mesh/App/MeshPy.xml

@@ -127,6 +127,16 @@ Example:
 				<UserDocu>Remove a list of facet indices from the mesh</UserDocu>
 			</Documentation>
 		</Methode>
+		<Methode Name="getInternalFacets" Const="true">
+			<Documentation>
+				<UserDocu>Builds a list of facet indices with triangles that are inside a volume mesh</UserDocu>
+			</Documentation>
+		</Methode>
+		<Methode Name="rebuildNeighbourHood">
+			<Documentation>
+				<UserDocu>Repairs the neighbourhood which might be broken</UserDocu>
+			</Documentation>
+		</Methode>
 		<Methode Name="addMesh">
 			<Documentation>
 				<UserDocu>Combine this mesh with another mesh.</UserDocu>

src/Mod/Mesh/App/MeshPyImp.cpp

@@ -36,6 +36,7 @@
 #include "MeshProperties.h"
 #include "Core/Algorithm.h"
 #include "Core/Iterator.h"
+#include "Core/Degeneration.h"
 #include "Core/Elements.h"
 #include "Core/Grid.h"
 #include "Core/MeshKernel.h"
@@ -581,6 +582,35 @@ PyObject* MeshPy::removeFacets(PyObject *args)
     Py_Return;
 }
 
+PyObject* MeshPy::getInternalFacets(PyObject *args)
+{
+    if (!PyArg_ParseTuple(args, ""))
+        return 0;
+
+    const MeshCore::MeshKernel& kernel = getMeshObjectPtr()->getKernel();
+    MeshCore::MeshEvalInternalFacets eval(kernel);
+    eval.Evaluate();
+
+    const std::vector<unsigned long>& indices = eval.GetIndices();
+    Py::List ary(indices.size());
+    Py::List::size_type pos=0;
+    for (std::vector<unsigned long>::const_iterator it = indices.begin(); it != indices.end(); ++it) {
+        ary[pos++] = Py::Long(*it);
+    }
+
+    return Py::new_reference_to(ary);
+}
+
+PyObject* MeshPy::rebuildNeighbourHood(PyObject *args)
+{
+    if (!PyArg_ParseTuple(args, ""))
+        return 0;
+
+    MeshCore::MeshKernel& kernel = getMeshObjectPtr()->getKernel();
+    kernel.RebuildNeighbours();
+    Py_Return;
+}
+
 PyObject*  MeshPy::addMesh(PyObject *args)
 {
     PyObject* mesh;

src/Mod/TemplatePyMod/MengerSponge.py

@@ -0,0 +1,120 @@
+# Script to create a Menger sponge
+# (c) 2012 Werner Mayer LGPL
+
+# The script is based on the work of daxmick at
+# https://sourceforge.net/apps/phpbb/free-cad/viewtopic.php?f=3&t=2307
+
+import threading
+import Mesh, MeshGui
+from FreeCAD import Base
+
+# Create a Box and Place it a coords (x,y,z)
+def PlaceBox(x,y,z):
+   box = Mesh.createBox(1,1,1)
+   box.translate(x,y,z)
+   return box
+
+def Sierpinski(level,x0,y0,z0):
+   #print threading.current_thread().name
+   boxnums = pow(3,level)
+   thirds = boxnums / 3
+   twothirds = thirds * 2
+   if(level == 0):
+      rangerx = [x0]
+      rangery = [y0]
+      rangerz = [z0]
+   else:
+      rangerx = [ x0, x0 + thirds, x0 + twothirds ]
+      rangery = [ y0, y0 + thirds, y0 + twothirds ]
+      rangerz = [ z0, z0 + thirds, z0 + twothirds ]
+   block = 1
+   skip=[5,11,13,14,15,17,23]
+   mesh=Mesh.Mesh()
+   for i in rangerx:
+      for j in rangery:
+         for k in rangerz:
+            if block not in skip:
+               if(level > 0):
+                  mesh.addMesh(Sierpinski(level-1,i,j,k))
+               else:
+                  mesh.addMesh(PlaceBox(i,j,k))
+            block+=1
+   return mesh
+
+### Multi-threaded ###
+
+class MengerThread(threading.Thread):
+    def __init__(self,args):
+        self.args=args
+        self.mesh=Mesh.Mesh()
+        threading.Thread.__init__(self)
+    def run(self):
+        for i in self.args:
+            self.mesh.addMesh(Sierpinski(*i))
+
+def makeMengerSponge_mt(level=3,x0=0,y0=0,z0=0):
+    """
+    Is much slower than makeMengerSponge!!! :(
+    """
+    if level == 0:
+        mesh=Sierpinski(level,x0,y0,z0)
+        Mesh.show(mesh)
+        return
+
+    boxnums = pow(3,level)
+    thirds = boxnums / 3
+    twothirds = thirds * 2
+    
+    rangerx = [ x0, x0 + thirds, x0 + twothirds ]
+    rangery = [ y0, y0 + thirds, y0 + twothirds ]
+    rangerz = [ z0, z0 + thirds, z0 + twothirds ]
+    block = 1
+    skip=[5,11,13,14,15,17,23]
+
+    # collect the arguments for the algorithm in a list
+    args=[]
+    for i in rangerx:
+        for j in rangery:
+            for k in rangerz:
+                if block not in skip:
+                    args.append((level-1,i,j,k))
+                block+=1
+
+    numJobs = 4
+    threads=[]
+    while numJobs > 0:
+        size = len(args)
+        count = size / numJobs
+        numJobs-=1
+        thr=MengerThread(args[:count])
+        threads.append(thr)
+        args=args[count:]
+
+    print "Number of threads: %i" % (len(threads))
+    for thr in threads:
+        thr.start()
+    for thr in threads:
+        thr.join()
+
+    mesh=Mesh.Mesh()
+    for thr in threads:
+        mesh.addMesh(thr.mesh)
+        del thr.mesh
+        
+    print mesh
+    mesh.removeDuplicatedPoints()
+    mesh.removeFacets(mesh.getInternalFacets())
+    mesh.rebuildNeighbourHood()
+    print "Mesh is solid: %s" % (mesh.isSolid())
+    Mesh.show(mesh)
+
+
+### Single-threaded ###
+
+def makeMengerSponge(level=3,x0=0,y0=0,z0=0):
+    mesh=Sierpinski(level,x0,y0,z0)
+    mesh.removeDuplicatedPoints()
+    mesh.removeFacets(mesh.getInternalFacets())
+    mesh.rebuildNeighbourHood()
+    print "Mesh is solid: %s" % (mesh.isSolid())
+    Mesh.show(mesh)