MeshFlattening: add function to unwrap face

CMakeLists.txt

@@ -587,6 +587,20 @@ endif(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
     if(FREECAD_USE_PCL)
         find_package(PCL REQUIRED COMPONENTS common kdtree features surface io filters segmentation sample_consensus)
     endif(FREECAD_USE_PCL)
+
+# -------------------------------- PyBind11 -----------------------------
+EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -c 
+  "import pybind11; print(bool(pybind11))"
+  OUTPUT_VARIABLE PYBIND11_FOUND OUTPUT_STRIP_TRAILING_WHITESPACE )
+
+EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -c 
+  "import pybind11; print(pybind11.get_include())"
+  OUTPUT_VARIABLE PYBIND11_INCLUDE_DIR OUTPUT_STRIP_TRAILING_WHITESPACE )
+
+if (PYBIND11_FOUND)
+    message(STATUS "successfully found pybind11")
+    message(STATUS "pybind11 include dir is: " ${PYBIND11_INCLUDE_DIR})
+endif()
 
 # -------------------------------- Boost --------------------------------
 

src/Mod/Mesh/Gui/Workbench.cpp

@@ -193,7 +193,7 @@ Gui::MenuItem* Workbench::setupMenuBar() const
           << "Mesh_Merge" << "Mesh_PolySelect" << "Mesh_PolyCut"
           << "Mesh_PolySplit" << "Mesh_PolySegm" << "Mesh_PolyTrim" << "Separator"
           << "Mesh_TrimByPlane" << "Mesh_SectionByPlane" << "Mesh_Segmentation"
-          << "Mesh_VertexCurvature" << "CreateFlatMesh";
+          << "Mesh_VertexCurvature" << "CreateFlatMesh" << "CreateFlatFace";
     return root;
 }
 

src/Mod/Mesh/InitGui.py

@@ -39,8 +39,12 @@ def __init__(self):
     def Initialize(self):
         import Mesh
         import MeshGui
-        import MeshFlatteningCommand
-
+        try:
+            import flatmesh
+            import MeshFlatteningCommand
+        except ImportError as e:
+            import FreeCAD
+            FreeCAD.Console.PrintLog((str(e)))
     def GetClassName(self):
         return "MeshGui::Workbench"
 

src/Mod/MeshPart/App/CMakeLists.txt

@@ -81,24 +81,25 @@ SET_PYTHON_PREFIX_SUFFIX(MeshPart)
 INSTALL(TARGETS MeshPart DESTINATION ${CMAKE_INSTALL_LIBDIR})
 
 
-
-################################ flat mesh ###############################
-SET(FLATMESH_SRCS
-    MeshFlatteningPy.cpp
-    MeshFlattening.cpp
-    MeshFlattening.h
-    MeshFlatteningNurbs.h
-    MeshFlatteningNurbs.cpp
-    MeshFlatteningLscmRelax.h
-    MeshFlatteningLscmRelax.cpp
-)
+if (PYBIND11_FOUND)
+    ################################ flat mesh ###############################
+    SET(FLATMESH_SRCS
+        MeshFlatteningPy.cpp
+        MeshFlattening.cpp
+        MeshFlattening.h
+        MeshFlatteningNurbs.h
+        MeshFlatteningNurbs.cpp
+        MeshFlatteningLscmRelax.h
+        MeshFlatteningLscmRelax.cpp
+    )
 
 
 
-add_library(flatmesh SHARED ${FLATMESH_SRCS})
-SET_PYTHON_PREFIX_SUFFIX(flatmesh)
-target_link_libraries(flatmesh ${PYTHON_LIBRARIES} ${MeshPart_LIBS})
+    add_library(flatmesh SHARED ${FLATMESH_SRCS})
+    SET_PYTHON_PREFIX_SUFFIX(flatmesh)
+    target_link_libraries(flatmesh ${PYTHON_LIBRARIES} ${MeshPart_LIBS})
 
-SET_BIN_DIR(flatmesh flatmesh /Mod/MeshPart)
-install(TARGETS flatmesh DESTINATION ${CMAKE_INSTALL_LIBDIR})
-############################################################################
+    SET_BIN_DIR(flatmesh flatmesh /Mod/MeshPart)
+    install(TARGETS flatmesh DESTINATION ${CMAKE_INSTALL_LIBDIR})
+    ############################################################################
+endif()

src/Mod/MeshPart/App/MeshFlattening.cpp

@@ -169,13 +169,13 @@ FaceUnwrapper::FaceUnwrapper(const TopoDS_Face& face)
     }
 }
 
