diff --git a/src/Mod/MeshPart/App/AppMeshPartPy.cpp b/src/Mod/MeshPart/App/AppMeshPartPy.cpp
index e09f6e2ae..e0b932400 100644
--- a/src/Mod/MeshPart/App/AppMeshPartPy.cpp
+++ b/src/Mod/MeshPart/App/AppMeshPartPy.cpp
@@ -146,6 +146,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
         if (PyArg_ParseTupleAndKeywords(args, kwds, "O!d", kwds_maxLength,
                                         &(Part::TopoShapePy::Type), &shape, &maxLength)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setMaxLength(maxLength);
             mesher.setRegular(true);
             return new Mesh::MeshPy(mesher.createMesh());
@@ -157,6 +158,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
         if (PyArg_ParseTupleAndKeywords(args, kwds, "O!d", kwds_maxArea,
                                         &(Part::TopoShapePy::Type), &shape, &maxArea)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setMaxArea(maxArea);
             mesher.setRegular(true);
             return new Mesh::MeshPy(mesher.createMesh());
@@ -168,6 +170,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
         if (PyArg_ParseTupleAndKeywords(args, kwds, "O!d", kwds_localLen,
                                         &(Part::TopoShapePy::Type), &shape, &localLen)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setLocalLength(localLen);
             mesher.setRegular(true);
             return new Mesh::MeshPy(mesher.createMesh());
@@ -179,6 +182,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
         if (PyArg_ParseTupleAndKeywords(args, kwds, "O!d", kwds_deflection,
                                         &(Part::TopoShapePy::Type), &shape, &deflection)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setDeflection(deflection);
             mesher.setRegular(true);
             return new Mesh::MeshPy(mesher.createMesh());
@@ -190,11 +194,13 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
         if (PyArg_ParseTupleAndKeywords(args, kwds, "O!dd", kwds_minmaxLen,
                                         &(Part::TopoShapePy::Type), &shape, &minLen, &maxLen)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setMinMaxLengths(minLen, maxLen);
             mesher.setRegular(true);
             return new Mesh::MeshPy(mesher.createMesh());
         }
 
+#if defined (HAVE_NETGEN)
         static char* kwds_fineness[] = {"Shape", "Fineness", "SecondOrder", "Optimize", "AllowQuad",NULL};
         PyErr_Clear();
         int fineness=0, secondOrder=0, optimize=1, allowquad=0;
@@ -202,6 +208,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
                                         &(Part::TopoShapePy::Type), &shape, &fineness,
                                         &secondOrder, &optimize, &allowquad)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Netgen);
             mesher.setFineness(fineness);
             mesher.setSecondOrder(secondOrder > 0);
             mesher.setOptimize(optimize > 0);
@@ -217,6 +224,7 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
                                         &growthRate, &nbSegPerEdge, &nbSegPerRadius,
                                         &secondOrder, &optimize, &allowquad)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+            mesher.setMethod(MeshPart::Mesher::Netgen);
             mesher.setGrowthRate(growthRate);
             mesher.setNbSegPerEdge(nbSegPerEdge);
             mesher.setNbSegPerRadius(nbSegPerRadius);
@@ -225,11 +233,17 @@ meshFromShape(PyObject *self, PyObject *args, PyObject* kwds)
             mesher.setQuadAllowed(allowquad > 0);
             return new Mesh::MeshPy(mesher.createMesh());
         }
+#endif
 
         PyErr_Clear();
         if (PyArg_ParseTuple(args, "O!", &(Part::TopoShapePy::Type), &shape)) {
             MeshPart::Mesher mesher(static_cast<Part::TopoShapePy*>(shape)->getTopoShapePtr()->_Shape);
+#if defined (HAVE_NETGEN)
+            mesher.setMethod(MeshPart::Mesher::Netgen);
+#else
+            mesher.setMethod(MeshPart::Mesher::Mefisto);
             mesher.setRegular(true);
+#endif
             return new Mesh::MeshPy(mesher.createMesh());
         }
     }
diff --git a/src/Mod/MeshPart/App/CMakeLists.txt b/src/Mod/MeshPart/App/CMakeLists.txt
index 08cdc91fe..94b35fc57 100644
--- a/src/Mod/MeshPart/App/CMakeLists.txt
+++ b/src/Mod/MeshPart/App/CMakeLists.txt
@@ -7,6 +7,10 @@ endif(MSVC)
 
 add_definitions(-DHAVE_SMESH)
 
