diff --git a/src/Mod/MeshPart/App/MeshFlattening.cpp b/src/Mod/MeshPart/App/MeshFlattening.cpp index 2b0cc22792..a185d2b16b 100644 --- a/src/Mod/MeshPart/App/MeshFlattening.cpp +++ b/src/Mod/MeshPart/App/MeshFlattening.cpp @@ -176,7 +176,7 @@ void FaceUnwrapper::findFlatNodes(int steps, double val) lscmrelax::LscmRelax mesh_flattener(this->xyz_nodes.transpose(), this->tris.transpose(), fixed_pins); mesh_flattener.lscm(); for (int j=0; jze_nodes = mesh_flattener.flat_vertices.transpose(); } diff --git a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp index 8e1225c16b..bcce25e4ba 100644 --- a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp +++ b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp @@ -295,11 +295,13 @@ void LscmRelax::relax(double weight) void LscmRelax::area_relax(double weight) { // TODO: doesn't work so far + if (this->sol.size() == 0) + this->sol.Zero(this->vertices.cols()); std::vector K_g_triplets; - spMat K_g(this->vertices.cols() * 2 + 3, this->vertices.cols() * 2 + 3); - spMat K_g_lsq(this->triangles.cols(), this->vertices.cols() * 2 + 3); + spMat K_g(this->vertices.cols() * 2, this->vertices.cols() * 2); + spMat K_g_lsq(this->triangles.cols(), this->vertices.cols() * 2); Eigen::VectorXd rhs_lsq(this->triangles.cols()); - Eigen::VectorXd rhs(this->vertices.cols() * 2 + 3); + Eigen::VectorXd rhs(this->vertices.cols() * 2); rhs.setZero(); Eigen::Matrix B; @@ -337,29 +339,11 @@ void LscmRelax::area_relax(double weight) K_g_lsq.setFromTriplets(K_g_triplets.begin(), K_g_triplets.end()); K_g_triplets.clear(); - - for (long i=0; i < this->flat_vertices.cols() ; i++) - { - // fixing total ux - K_g_triplets.push_back(trip(i * 2, this->flat_vertices.cols() * 2, 1)); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2, i * 2, 1)); - // fixing total uy - K_g_triplets.push_back(trip(i * 2 + 1, this->flat_vertices.cols() * 2 + 1, 1)); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 1, i * 2 + 1, 1)); - // fixing ux*y-uy*x - K_g_triplets.push_back(trip(i * 2, this->flat_vertices.cols() * 2 + 2, - this->flat_vertices(1, i))); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 2, i * 2, - this->flat_vertices(1, i))); - K_g_triplets.push_back(trip(i * 2 + 1, this->flat_vertices.cols() * 2 + 2, this->flat_vertices(0, i))); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 2, i * 2 + 1, this->flat_vertices(0, i))); - } K_g.setFromTriplets(K_g_triplets.begin(), K_g_triplets.end()); - K_g += K_g_lsq.transpose() * K_g_lsq; - rhs = K_g_lsq.transpose() * rhs_lsq; - Eigen::SimplicialLDLT solver; + Eigen::ConjugateGradient solver; solver.compute(K_g); this->sol = solver.solve(-rhs); - this->set_shift(this->sol.head(this->vertices.cols() * 2) * weight); - this->set_q_l_m(); + this->set_shift(this->sol * weight); } @@ -383,10 +367,13 @@ void LscmRelax::edge_relax(double weight) } } // 2. create system - + + if (this->sol.size() == 0) + this->sol.Zero(this->vertices.cols()); + std::vector K_g_triplets; - spMat K_g(this->vertices.cols() * 2 + 3, this->vertices.cols() * 2 + 3); - Eigen::VectorXd rhs(this->vertices.cols() * 2 + 3); + spMat K_g(this->vertices.cols() * 2, this->vertices.cols() * 2); + Eigen::VectorXd rhs(this->vertices.cols() * 2); rhs.setZero(); for(auto edge : edges) { @@ -420,25 +407,13 @@ void LscmRelax::edge_relax(double weight) rhs(indices[row]) += rhs_m[row]; } } - for (long i=0; i < this->flat_vertices.cols() ; i++) - { - // fixing total ux - K_g_triplets.push_back(trip(i * 2, this->flat_vertices.cols() * 2, 1)); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2, i * 2, 1)); - // fixing total uy - K_g_triplets.push_back(trip(i * 2 + 1, this->flat_vertices.cols() * 2 + 1, 1)); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 1, i * 2 + 1, 1)); - // fixing ux*y-uy*x - K_g_triplets.push_back(trip(i * 2, this->flat_vertices.cols() * 2 + 2, - this->flat_vertices(1, i))); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 2, i * 2, - this->flat_vertices(1, i))); - K_g_triplets.push_back(trip(i * 2 + 1, this->flat_vertices.cols() * 2 + 2, this->flat_vertices(0, i))); - K_g_triplets.push_back(trip(this->flat_vertices.cols() * 2 + 2, i * 2 + 1, this->flat_vertices(0, i))); - } + K_g.setFromTriplets(K_g_triplets.begin(), K_g_triplets.end()); - Eigen::SimplicialLDLT solver; + Eigen::ConjugateGradient solver; + solver.preconditioner().setNullSpace(this->get_nullspace()); solver.compute(K_g); this->sol = solver.solve(-rhs); - this->set_shift(this->sol.head(this->vertices.cols() * 2) * weight); + this->set_shift(this->sol * weight); } @@ -721,4 +696,21 @@ std::vector LscmRelax::get_fem_fixed_pins() return std::vector{min_x_index * 2, min_x_index * 2 + 1, max_x_index * 2 + 1}; } +Eigen::MatrixXd LscmRelax::get_nullspace() +{ + Eigen::MatrixXd null_space; + null_space.setZero(this->flat_vertices.size() * 2, 3); + + for (int i=0; iflat_vertices.cols(); i++) + { + null_space(i * 2, 0) = 1; // x-translation + null_space(i * 2 + 1, 1) = 1; // y-translation + null_space(i * 2, 2) = - this->flat_vertices(1, i); // z-rotation + null_space(i * 2 + 1, 2) = this->flat_vertices(0, i); // z-rotation + + } + return null_space; +} + + } diff --git a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h index 87b7c31735..fee5c04727 100644 --- a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h +++ b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h @@ -48,6 +48,24 @@ typedef Eigen::SparseMatrix spMat; namespace lscmrelax { + +class NullSpaceProjector: public Eigen::IdentityPreconditioner +{ + public: + Eigen::MatrixXd null_space_1; + Eigen::MatrixXd null_space_2; + + template + inline Rhs solve(Rhs& b) const { + return b - this->null_space_1 * (this->null_space_2 * b); + } + + void setNullSpace(Eigen::MatrixXd null_space) { + // normalize + orthogonalize the nullspace + this->null_space_1 = null_space * ((null_space.transpose() * null_space).inverse()); + this->null_space_2 = null_space.transpose(); + } +}; typedef Eigen::Vector3d Vector3; typedef Eigen::Vector2d Vector2; @@ -70,6 +88,7 @@ private: Eigen::VectorXd sol; std::vector get_fem_fixed_pins(); + Eigen::MatrixXd get_nullspace(); public: LscmRelax( diff --git a/src/Mod/MeshPart/Gui/MeshFlatteningCommand.py b/src/Mod/MeshPart/Gui/MeshFlatteningCommand.py index c96381f6d9..7283d2edd2 100644 --- a/src/Mod/MeshPart/Gui/MeshFlatteningCommand.py +++ b/src/Mod/MeshPart/Gui/MeshFlatteningCommand.py @@ -30,7 +30,7 @@ class CreateFlatMesh(BaseCommand): points = np.array([[i.x, i.y, i.z] for i in obj.Mesh.Points]) faces = np.array([list(i) for i in obj.Mesh.Topology[1]]) flattener = flatmesh.FaceUnwrapper(points, faces) - flattener.findFlatNodes(5, 0.99) + flattener.findFlatNodes(5, 0.95) boundaries = flattener.getFlatBoundaryNodes() print('number of nodes: {}'.format(len(flattener.ze_nodes))) print('number of faces: {}'.format(len(flattener.tris)))