-void FaceUnwrapper::findFlatNodes()
+void FaceUnwrapper::findFlatNodes(int steps, double val)
 {
     std::vector<long> fixed_pins;  //TODO: INPUT
     LscmRelax mesh_flattener(this->xyz_nodes.transpose(), this->tris.transpose(), fixed_pins);
     mesh_flattener.lscm();
-    for (int j=0; j<9; j++)
-        mesh_flattener.relax(0.9);
+    for (int j=0; j<steps; j++)
+        mesh_flattener.relax(0.95);
     this->ze_nodes = mesh_flattener.flat_vertices.transpose();
 }
 

src/Mod/MeshPart/App/MeshFlattening.h

@@ -62,7 +62,7 @@ class FaceUnwrapper{
 public:
 	FaceUnwrapper(const TopoDS_Face & face);
         FaceUnwrapper(ColMat<double, 3> xyz_nodes, ColMat<long, 3> tris);
-	void findFlatNodes();
+	void findFlatNodes(int steps, double val);
 	ColMat<double, 3> interpolateFlatFace(const TopoDS_Face& face);
         std::vector<ColMat<double, 3>> getFlatBoundaryNodes();
 

src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp

@@ -23,6 +23,7 @@
 
 #include "MeshFlatteningLscmRelax.h"
 #include <Eigen/IterativeLinearSolvers>
+#include<Eigen/SparseCholesky>
 #include <Eigen/SVD>
 #include <iostream>
 #include <algorithm>
@@ -278,10 +279,9 @@ void LscmRelax::relax(double weight)
     // rhs +=  K_g * Eigen::VectorXd::Ones(K_g.rows());
 
     // solve linear system (privately store the value for guess in next step)
-    Eigen::ConjugateGradient<spMat, Eigen::Lower> solver;
-    solver.setTolerance(0.0000001);
+    Eigen::SimplicialLDLT<spMat, Eigen::Lower> solver;
     solver.compute(K_g);
-    this->sol = solver.solveWithGuess(-rhs, this->sol);
+    this->sol = solver.solve(-rhs);
     this->set_shift(this->sol.head(this->vertices.cols() * 2) * weight);
     this->set_q_l_m();
 }

src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h

@@ -88,7 +88,7 @@ class LscmRelax{
     ColMat<double, 1> rhs;
     Eigen::MatrixXd MATRIX;
 
-    double nue=0.0;
+    double nue=0.9;
     double elasticity=1.;
 
     void lscm();

src/Mod/MeshPart/App/MeshFlatteningNurbs.cpp

@@ -25,6 +25,7 @@
 #include <iostream>
 #include "math.h"
 
+
 namespace nurbs{
 
 double divide(double a, double b)
@@ -35,6 +36,28 @@ double divide(double a, double b)
         return a / b;
 }
 
+Eigen::VectorXd NurbsBase1D::getKnotSequence(double u_min, double u_max, int u_deg, int num_u_poles)
+{
+    // boarder poles are on the surface
+    std::vector<double> u_knots;
+    for (int i=0; i < u_deg; i++)
+        u_knots.push_back(u_min);
+    for (int i=0; i < num_u_poles; i++)
+        u_knots.push_back(u_min + (u_max - u_min) * i / (num_u_poles - 1));
+    for (int i=0; i < u_deg; i++)
+        u_knots.push_back(u_max);
+    return Eigen::Map<Eigen::VectorXd>(u_knots.data(), u_knots.size());
+}
+
+Eigen::VectorXd NurbsBase1D::getWeightList(Eigen::VectorXd knots, int u_deg)
+{
+    Eigen::VectorXd weights;
+    weights.resize(knots.rows() - u_deg - 1);
+    weights.setOnes();
+    return weights;
+}
+
+
 // DE BOOR ALGORITHM FROM OPENGLIDER
 std::function<double(double)> get_basis(int degree, int i, Eigen::VectorXd knots)
     // Return a basis_function for the given degree """
@@ -180,7 +203,7 @@ void add_triplets(Eigen::VectorXd values, double row, std::vector<trip> &triplet
 spMat NurbsBase2D::getInfluenceMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U)
 {
     std::vector<trip> triplets;
-    for (int row_index; row_index < U.rows(); row_index++)
+    for (unsigned int row_index = 0; row_index < U.rows(); row_index++)
         add_triplets(this->getInfluenceVector(U.row(row_index)), row_index, triplets);
     spMat mat(U.rows(), this->u_functions.size() * this->v_functions.size());
     mat.setFromTriplets(triplets.begin(), triplets.end());
@@ -288,7 +311,7 @@ Eigen::VectorXd NurbsBase2D::getDvVector(Eigen::Vector2d u)
 spMat NurbsBase2D::getDuMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U)
 {
     std::vector<trip> triplets;
-    for (int row_index; row_index < U.rows(); row_index++)
+    for (unsigned int row_index = 0; row_index < U.rows(); row_index++)
         add_triplets(this->getDuVector(U.row(row_index)), row_index, triplets);
     spMat mat(U.rows(), this->u_functions.size() * this->v_functions.size());
     mat.setFromTriplets(triplets.begin(), triplets.end());
@@ -298,13 +321,66 @@ spMat NurbsBase2D::getDuMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U)
 spMat NurbsBase2D::getDvMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U)
 {
     std::vector<trip> triplets;
-    for (int row_index; row_index < U.rows(); row_index++)
+    for (unsigned int row_index = 0; row_index < U.rows(); row_index++)
         add_triplets(this->getDvVector(U.row(row_index)), row_index, triplets);
     spMat mat(U.rows(), this->u_functions.size() * this->v_functions.size());
     mat.setFromTriplets(triplets.begin(), triplets.end());
     return mat;
 }
 