+if(FREECAD_BUILD_FEM_NETGEN)
+    add_definitions(-DHAVE_NETGEN)
+endif(FREECAD_BUILD_FEM_NETGEN)
+
 include_directories(
     ${CMAKE_SOURCE_DIR}/src
     ${CMAKE_SOURCE_DIR}/src/3rdParty/salomesmesh/inc
diff --git a/src/Mod/MeshPart/App/Mesher.cpp b/src/Mod/MeshPart/App/Mesher.cpp
index 2ff6adda8..cdfb7684e 100644
--- a/src/Mod/MeshPart/App/Mesher.cpp
+++ b/src/Mod/MeshPart/App/Mesher.cpp
@@ -89,10 +89,24 @@ int MeshingOutput::sync()
 }
 
 Mesher::Mesher(const TopoDS_Shape& s)
-  : shape(s), maxLength(0), maxArea(0), localLength(0),
-    deflection(0), minLen(0), maxLen(0), regular(false),
-    fineness(5), growthRate(0), nbSegPerEdge(0), nbSegPerRadius(0),
-    secondOrder(false), optimize(true), allowquad(false)
+  : shape(s)
+  , method(None)
+  , maxLength(0)
+  , maxArea(0)
+  , localLength(0)
+  , deflection(0)
+  , minLen(0)
+  , maxLen(0)
+  , regular(false)
+#if defined (HAVE_NETGEN)
+  , fineness(5)
+  , growthRate(0)
+  , nbSegPerEdge(0)
+  , nbSegPerRadius(0)
+  , secondOrder(false)
+  , optimize(true)
+  , allowquad(false)
+#endif
 {
 }
 
@@ -111,77 +125,86 @@ Mesh::MeshObject* Mesher::createMesh() const
     SMESH_Mesh* mesh = meshgen->CreateMesh(0, true);
     int hyp=0;
 
-#if defined(_MSC_VER)
-    NETGENPlugin_Hypothesis_2D* hyp2d = new NETGENPlugin_Hypothesis_2D(hyp++,0,meshgen);
+    switch (method) {
+#if defined (HAVE_NETGEN)
+    case Netgen: {
+        NETGENPlugin_Hypothesis_2D* hyp2d = new NETGENPlugin_Hypothesis_2D(hyp++,0,meshgen);
 
-    if (fineness >=0 && fineness < 5) {
-        hyp2d->SetFineness(NETGENPlugin_Hypothesis_2D::Fineness(fineness));
-    }
-    // user defined values
-    else {
-        if (growthRate > 0)
-            hyp2d->SetGrowthRate(growthRate);
-        if (nbSegPerEdge > 0)
-            hyp2d->SetNbSegPerEdge(nbSegPerEdge);
-        if (nbSegPerRadius > 0)
-            hyp2d->SetNbSegPerRadius(nbSegPerRadius);
-    }
+        if (fineness >=0 && fineness < 5) {
+            hyp2d->SetFineness(NETGENPlugin_Hypothesis_2D::Fineness(fineness));
+        }
+        // user defined values
+        else {
+            if (growthRate > 0)
+                hyp2d->SetGrowthRate(growthRate);
+            if (nbSegPerEdge > 0)
+                hyp2d->SetNbSegPerEdge(nbSegPerEdge);
+            if (nbSegPerRadius > 0)
+                hyp2d->SetNbSegPerRadius(nbSegPerRadius);
+        }
 
-    hyp2d->SetQuadAllowed(allowquad);
-    hyp2d->SetOptimize(optimize);
-    hyp2d->SetSecondOrder(secondOrder); // apply bisecting to create four triangles out of one
-    hypoth.push_back(hyp2d);
-
-    NETGENPlugin_NETGEN_2D* alg2d = new NETGENPlugin_NETGEN_2D(hyp++,0,meshgen);
-    hypoth.push_back(alg2d);
-#else
-    if (maxLength > 0) {
-        StdMeshers_MaxLength* hyp1d = new StdMeshers_MaxLength(hyp++, 0, meshgen);
-        hyp1d->SetLength(maxLength);
-        hypoth.push_back(hyp1d);
-    }
-    else if (localLength > 0) {
-        StdMeshers_LocalLength* hyp1d = new StdMeshers_LocalLength(hyp++,0,meshgen);
-        hyp1d->SetLength(localLength);
-        hypoth.push_back(hyp1d);
-    }
-    else if (maxArea > 0) {
-        StdMeshers_MaxElementArea* hyp2d = new StdMeshers_MaxElementArea(hyp++,0,meshgen);
-        hyp2d->SetMaxArea(maxArea);
+        hyp2d->SetQuadAllowed(allowquad);
+        hyp2d->SetOptimize(optimize);
+        hyp2d->SetSecondOrder(secondOrder); // apply bisecting to create four triangles out of one
         hypoth.push_back(hyp2d);
-    }
-    else if (deflection > 0) {
-        StdMeshers_Deflection1D* hyp1d = new StdMeshers_Deflection1D(hyp++,0,meshgen);
-        hyp1d->SetDeflection(deflection);
-        hypoth.push_back(hyp1d);
-    }
-    else if (minLen > 0 && maxLen > 0) {
-        StdMeshers_Arithmetic1D* hyp1d = new StdMeshers_Arithmetic1D(hyp++,0,meshgen);
-        hyp1d->SetLength(minLen, false);
-        hyp1d->SetLength(maxLen, true);
-        hypoth.push_back(hyp1d);
-    }
-    else {
-        StdMeshers_AutomaticLength* hyp1d = new StdMeshers_AutomaticLength(hyp++,0,meshgen);
-        hypoth.push_back(hyp1d);
-    }
 
-    {
-        StdMeshers_NumberOfSegments* hyp1d = new StdMeshers_NumberOfSegments(hyp++,0,meshgen);
-        hyp1d->SetNumberOfSegments(1);
-        hypoth.push_back(hyp1d);
-    }
-
-    if (regular) {
-        StdMeshers_Regular_1D* hyp1d = new StdMeshers_Regular_1D(hyp++,0,meshgen);
-        hypoth.push_back(hyp1d);
-    }
-
-    StdMeshers_TrianglePreference* hyp2d_1 = new StdMeshers_TrianglePreference(hyp++,0,meshgen);
-    hypoth.push_back(hyp2d_1);
-    StdMeshers_MEFISTO_2D* alg2d = new StdMeshers_MEFISTO_2D(hyp++,0,meshgen);
-    hypoth.push_back(alg2d);
+        NETGENPlugin_NETGEN_2D* alg2d = new NETGENPlugin_NETGEN_2D(hyp++,0,meshgen);
+        hypoth.push_back(alg2d);
+    } break;
 #endif
+#if !defined (_MSC_VER)
+    case Mefisto: {
+        if (maxLength > 0) {
+            StdMeshers_MaxLength* hyp1d = new StdMeshers_MaxLength(hyp++, 0, meshgen);
+            hyp1d->SetLength(maxLength);
+            hypoth.push_back(hyp1d);
+        }
+        else if (localLength > 0) {
+            StdMeshers_LocalLength* hyp1d = new StdMeshers_LocalLength(hyp++,0,meshgen);
+            hyp1d->SetLength(localLength);
+            hypoth.push_back(hyp1d);
+        }
+        else if (maxArea > 0) {
+            StdMeshers_MaxElementArea* hyp2d = new StdMeshers_MaxElementArea(hyp++,0,meshgen);
+            hyp2d->SetMaxArea(maxArea);
+            hypoth.push_back(hyp2d);
+        }
+        else if (deflection > 0) {
+            StdMeshers_Deflection1D* hyp1d = new StdMeshers_Deflection1D(hyp++,0,meshgen);
+            hyp1d->SetDeflection(deflection);
+            hypoth.push_back(hyp1d);
+        }
+        else if (minLen > 0 && maxLen > 0) {
+            StdMeshers_Arithmetic1D* hyp1d = new StdMeshers_Arithmetic1D(hyp++,0,meshgen);
+            hyp1d->SetLength(minLen, false);
+            hyp1d->SetLength(maxLen, true);
+            hypoth.push_back(hyp1d);
+        }
+        else {
+            StdMeshers_AutomaticLength* hyp1d = new StdMeshers_AutomaticLength(hyp++,0,meshgen);
+            hypoth.push_back(hyp1d);
+        }
+
+        {
+            StdMeshers_NumberOfSegments* hyp1d = new StdMeshers_NumberOfSegments(hyp++,0,meshgen);
+            hyp1d->SetNumberOfSegments(1);
+            hypoth.push_back(hyp1d);
+        }
+
+        if (regular) {
+            StdMeshers_Regular_1D* hyp1d = new StdMeshers_Regular_1D(hyp++,0,meshgen);
+            hypoth.push_back(hyp1d);
+        }
+
+        StdMeshers_TrianglePreference* hyp2d_1 = new StdMeshers_TrianglePreference(hyp++,0,meshgen);
+        hypoth.push_back(hyp2d_1);
+        StdMeshers_MEFISTO_2D* alg2d = new StdMeshers_MEFISTO_2D(hyp++,0,meshgen);
+        hypoth.push_back(alg2d);
+    } break;
+#endif
+    default:
+        break;
+    }
 
     // Set new cout
     MeshingOutput stdcout;
@@ -272,7 +295,6 @@ Mesh::MeshObject* Mesher::createMesh() const
             faces.push_back(f3);
             faces.push_back(f4);
         }