+
+std::tuple<NurbsBase2D, Eigen::MatrixXd> NurbsBase2D::interpolateUBS(
+    Eigen::Matrix<double, Eigen::Dynamic, 3> poles,
+    int degree_u,
+    int degree_v,
+    int num_u_poles,
+    int num_v_poles,
+    int num_u_points,
+    int num_v_points)
+{
+    double u_min = this->u_knots(0);
+    double u_max = this->u_knots(this->u_knots.size() - 1);
+    double v_min = this->v_knots(0);
+    double v_max = this->v_knots(this->v_knots.size() - 1);
+    Eigen::VectorXd weights, u_knots, v_knots;
+    u_knots = NurbsBase1D::getKnotSequence(u_min, u_max, degree_u, num_u_poles);
+    v_knots = NurbsBase1D::getKnotSequence(v_min, v_max, degree_v, num_v_poles);
+
+    weights.resize((u_knots.rows() - degree_u - 1) * (v_knots.rows() - degree_v - 1));
+    weights.setOnes();
+    NurbsBase2D new_base(u_knots, v_knots, weights, degree_u, degree_v);
+    Eigen::Matrix<double, Eigen::Dynamic, 2> uv_points = this->getUVMesh(num_u_points, num_v_points);
+    Eigen::Matrix<double, Eigen::Dynamic, 3> xyz_points = this->getInfluenceMatrix(uv_points) * poles;
+    spMat A = new_base.getInfluenceMatrix(uv_points);
+    Eigen::LeastSquaresConjugateGradient<spMat > solver;
+    solver.compute(A);
+    Eigen::Matrix<double, Eigen::Dynamic, 3> new_poles = solver.solve(xyz_points);
+    return std::tuple<NurbsBase2D, Eigen::MatrixXd >(new_base, new_poles);
+}
+
+Eigen::Matrix<double, Eigen::Dynamic, 2> NurbsBase2D::getUVMesh(int num_u_points, int num_v_points)
+{
+    double u_min = this->u_knots(0);
+    double u_max = this->u_knots(this->u_knots.size() - 1);
+    double v_min = this->v_knots(0);
+    double v_max = this->v_knots(this->v_knots.size() - 1);
+    Eigen::Matrix<double, Eigen::Dynamic, 2> uv_points;
+    uv_points.resize(num_u_points * num_v_points, 2);
+    int i = 0;
+    for (int u = 0; u < num_u_points; u++)
+    {
+        for (int v = 0; v < num_v_points; v++)
+        {
+            uv_points(i, 0) = u_min + (u_max - u_min) * u / (num_u_points - 1);
+            uv_points(i, 1) = v_min + (v_max - v_min) * v / (num_v_points - 1);
+            i++;
+        }
+
+    }
+    return uv_points;
+}
+
+
 NurbsBase1D::NurbsBase1D(Eigen::VectorXd u_knots, Eigen::VectorXd weights, int degree_u)
 {
     this->u_knots = u_knots;
@@ -336,7 +412,7 @@ Eigen::VectorXd NurbsBase1D::getInfluenceVector(double u)
 spMat NurbsBase1D::getInfluenceMatrix(Eigen::VectorXd u)
 {
     std::vector<trip> triplets;
-    for (int row_index; row_index < u.size(); row_index++)
+    for (unsigned int row_index = 0; row_index < u.size(); row_index++)
         add_triplets(this->getInfluenceVector(u[row_index]), row_index, triplets);
     spMat mat(u.size(), this->u_functions.size());
     mat.setFromTriplets(triplets.begin(), triplets.end());
@@ -393,11 +469,44 @@ Eigen::VectorXd NurbsBase1D::getDuVector(double u)
 spMat NurbsBase1D::getDuMatrix(Eigen::VectorXd U)
 {
     std::vector<trip> triplets;
-    for (int row_index; row_index < U.size(); row_index++)
+    for (unsigned int row_index = 0; row_index < U.size(); row_index++)
         add_triplets(this->getDuVector(U[row_index]), row_index, triplets);
     spMat mat(U.size(), this->u_functions.size());
     mat.setFromTriplets(triplets.begin(), triplets.end());
     return mat;
 }
 
+std::tuple<NurbsBase1D, Eigen::Matrix<double, Eigen::Dynamic, 3>> NurbsBase1D::interpolateUBS(
+    Eigen::Matrix<double,
+    Eigen::Dynamic, 3> poles,
+    int degree,
+    int num_poles,
+    int num_points)
+{
+    double u_min = this->u_knots(0);
+    double u_max = this->u_knots(this->u_knots.size() - 1);
+    Eigen::VectorXd u_knots, weights;
+    u_knots = NurbsBase1D::getKnotSequence(u_min, u_max, degree, num_poles);
+    weights = NurbsBase1D::getWeightList(u_knots, degree);
+    NurbsBase1D new_base(u_knots, weights, degree);
+    Eigen::Matrix<double, Eigen::Dynamic, 1> u_points = this->getUMesh(num_points);
+    Eigen::Matrix<double, Eigen::Dynamic, 3> xyz_points;
+    xyz_points = this->getInfluenceMatrix(u_points) * poles;
+    spMat A = new_base.getInfluenceMatrix(u_points);
+    Eigen::LeastSquaresConjugateGradient<spMat > solver;
+    solver.compute(A);
+    Eigen::Matrix<double, Eigen::Dynamic, 3> new_poles = solver.solve(xyz_points);
+    return std::tuple<NurbsBase1D, Eigen::Matrix<double, Eigen::Dynamic, 3> >(new_base, new_poles);
+}
+
+Eigen::VectorXd NurbsBase1D::getUMesh(int num_u_points)
+{
+    double u_min = this->u_knots(0);
+    double u_max = this->u_knots(this->u_knots.size() - 1);
+    Eigen::Matrix<double, Eigen::Dynamic, 1> u_points;
+    u_points.setLinSpaced(num_u_points, u_min, u_max);
+    return u_points;    
+}
+
+
 }

src/Mod/MeshPart/App/MeshFlatteningNurbs.h

@@ -28,6 +28,7 @@
 #include <Eigen/Geometry>
 #include <Eigen/IterativeLinearSolvers>
 #include <Eigen/SparseCore>
+#include <tuple>
 
 namespace nurbs{
 
@@ -67,13 +68,24 @@ struct NurbsBase2D
 
     Eigen::VectorXd getDvVector(Eigen::Vector2d u);
     spMat getDvMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U);
+
+    Eigen::Matrix<double, Eigen::Dynamic, 2> getUVMesh(int num_u_points, int num_v_points);
+
+    std::tuple<NurbsBase2D, Eigen::MatrixXd> interpolateUBS(
+        Eigen::Matrix<double, Eigen::Dynamic, 3> poles,
+        int degree_u,
+        int degree_v,
+        int num_u_poles,
+        int num_v_poles,
+        int num_u_points,
+        int num_v_points);
 };
 
 struct NurbsBase1D
 {
     NurbsBase1D(){;};
     NurbsBase1D(Eigen::VectorXd u_knots, Eigen::VectorXd weights, int degree_u=3);
-    int degree_u = degree_u;
+    int degree_u;
     Eigen::VectorXd u_knots;
     Eigen::VectorXd weights;
     std::vector<std::function<double(double)>> u_functions;
@@ -88,6 +100,17 @@ struct NurbsBase1D
 
     Eigen::VectorXd getDuVector(double u);
     spMat getDuMatrix(Eigen::VectorXd u);
+
+    static Eigen::VectorXd getKnotSequence(double u_min, double u_max, int deg, int num_poles);
+    static Eigen::VectorXd getWeightList(Eigen::VectorXd knots, int u_deg);
+
+    Eigen::VectorXd getUMesh(int num_u_points);
+
+    std::tuple<NurbsBase1D, Eigen::Matrix<double, Eigen::Dynamic, 3>> interpolateUBS(
+        Eigen::Matrix<double, Eigen::Dynamic, 3> poles, 
+        int degree,
+        int num_u_poles,
+        int num_u_points);
 
 };