-#if 0 // FIXME: how does the structure look like?
         else if (aFace->NbNodes() == 8) {
             MeshCore::MeshFacet f1, f2, f3, f4, f5, f6;
             const SMDS_MeshNode* node0 = aFace->GetNode(0);
@@ -281,8 +303,8 @@ Mesh::MeshObject* Mesher::createMesh() const
             const SMDS_MeshNode* node3 = aFace->GetNode(3);
             const SMDS_MeshNode* node4 = aFace->GetNode(4);
             const SMDS_MeshNode* node5 = aFace->GetNode(5);
-            const SMDS_MeshNode* node6 = aFace->GetNode(5);
-            const SMDS_MeshNode* node7 = aFace->GetNode(5);
+            const SMDS_MeshNode* node6 = aFace->GetNode(6);
+            const SMDS_MeshNode* node7 = aFace->GetNode(7);
 
             f1._aulPoints[0] = mapNodeIndex[node0];
             f1._aulPoints[1] = mapNodeIndex[node4];
@@ -300,13 +322,32 @@ Mesh::MeshObject* Mesher::createMesh() const
             f4._aulPoints[1] = mapNodeIndex[node7];
             f4._aulPoints[2] = mapNodeIndex[node6];
 
-            f5._aulPoints[0] = mapNodeIndex[node4];
-            f5._aulPoints[1] = mapNodeIndex[node6];
-            f5._aulPoints[2] = mapNodeIndex[node7];
+            // Two solutions are possible:
+            // <4,6,7>, <4,5,6> or <4,5,7>, <5,6,7>
+            Base::Vector3d v4(node4->X(),node4->Y(),node4->Z());
+            Base::Vector3d v5(node5->X(),node5->Y(),node5->Z());
+            Base::Vector3d v6(node6->X(),node6->Y(),node6->Z());
+            Base::Vector3d v7(node7->X(),node7->Y(),node7->Z());
+            double dist46 = Base::DistanceP2(v4,v6);
+            double dist57 = Base::DistanceP2(v5,v7);
+            if (dist46 > dist57) {
+                f5._aulPoints[0] = mapNodeIndex[node4];
+                f5._aulPoints[1] = mapNodeIndex[node6];
+                f5._aulPoints[2] = mapNodeIndex[node7];
 
-            f6._aulPoints[0] = mapNodeIndex[node4];
-            f6._aulPoints[1] = mapNodeIndex[node5];
-            f6._aulPoints[2] = mapNodeIndex[node6];
+                f6._aulPoints[0] = mapNodeIndex[node4];
+                f6._aulPoints[1] = mapNodeIndex[node5];
+                f6._aulPoints[2] = mapNodeIndex[node6];
+            }
+            else {
+                f5._aulPoints[0] = mapNodeIndex[node4];
+                f5._aulPoints[1] = mapNodeIndex[node5];
+                f5._aulPoints[2] = mapNodeIndex[node7];
+
+                f6._aulPoints[0] = mapNodeIndex[node5];
+                f6._aulPoints[1] = mapNodeIndex[node6];
+                f6._aulPoints[2] = mapNodeIndex[node7];
+            }
 
             faces.push_back(f1);
             faces.push_back(f2);
@@ -315,7 +356,6 @@ Mesh::MeshObject* Mesher::createMesh() const
             faces.push_back(f5);
             faces.push_back(f6);
         }
-#endif
         else {
             Base::Console().Warning("Face with %d nodes ignored\n", aFace->NbNodes());
         }
diff --git a/src/Mod/MeshPart/App/Mesher.h b/src/Mod/MeshPart/App/Mesher.h
index f3ac74f2b..5b1a8d119 100644
--- a/src/Mod/MeshPart/App/Mesher.h
+++ b/src/Mod/MeshPart/App/Mesher.h
@@ -34,9 +34,22 @@ namespace MeshPart {
 class Mesher
 {
 public:
+    enum Method {
+        None = 0,
+        Mefisto = 1,
+#if defined (HAVE_NETGEN)
+        Netgen = 2,
+#endif
+    };
+
     Mesher(const TopoDS_Shape&);
     ~Mesher();
 
+    void setMethod(Method m)
+    { method = m; }
+    Method getMethod() const
+    { return method; }
+
     /** @name Mefisto settings */
     //@{
     void setMaxLength(double s)
@@ -65,6 +78,7 @@ public:
     { return regular; }
     //@}
 
+#if defined (HAVE_NETGEN)
     /** @name Netgen settings */
     //@{
     void setFineness(int s)
@@ -96,17 +110,20 @@ public:
     bool isQuadAllowed() const
     { return allowquad; }
     //@}
+#endif
 
     Mesh::MeshObject* createMesh() const;
 
 private:
     const TopoDS_Shape& shape;
+    Method method;
     double maxLength;
     double maxArea;
     double localLength;
     double deflection;
     double minLen, maxLen;
     bool regular;
+#if defined (HAVE_NETGEN)
     int fineness;
     double growthRate;
     double nbSegPerEdge;
@@ -114,6 +131,7 @@ private:
     bool secondOrder;
     bool optimize;
     bool allowquad;
+#endif
 };
 
 class MeshingOutput : public std::streambuf