diff --git a/src/Mod/MeshPart/App/CMakeLists.txt b/src/Mod/MeshPart/App/CMakeLists.txt index d3efec98f8..212a41bbd4 100644 --- a/src/Mod/MeshPart/App/CMakeLists.txt +++ b/src/Mod/MeshPart/App/CMakeLists.txt @@ -21,6 +21,7 @@ include_directories( ${XercesC_INCLUDE_DIRS} ${SMESH_INCLUDE_DIR} ${VTK_INCLUDE_DIRS} + ${EIGEN3_INCLUDE_DIR} ) @@ -78,3 +79,21 @@ SET_BIN_DIR(MeshPart MeshPart /Mod/MeshPart) SET_PYTHON_PREFIX_SUFFIX(MeshPart) INSTALL(TARGETS MeshPart DESTINATION ${CMAKE_INSTALL_LIBDIR}) + +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}) + +SET_BIN_DIR(flatmesh flatmesh /Mod/MeshPart) +install(TARGETS flatmesh DESTINATION ${CMAKE_INSTALL_LIBDIR}) diff --git a/src/Mod/MeshPart/App/MeshFlattening.cpp b/src/Mod/MeshPart/App/MeshFlattening.cpp new file mode 100644 index 0000000000..5458d6b3f8 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlattening.cpp @@ -0,0 +1,142 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#include "MeshFlattening.h" +#include "MeshFlatteningLscmRelax.h" +#include +#include +#include +#include + +FaceUnwrapper::FaceUnwrapper(const TopoDS_Face& face) +{ + int i = 0; +// transform to nurbs: + TopLoc_Location location; + +// triangulate: + const Handle(Poly_Triangulation) &triangulation = BRep_Tool::Triangulation(face, location); + + if (triangulation.IsNull()) + throw std::runtime_error("null triangulation in face construction"); + +// compute uv coordinates + if (triangulation->HasUVNodes()) + { + const TColgp_Array1OfPnt2d &_uv_nodes = triangulation->UVNodes(); + this->uv_nodes.resize(triangulation->NbNodes(), 2); + i = 0; + for (gp_Pnt2d _uv_node: _uv_nodes) + { + this->uv_nodes.row(i) << _uv_node.X(), _uv_node.Y(); + i++; + } + } +// + const TColgp_Array1OfPnt &_nodes = triangulation->Nodes(); + this->xyz_nodes.resize(triangulation->NbNodes(), 3); + i = 0; + for (gp_Pnt _node: _nodes) + { + this->xyz_nodes.row(i) << _node.X(), _node.Y(), _node.Z(); + i++; + } + + const Poly_Array1OfTriangle &_tris = triangulation->Triangles(); + this->tris.resize(triangulation->NbTriangles(), 3); + i = 0; + for (Poly_Triangle _tri: _tris) + { + int n1, n2, n3; + _tri.Get(n1, n2, n3); + this->tris.row(i) << n1-1, n2-1, n3-1; + i++; + } +} + +void FaceUnwrapper::findFlatNodes() +{ + std::vector fixed_pins; + LscmRelax mesh_flattener(this->xyz_nodes.transpose(), this->tris.transpose(), fixed_pins); + mesh_flattener.lscm(); + mesh_flattener.relax(0.9); + this->ze_nodes = mesh_flattener.flat_vertices.transpose(); +} + +ColMat FaceUnwrapper::interpolateNurbsFace(const TopoDS_Face& face) +{ + // extract xyz poles, knots, weights, degree + const Handle(Geom_Surface) &_surface = BRep_Tool::Surface(face); + const Handle(Geom_BSplineSurface) &_bspline = Handle(Geom_BSplineSurface)::DownCast(_surface); + const TColStd_Array1OfReal &_uknots = _bspline->UKnotSequence(); + const TColStd_Array1OfReal &_vknots = _bspline->VKnotSequence(); + + Eigen::VectorXd weights; + weights.resize(_bspline->NbUPoles() * _bspline->NbVPoles()); + int i = 0; + for (int u=1; u <= _bspline->NbUPoles(); u++) + { + for (int v=1; v <= _bspline->NbVPoles(); v++) + { + weights[i] = _bspline->Weight(u, v); + i++; + } + } + + Eigen::VectorXd u_knots; + Eigen::VectorXd v_knots; + u_knots.resize(_uknots.Size()); + v_knots.resize(_vknots.Size()); + for (int u=1; u <= _uknots.Size(); u++) + { + u_knots[u - 1] = _uknots.Value(u); + } + for (int v=1; v <= _vknots.Size(); v++) + { + v_knots[v - 1] = _vknots.Value(v); + } + + + nu = nurbs::NurbsBase2D(u_knots, v_knots, weights, _bspline->UDegree(), _bspline->VDegree()); + A = nu.getInfluenceMatrix(this->uv_nodes); + + Eigen::LeastSquaresConjugateGradient solver; + solver.compute(A); + ColMat ze_poles; + ColMat flat_poles; + ze_poles.resize(weights.rows(), 2); + flat_poles.resize(weights.rows(), 3); + flat_poles.setZero(); + ze_poles = solver.solve(ze_nodes); + flat_poles.col(0) << ze_poles.col(0); + flat_poles.col(1) << ze_poles.col(1); + return flat_poles; +} + + + + + + + + diff --git a/src/Mod/MeshPart/App/MeshFlattening.h b/src/Mod/MeshPart/App/MeshFlattening.h new file mode 100644 index 0000000000..430190f669 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlattening.h @@ -0,0 +1,73 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + +#ifndef MESHFLATTENING +#define MESHFLATTENING + +#include "MeshFlatteningNurbs.h" +#include +#include +#include + +#include +#include +#include +#include + + +typedef Eigen::Vector3d Vector3; +typedef Eigen::Vector2d Vector2; + +template +using ColMat = Eigen::Matrix; + +template +using RowMat = Eigen::Matrix; + + +typedef Eigen::Triplet trip; +typedef Eigen::SparseMatrix spMat; + +class FaceUnwrapper{ + nurbs::NurbsBase2D nu; +public: + FaceUnwrapper(const TopoDS_Face & face); + void findFlatNodes(); + ColMat interpolateNurbsFace(const TopoDS_Face& face); + + bool use_nurbs = true; + // the mesh + ColMat tris; // input + ColMat fixed_nodes; // input + ColMat xyz_nodes; // compute from uv_mesh (xyz = A * poles) + ColMat uv_nodes; // input + ColMat ze_nodes; // copute + + // nurbs + ColMat ze_poles; // compute + spMat A; // mapping between nurbs(poles) and mesh(vertices) computed with nurbs-basis-functions and uv_mesh + + // + +}; + +#endif // MESHFLATTENING diff --git a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp new file mode 100644 index 0000000000..471d720b70 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.cpp @@ -0,0 +1,567 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#include "MeshFlatteningLscmRelax.h" +#include +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846f +#endif + +// TODO: +// vertices, flat vertices = EIgen::MatrixXd (better python conversation, parallell) +// make it a own library +// area constrained (scale the final unwrapped mesh to the original area) +// FEM approach + + +typedef Eigen::Triplet trip; +typedef Eigen::SparseMatrix spMat; + + + +ColMat map_to_2D(ColMat points) +{ + ColMat mat(points.size(), 4); + mat << points, ColMat::Ones(points.size()); + Eigen::JacobiSVD svd(mat, Eigen::ComputeThinU | Eigen::ComputeThinV); + Vector3 n = svd.matrixV().row(2); + + n.normalize(); + Vector3 y(0, 1, 0); + if ((n - y).norm() < 0.0001) + y = Vector3(1, 0, 0); + Vector3 x = n.cross(y); + x.normalize(); + y = n.cross(x); + Eigen::Matrix transform; + transform.col(0) = x; + transform.col(1) = y; + return points * transform; +} + + +unsigned int get_max_distance(Vector3 point, RowMat vertices, double & max_dist) +{ + max_dist = 0; + double dist; + unsigned long max_dist_index; + long j = 0; + for (j=0; j < vertices.cols(); j++) + { + // debugging + dist = (point - vertices.col(j)).norm(); + if (dist > max_dist) + { + max_dist = dist; + max_dist_index = j; + } + } + return max_dist_index; +} + + +LscmRelax::LscmRelax( + RowMat vertices, + RowMat triangles, + std::vector fixed_pins) +{ + this->vertices = vertices; + this->triangles = triangles; + this->flat_vertices.resize(2, this->vertices.cols()); + this->fixed_pins = fixed_pins; + + // set the fixed pins of the flat-mesh: + this->set_fixed_pins(); + + + unsigned int fixed_count = 0; + for (long i=0; i < this->vertices.cols(); i++) + { + if (fixed_count < this->fixed_pins.size()) + { + if (i == this->fixed_pins[fixed_count]) + fixed_count ++; + else + this->old_order.push_back(i); + } + else + this->old_order.push_back(i); + } + + for (auto fixed_index: this->fixed_pins) + this->old_order.push_back(fixed_index); + + // get the reversed map: + this->new_order.resize(this->old_order.size()); + long j = 0; + for (auto index: this->old_order) + { + this->new_order[index] = j; + j++; + } + + // set the C-mat + this->C << 1, nue, 0, + nue, 1, 0, + 0, 0, (1 - nue) / 2; + this->C *= elasticity / (1 - nue * nue); +} + + + +////////////////////////////////////////////////////////////////////////// +///////////////// F.E.M ///////////// +////////////////////////////////////////////////////////////////////////// +void LscmRelax::relax(double weight) +{ + ColMat d_q_l_g = this->q_l_m - this->q_l_g; + // for every triangle + Eigen::Matrix B; + Eigen::Matrix T; + Eigen::Matrix K_m; + Eigen::Matrix u_m, rhs_m; + Eigen::VectorXd rhs(this->vertices.cols() * 2 + 3); + if (this->sol.size() == 0) + this->sol.Zero(this->vertices.cols() * 2 + 3); + spMat K_g(this->vertices.cols() * 2 + 3, this->vertices.cols() * 2 + 3); + std::vector K_g_triplets; + Vector2 v1, v2, v3, v12, v23, v31; + long row_pos, col_pos; + double A; + + rhs.setZero(); + + for (long i=0; itriangles.cols(); i++) + { + // 1: construct B-mat in m-system + v1 = this->flat_vertices.col(this->triangles(0, i)); + v2 = this->flat_vertices.col(this->triangles(1, i)); + v3 = this->flat_vertices.col(this->triangles(2, i)); + v12 = v2 - v1; + v23 = v3 - v2; + v31 = v1 - v3; + B << -v23.y(), 0, -v31.y(), 0, -v12.y(), 0, + 0, v23.x(), 0, v31.x(), 0, v12.x(), + -v23.x(), v23.y(), -v31.x(), v31.y(), -v12.x(), v12.y(); + T << v12.x(), -v12.y(), + v12.y(), v12.x(); + T /= v12.norm(); + A = std::abs(this->q_l_m(i, 0) * this->q_l_m(i, 2) / 2); + B /= A * 2; // (2*area) + + // 2: sigma due dqlg in m-system + u_m << Vector2(0, 0), T * Vector2(d_q_l_g(i, 0), 0), T * Vector2(d_q_l_g(i, 1), d_q_l_g(i, 2)); + + // 3: rhs_m = B.T * C * B * dqlg_m + // K_m = B.T * C * B + rhs_m = B.transpose() * this->C * B * u_m * A; + K_m = B.transpose() * this->C * B * A; + + // 5: add to rhs_g, K_g + for (int j=0; j < 3; j++) + { + row_pos = this->triangles(j, i); + rhs[row_pos * 2] += rhs_m[j * 2]; + rhs[row_pos * 2 + 1] += rhs_m[j * 2 +1]; + for (int k=0; k < 3; k++) + { + col_pos = this->triangles(k, i); + K_g_triplets.push_back(trip(row_pos * 2, col_pos * 2, K_m(j * 2, k * 2))); + K_g_triplets.push_back(trip(row_pos * 2 + 1, col_pos * 2, K_m(j * 2 + 1, k * 2))); + K_g_triplets.push_back(trip(row_pos * 2 + 1, col_pos * 2 + 1, K_m(j * 2 + 1, k * 2 + 1))); + K_g_triplets.push_back(trip(row_pos * 2, col_pos * 2 + 1, K_m(j * 2, k * 2 + 1))); + // we don't have to fill all because the matrix is symetric. + } + } + } + // FIXING SOME PINS: + // - if there are no pins (or only one pin) selected solve the system without the nullspace solution. + // - if there are some pins selected, delete all colums, rows that refer to this pins + // set the diagonal element of these pins to 1 + the rhs to zero + // (? is it possible to fix in the inner of the face? for sure for fem, but lscm could have some problems) + // (we also need some extra variables to see if the pins come from user) + + // fixing some points + // allthough only internal forces are applied there has to be locked + // at least 3 degrees of freedom to stop the mesh from pure rotation and pure translation + // std::vector fixed_dof; + // fixed_dof.push_back(this->triangles(0, 0) * 2); //x0 + // fixed_dof.push_back(this->triangles(0, 0) * 2 + 1); //y0 + // fixed_dof.push_back(this->triangles(1, 0) * 2 + 1); // y1 + + // align flat mesh to fixed edge + // Vector2 edge = this->flat_vertices.col(this->triangles(1, 0)) - + // this->flat_vertices.col(this->triangles(0, 0)); + // edge.normalize(); + // Eigen::Matrix rot; + // rot << edge.x(), edge.y(), -edge.y(), edge.x(); + // this->flat_vertices = rot * this->flat_vertices; + + // // return true if triplet row / col is in fixed_dof + // auto is_in_fixed_dof = [fixed_dof](const trip & element) -> bool { + // return ( + // (std::find(fixed_dof.begin(), fixed_dof.end(), element.row()) != fixed_dof.end()) or + // (std::find(fixed_dof.begin(), fixed_dof.end(), element.col()) != fixed_dof.end())); + // }; + // std::cout << "size of triplets: " << K_g_triplets.size() << std::endl; + // K_g_triplets.erase( + // std::remove_if(K_g_triplets.begin(), K_g_triplets.end(), is_in_fixed_dof), + // K_g_triplets.end()); + // std::cout << "size of triplets: " << K_g_triplets.size() << std::endl; + // for (long fixed: fixed_dof) + // { + // K_g_triplets.push_back(trip(fixed, fixed, 1.)); + // rhs[fixed] = 0; + // } + + // for (long i=0; i< this->vertices.cols() * 2; i++) + // K_g_triplets.push_back(trip(i, i, 0.01)); + + // lagrange multiplier + 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))); + } + + // project out the nullspace solution: + + // Eigen::VectorXd nullspace1(this->flat_vertices.cols() * 2); + // Eigen::VectorXd nullspace2(this->flat_vertices.cols() * 2); + // nullspace1.setZero(); + // nullspace2.setOnes(); + // for (long i= 0; i < this->flat_vertices.cols(); i++) + // { + // nullspace1(i) = 1; + // nullspace2(i) = 0; + // } + // nullspace1.normalize(); + // nullspace2.normalize(); + // rhs -= nullspace1.dot(rhs) * nullspace1; + // rhs -= nullspace2.dot(rhs) * nullspace2; + + K_g.setFromTriplets(K_g_triplets.begin(), K_g_triplets.end()); + // rhs += K_g * Eigen::VectorXd::Ones(K_g.rows()); + + // solve linear system (privately store the value for guess in next step) + Eigen::ConjugateGradient solver; + solver.setTolerance(0.0000001); + solver.compute(K_g); + this->sol = solver.solveWithGuess(-rhs, this->sol); + this->set_shift(this->sol.head(this->vertices.cols() * 2) * weight); + this->set_q_l_m(); +} + + +////////////////////////////////////////////////////////////////////////// +///////////////// L.S.C.M ///////////// +////////////////////////////////////////////////////////////////////////// +void LscmRelax::lscm() +{ + this->set_q_l_g(); + std::vector triple_list; + long i; + double x21, x31, y31, x32; + + // 1. create the triplet list (t * 2, v * 2) + for(i=0; itriangles.cols(); i++) + { + x21 = this->q_l_g(i, 0); + x31 = this->q_l_g(i, 1); + y31 = this->q_l_g(i, 2); + x32 = x31 - x21; + + triple_list.push_back(trip(2 * i, this->new_order[this->triangles(0, i)] * 2, x32)); + triple_list.push_back(trip(2 * i, this->new_order[this->triangles(0, i)] * 2 + 1, -y31)); + triple_list.push_back(trip(2 * i, this->new_order[this->triangles(1, i)] * 2, -x31)); + triple_list.push_back(trip(2 * i, this->new_order[this->triangles(1, i)] * 2 + 1, y31)); + triple_list.push_back(trip(2 * i, this->new_order[this->triangles(2, i)] * 2, x21)); + + triple_list.push_back(trip(2 * i + 1, this->new_order[this->triangles(0, i)] * 2, y31)); + triple_list.push_back(trip(2 * i + 1, this->new_order[this->triangles(0, i)] * 2 + 1, x32)); + triple_list.push_back(trip(2 * i + 1, this->new_order[this->triangles(1, i)] * 2, -y31)); + triple_list.push_back(trip(2 * i + 1, this->new_order[this->triangles(1, i)] * 2 + 1, -x31)); + triple_list.push_back(trip(2 * i + 1, this->new_order[this->triangles(2, i)] * 2 + 1, x21)); + + } + // 2. divide the triplets in matrix(unknown part) and rhs(known part) and reset the position + std::vector rhs_triplets; + std::vector mat_triplets; + for (auto triplet: triple_list) + { + if (triplet.col() > (this->vertices.cols() - this->fixed_pins.size()) * 2 - 1) + rhs_triplets.push_back(triplet); + else + mat_triplets.push_back(triplet); + } + + // 3. create a rhs_pos vector + Eigen::VectorXd rhs_pos(this->vertices.cols() * 2); + rhs_pos.setZero(); + for (auto index: this->fixed_pins) + { + rhs_pos[this->new_order[index] * 2] = this->flat_vertices(0, index); //TODO: not yet set + rhs_pos[this->new_order[index] * 2 + 1] = this->flat_vertices(1, index); + } + + // 4. fill a sparse matrix and calculdate the rhs + Eigen::VectorXd rhs(this->triangles.cols() * 2); // maybe use a sparse vector + spMat B(this->triangles.cols() * 2, this->vertices.cols() * 2); + B.setFromTriplets(rhs_triplets.begin(), rhs_triplets.end()); + rhs = B * rhs_pos; + + // 5. create the lhs matrix + spMat A(this->triangles.cols() * 2, (this->vertices.cols() - this->fixed_pins.size()) * 2); + A.setFromTriplets(mat_triplets.begin(), mat_triplets.end()); + + // 6. solve the system and set the flatted coordinates + // Eigen::SparseQR > solver; + Eigen::LeastSquaresConjugateGradient solver; + Eigen::VectorXd sol(this->vertices.size() * 2); + solver.compute(A); + sol = solver.solve(-rhs); + + // TODO: create function, is needed also in the fem step + this->set_position(sol); + this->set_q_l_m(); + this->transform(true); + // this->rotate_by_min_bound_area(); + this->set_q_l_m(); + +} + +void LscmRelax::set_q_l_g() +{ + // get the coordinates of a triangle in local coordinates from the 3d mesh + // x1, y1, y2 = 0 + // -> vector + this->q_l_g.resize(this->triangles.cols(), 3); + for (long i = 0; i < this->triangles.cols(); i++) + { + Vector3 r1 = this->vertices.col(this->triangles(0, i)); + Vector3 r2 = this->vertices.col(this->triangles(1, i)); + Vector3 r3 = this->vertices.col(this->triangles(2, i)); + Vector3 r21 = r2 - r1; + Vector3 r31 = r3 - r1; + double r21_norm = r21.norm(); + r21.normalize(); + // if triangle is fliped this gives wrong results? + this->q_l_g.row(i) << r21_norm, r31.dot(r21), r31.cross(r21).norm(); + } +} + +void LscmRelax::set_q_l_m() +{ + // get the coordinates of a triangle in local coordinates from the 2d map + // x1, y1, y2 = 0 + // -> vector + this->q_l_m.resize(this->triangles.cols(), 3); + for (long i = 0; i < this->triangles.cols(); i++) + { + Vector2 r1 = this->flat_vertices.col(this->triangles(0, i)); + Vector2 r2 = this->flat_vertices.col(this->triangles(1, i)); + Vector2 r3 = this->flat_vertices.col(this->triangles(2, i)); + Vector2 r21 = r2 - r1; + Vector2 r31 = r3 - r1; + double r21_norm = r21.norm(); + r21.normalize(); + // if triangle is fliped this gives wrong results! + this->q_l_m.row(i) << r21_norm, r31.dot(r21), -(r31.x() * r21.y() - r31.y() * r21.x()); + } +} + +void LscmRelax::set_fixed_pins() +{ + // if less then one fixed pin is set find two by an automated algorithm and align them to y = 0 + // if more then two pins are choosen find a leastsquare-plane and project the points on it + // insert the points in the flat-vertices vector + if (this->fixed_pins.size() == 0) + this->fixed_pins.push_back(0); + if (this->fixed_pins.size() == 1) + { + + double dist; + this->fixed_pins.push_back(get_max_distance(this->vertices.col(this->fixed_pins[0]), this->vertices, dist)); + this->flat_vertices.col(this->fixed_pins[0]) = Vector2(0, 0); + this->flat_vertices.col(this->fixed_pins[1]) = Vector2(dist, 0); + } + std::sort(this->fixed_pins.begin(), this->fixed_pins.end()); + // not yet working + // if (this->fixed_pins.size() > 2) + // { + // std::vector fixed_3d_points; + // for (unsigned int index: this->fixed_pins) + // fixed_3d_points.push_back(this->vertices[index]); + // std::vector maped_points = map_to_2D(fixed_3d_points); + // unsigned int i = 0; + // for (unsigned int index: this->fixed_pins) + // { + // this->flat_vertices[i] = maped_points[i]; + // i++; + // } + // } +} + + +ColMat LscmRelax::get_flat_vertices_3D() +{ + ColMat mat = this->flat_vertices.transpose(); + ColMat mat3d(mat.rows(), 3); + mat3d << mat, ColMat::Zero(mat.rows()); + return mat3d; +} + +void LscmRelax::set_position(Eigen::VectorXd sol) +{ + for (long i=0; i < this->vertices.size(); i++) + { + if (sol.size() > i * 2 + 1) + this->flat_vertices.col(this->old_order[i]) << sol[i * 2], sol[i * 2 + 1]; + } +} + +void LscmRelax::set_shift(Eigen::VectorXd sol) +{ + for (long i=0; i < this->vertices.size(); i++) + { + if (sol.size() > i * 2 + 1) + this->flat_vertices.col(i) += Vector2(sol[i * 2], sol[i * 2 + 1]); + } +} + +double LscmRelax::get_area() +{ + double area = 0; + for(long i = 0; i < this->triangles.cols(); i++) + { + area += this->q_l_g(i, 0) * this->q_l_g(i, 2); + } + return area / 2; +} + +double LscmRelax::get_flat_area() +{ + double area = 0; + for(long i = 0; i < this->triangles.cols(); i++) + { + area += this->q_l_m(i, 0) * this->q_l_m(i, 2); + } + return area / 2; +} + +void LscmRelax::transform(bool scale) +{ + // assuming q_l_m and flat_vertices are set + Vector2 weighted_center, center; + weighted_center.setZero(); + double flat_area = 0; + double global_area = 0; + double element_area; + for(long i = 0; i < this->triangles.cols(); i++) + { + global_area += this->q_l_g(i, 0) * this->q_l_g(i, 2) / 2; + element_area = this->q_l_m(i, 0) * this->q_l_m(i, 2) / 2; + for (int j=0; j < 3; j++) + weighted_center += this->flat_vertices.col(this->triangles(j, i)) * element_area / 3; + flat_area += element_area; + } + center = weighted_center / flat_area; + for (long i = 0; i < this->flat_vertices.cols(); i++) + this->flat_vertices.col(i) -= center; + if (scale) + this->flat_vertices *= std::pow(global_area / flat_area, 0.5); + this->set_q_l_m(); +} + +void LscmRelax::rotate_by_min_bound_area() +{ + int n = 100; + double phi; + double min_phi = 0; + double min_area = 0; + bool x_dominant; + // rotate vector by 90 degree and find min area + for (int i = 0; i < n + 1; i++ ) + { + phi = i * M_PI / n; + Eigen::VectorXd x_proj = this->flat_vertices.transpose() * Vector2(std::cos(phi), std::sin(phi)); + Eigen::VectorXd y_proj = this->flat_vertices.transpose() * Vector2(-std::sin(phi), std::cos(phi)); + double x_distance = x_proj.maxCoeff() - x_proj.minCoeff(); + double y_distance = y_proj.maxCoeff() - y_proj.minCoeff(); + double area = x_distance * y_distance; + if (min_area == 0 or area < min_area) + { + min_area = area; + min_phi = phi; + x_dominant = x_distance > y_distance; + } + Eigen::Matrix rot; + min_phi += x_dominant * M_PI / 2; + rot << std::cos(min_phi), std::sin(min_phi), -std::sin(min_phi), std::cos(min_phi); + this->flat_vertices = rot * this->flat_vertices; + } +} + +std::vector LscmRelax::get_fem_fixed_pins() +{ + long min_x_index = 0; + double min_x = this->vertices(0, 0); + for (long i=0; iflat_vertices.cols(); i++) + { + // get index with min x + if (this->flat_vertices(0, i) < min_x) + { + min_x = this->flat_vertices(0, i); + min_x_index = i; + } + } + double min_y = this->flat_vertices(1, min_x_index); + long max_x_index = 0; + double max_x = 0; + for (long i=0; iflat_vertices.cols(); i++) + { + // get index with min x + double d_x = (this->flat_vertices(0, i) - min_x); + double d_y = (this->flat_vertices(1, i) - min_y); + if (d_x * d_x - d_y * d_y > max_x) + { + max_x = d_x * d_x - d_y * d_y; + max_x_index = i; + } + } + return std::vector{min_x_index * 2, min_x_index * 2 + 1, max_x_index * 2 + 1}; +} diff --git a/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h new file mode 100644 index 0000000000..7795f6c7c0 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlatteningLscmRelax.h @@ -0,0 +1,108 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +// LeastSquareConformalMapping + fem relaxing +// ------------------------------------------ +// +#ifndef UNWRAP_H +#define UNWRAP_H + +// 1: local coordinates 2d representation q_l_0 +// 2: least square conformal map -> flat_vertices_0 +// 3: local coordinates of mapped mesh q_l_1 +// 4: diff in local coordinates -> forces R.B^T.B.(x1-x0) +// 5: stiffnes mat K +// 6: K.u=forces ->u +// 7: x1, y1 += w * u + +#include +#include +#include + +#include +#include + +#include "MeshFlatteningLscmRelax.h" + +typedef Eigen::SparseMatrix spMat; + + +typedef Eigen::Vector3d Vector3; +typedef Eigen::Vector2d Vector2; + +template +using ColMat = Eigen::Matrix; + +template +using RowMat = Eigen::Matrix; + +class LscmRelax{ +private: + ColMat q_l_g; // the position of the 3d triangles at there locale coord sys + ColMat q_l_m; // the mapped position in local coord sys + + void set_q_l_g(); + void set_q_l_m(); + void set_fixed_pins(); + void set_position(Eigen::VectorXd); + void set_shift(Eigen::VectorXd); + + std::vector new_order; + std::vector old_order; + + Eigen::Matrix C; + Eigen::VectorXd sol; + + std::vector get_fem_fixed_pins(); + +public: + LscmRelax( + RowMat vertices, + RowMat triangles, + std::vector fixed_pins); + + std::vector fixed_pins; + RowMat vertices; + RowMat triangles; + RowMat flat_vertices; + ColMat rhs; + Eigen::MatrixXd MATRIX; + + double nue=0.0; + double elasticity=1.; + + void lscm(); + void relax(double); + + ColMat get_flat_vertices_3D(); + + void rotate_by_min_bound_area(); + void transform(bool scale=false); + + double get_area(); + double get_flat_area(); + +}; + + +#endif diff --git a/src/Mod/MeshPart/App/MeshFlatteningNurbs.cpp b/src/Mod/MeshPart/App/MeshFlatteningNurbs.cpp new file mode 100644 index 0000000000..7604e5e71f --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlatteningNurbs.cpp @@ -0,0 +1,403 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#include "MeshFlatteningNurbs.h" +#include +#include "math.h" + +namespace nurbs{ + +double divide(double a, double b) +{ + if (fabs(b) < 10e-14) + return 0; + else + return a / b; +} + +// DE BOOR ALGORITHM FROM OPENGLIDER +std::function get_basis(int degree, int i, Eigen::VectorXd knots) + // Return a basis_function for the given degree """ +{ + if (degree == 0) + { + return [degree, i, knots](double t) + { + // The basis function for degree = 0 as per eq. 7 + double t_this = knots[i]; + double t_next = knots[i+1]; + if (t == knots[0]) + return (double)(t_next >= t and t >= t_this); + return (double)(t_next >= t and t > t_this); + }; + } + else + { + return [degree, i, knots](double t) + { + // The basis function for degree > 0 as per eq. 8 + if (t < knots[0]) + return get_basis(degree, i, knots)(knots[0]); + else if (t > knots[knots.rows() - 1]) + return get_basis(degree, i, knots)(knots[knots.rows() - 1]); + double out = 0.; + double t_this = knots[i]; + double t_next = knots[i + 1]; + double t_precog = knots[i + degree]; + double t_horizon = knots[i + degree + 1]; + if (t_this == t_horizon) + return 0.; + + double bottom = (t_precog - t_this); + out = divide(t - t_this, bottom) * get_basis(degree - 1, i, knots)(t); + + bottom = (t_horizon - t_next); + out += divide(t_horizon - t, bottom) * get_basis(degree - 1, i + 1, knots)(t); + return out; + }; + } +}; + + +std::function get_basis_derivative(int order, int degree, int i, Eigen::VectorXd knots) + // Return the derivation of the basis function + // order of basis function + // degree of basis function + // + // knots sequence +{ + if (order == 1) + { + return [degree, i, knots, order](double t) + { + double out = 0; + if (not (knots[i + degree] - knots[i] == 0)) + { + out += get_basis(degree - 1, i, knots)(t) * + degree / (knots[i + degree] - knots[i]); + } + if (not (knots[i + degree + 1] - knots[i + 1] == 0)) + { + out -= get_basis(degree - 1, i + 1, knots)(t) * + degree / (knots[i + degree + 1] - knots[i + 1]); + } + return out; + }; + } + else + { + return [degree, i, knots, order](double t) + { + double out = 0; + if (not (knots[i + degree] - knots[i] == 0)) + { + out += get_basis_derivative(order - 1, degree - 1, i, knots)(t) * + degree / (knots[i + degree] - knots[i]); + } + if (not (knots[i + degree + 1] - knots[i + 1] == 0)) + { + out -= get_basis_derivative(order - 1, degree - 1, i + 1, knots)(t) * + degree / (knots[i + degree + 1] - knots[i + 1]); + } + return out; + }; + } +} + + +NurbsBase2D::NurbsBase2D(Eigen::VectorXd u_knots, Eigen::VectorXd v_knots, + Eigen::VectorXd weights, + int degree_u, int degree_v) +{ + // assert(weights.size() == u_knots.size() * v_knots.size()); + this->u_knots = u_knots; + this->v_knots = v_knots; + this->weights = weights; + this->degree_u = degree_u; + this->degree_v = degree_v; + for (int u_i = 0; u_i < u_knots.size() - degree_u - 1; u_i ++) + this->u_functions.push_back(get_basis(degree_u, u_i, u_knots)); + for (int v_i = 0; v_i < v_knots.size() - degree_v - 1; v_i ++) + this->v_functions.push_back(get_basis(degree_v, v_i, v_knots)); +} + +Eigen::VectorXd NurbsBase2D::getInfluenceVector(Eigen::Vector2d u) +{ + Eigen::VectorXd n_u, n_v; + double sum_weights = 0; + Eigen::VectorXd infl(this->u_functions.size() * this->v_functions.size()); + int i = 0; + + n_u.resize(this->u_functions.size()); + n_v.resize(this->v_functions.size()); + + for (unsigned int i = 0; i < this->u_functions.size(); i ++) + n_u[i] = this->u_functions[i](u.x()); + for (unsigned int i = 0; i < this->v_functions.size(); i ++) + n_v[i] = this->v_functions[i](u.y()); + + for (unsigned int u_i = 0; u_i < this->u_functions.size(); u_i++) + { + for (unsigned int v_i = 0; v_i < this->v_functions.size(); v_i++) + { + sum_weights += weights[i] * n_u[u_i] * n_v[v_i]; + infl[i] = weights[i] * n_u[u_i] * n_v[v_i]; + i ++; + } + } + return infl / sum_weights; +} + +void add_triplets(Eigen::VectorXd values, double row, std::vector &triplets) +{ + for (unsigned int i=0; i < values.size(); i++) + { + if (values(i) != 0.) + triplets.push_back(trip(row, i, values(i))); + } +} + +spMat NurbsBase2D::getInfluenceMatrix(Eigen::Matrix U) +{ + std::vector triplets; + for (int row_index; 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()); + return mat; +} + +void NurbsBase2D::computeFirstDerivatives() +{ + for (unsigned int u_i = 0; u_i < u_functions.size(); u_i ++) + this->Du_functions.push_back(get_basis_derivative(1, this->degree_u, u_i, this->u_knots)); + for (unsigned int v_i = 0; v_i < v_functions.size(); v_i ++) + this->Dv_functions.push_back(get_basis_derivative(1, this->degree_v, v_i, this->v_knots)); +} + +void NurbsBase2D::computeSecondDerivatives() +{ + for (unsigned int u_i = 0; u_i < u_functions.size(); u_i ++) + this->DDu_functions.push_back(get_basis_derivative(2, this->degree_u, u_i, this->u_knots)); + for (unsigned int v_i = 0; v_i < v_functions.size(); v_i ++) + this->DDv_functions.push_back(get_basis_derivative(2, this->degree_v, v_i, this->v_knots)); +} + +Eigen::VectorXd NurbsBase2D::getDuVector(Eigen::Vector2d u) +{ + Eigen::VectorXd A1, A2; + double C1, C2; + A1.resize(this->u_functions.size() * this->v_functions.size()); + A2.resize(this->u_functions.size() * this->v_functions.size()); + double A3 = 0; + double A5 = 0; + int i = 0; + Eigen::VectorXd n_u, n_v, Dn_u; + n_u.resize(this->u_functions.size()); + Dn_u.resize(this->u_functions.size()); + n_v.resize(this->v_functions.size()); + for (unsigned int u_i=0; u_i < this->u_functions.size(); u_i++) + { + // std::cout << "u_i: " << u_i << " , n_u: " << n_u.size() + // << " , Dn_u: " << Dn_u.size() << std::endl; + n_u[u_i] = this->u_functions[u_i](u.x()); + Dn_u[u_i] = this->Du_functions[u_i](u.x()); + } + for (unsigned int v_i=0; v_i < this->v_functions.size(); v_i++) + { + n_v[v_i] = this->v_functions[v_i](u.y()); + // std::cout << v_i << std::endl; + } + + for (unsigned int u_i=0; u_i < this->u_functions.size(); u_i++) + { + for (unsigned int v_i=0; v_i < this->v_functions.size(); v_i++) + { + C1 = weights[i] * Dn_u[u_i] * n_v[v_i]; + C2 = weights[i] * n_u[u_i] * n_v[v_i]; + A1[i] = C1; + A2[i] = C2; + A3 += C2; + A5 += C1; + i ++; + } + } + return (A1 * A3 - A2 * A5) / A3 / A3; +} + +Eigen::VectorXd NurbsBase2D::getDvVector(Eigen::Vector2d u) +{ + Eigen::VectorXd A1, A2; + double C1, C2; + A1.resize(this->u_functions.size() * this->v_functions.size()); + A2.resize(this->u_functions.size() * this->v_functions.size()); + double A3 = 0; + double A5 = 0; + int i = 0; + Eigen::VectorXd n_u, n_v, Dn_v; + n_u.resize(this->u_functions.size()); + Dn_v.resize(this->v_functions.size()); + n_v.resize(this->v_functions.size()); + for (unsigned int u_i=0; u_i < this->u_functions.size(); u_i++) + { + n_u[u_i] = this->u_functions[u_i](u.x()); + } + for (unsigned int v_i=0; v_i < this->v_functions.size(); v_i++) + { + n_v[v_i] = this->v_functions[v_i](u.y()); + Dn_v[v_i] = this->Dv_functions[v_i](u.y()); + } + + for (unsigned int u_i=0; u_i < this->u_functions.size(); u_i++) + { + for (unsigned int v_i=0; v_i < this->v_functions.size(); v_i++) + { + C1 = weights[i] * Dn_v[v_i] * n_u[u_i]; + C2 = weights[i] * n_v[v_i] * n_u[u_i]; + A1[i] = C1; + A2[i] = C2; + A3 += C2; + A5 += C1; + i ++; + } + } + return (A1 * A3 - A2 * A5) / A3 / A3; +} + + +spMat NurbsBase2D::getDuMatrix(Eigen::Matrix U) +{ + std::vector triplets; + for (int row_index; 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()); + return mat; +} + +spMat NurbsBase2D::getDvMatrix(Eigen::Matrix U) +{ + std::vector triplets; + for (int row_index; 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; +} + +NurbsBase1D::NurbsBase1D(Eigen::VectorXd u_knots, Eigen::VectorXd weights, int degree_u) +{ + this->u_knots = u_knots; + this->weights = weights; + this->degree_u = degree_u; + for (int u_i = 0; u_i < u_knots.size() - degree_u - 1; u_i ++) + this->u_functions.push_back(get_basis(degree_u, u_i, u_knots)); +} + +Eigen::VectorXd NurbsBase1D::getInfluenceVector(double u) +{ + Eigen::VectorXd n_u; + double sum_weights = 0; + Eigen::VectorXd infl(this->u_functions.size()); + + n_u.resize(this->u_functions.size()); + + for (unsigned int i = 0; i < this->u_functions.size(); i ++) + n_u[i] = this->u_functions[i](u); + + for (unsigned int u_i = 0; u_i < this->u_functions.size(); u_i++) + { + sum_weights += weights[u_i] * n_u[u_i]; + infl[u_i] = weights[u_i] * n_u[u_i]; + } + return infl / sum_weights; +} + +spMat NurbsBase1D::getInfluenceMatrix(Eigen::VectorXd u) +{ + std::vector triplets; + for (int row_index; 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()); + return mat; +} + +void NurbsBase1D::computeFirstDerivatives() +{ + for (unsigned int u_i = 0; u_i < u_functions.size(); u_i ++) + this->Du_functions.push_back(get_basis_derivative(1, this->degree_u, u_i, this->u_knots)); +} + +void NurbsBase1D::computeSecondDerivatives() +{ + for (unsigned int u_i = 0; u_i < u_functions.size(); u_i ++) + this->DDu_functions.push_back(get_basis_derivative(2, this->degree_u, u_i, this->u_knots)); +} + + +Eigen::VectorXd NurbsBase1D::getDuVector(double u) +{ + Eigen::VectorXd A1, A2; + double C1, C2; + double C3 = 0; + double C4 = 0; + int i = 0; + A1.resize(this->u_functions.size()); + A2.resize(this->u_functions.size()); + Eigen::VectorXd n_u, Dn_u; + n_u.resize(this->u_functions.size()); + Dn_u.resize(this->u_functions.size()); + + for (unsigned u_i=0; u_i < this->u_functions.size(); u_i++) + { + n_u[u_i] = this->u_functions[u_i](u); + Dn_u[u_i] = this->Du_functions[u_i](u); + } + + for (unsigned int u_i=0; u_i < this->Du_functions.size(); u_i++) + { + C1 = weights[i] * Dn_u[u_i]; + C2 = weights[i] * n_u[u_i]; + C3 += C1; + C4 += C2; + + A1[i] = C1; + A2[i] = C2; + i ++; + } + return (A1 * C4 - A2 * C3) / C4 / C4 ; +} + + +spMat NurbsBase1D::getDuMatrix(Eigen::VectorXd U) +{ + std::vector triplets; + for (int row_index; 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; +} + +} diff --git a/src/Mod/MeshPart/App/MeshFlatteningNurbs.h b/src/Mod/MeshPart/App/MeshFlatteningNurbs.h new file mode 100644 index 0000000000..073b17b543 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlatteningNurbs.h @@ -0,0 +1,96 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + + +#ifndef NURBS_H +#define NURBS_H + +#include +#include +#include + +namespace nurbs{ + +typedef Eigen::Triplet trip; +typedef Eigen::SparseMatrix spMat; + +struct NurbsBase2D +{ + // + NurbsBase2D(){;}; + NurbsBase2D(Eigen::VectorXd u_knots, Eigen::VectorXd v_knots, + Eigen::VectorXd weights, + int degree_u=3, int degree_v=3); + int degree_u; + int degree_v; + Eigen::VectorXd u_knots; + Eigen::VectorXd v_knots; + Eigen::VectorXd weights; + + std::vector> u_functions; + std::vector> v_functions; + + std::vector> Du_functions; + std::vector> Dv_functions; + + std::vector> DDu_functions; + std::vector> DDv_functions; + + void computeFirstDerivatives(); + void computeSecondDerivatives(); + + Eigen::VectorXd getInfluenceVector(Eigen::Vector2d u); + spMat getInfluenceMatrix(Eigen::Matrix U); + + Eigen::VectorXd getDuVector(Eigen::Vector2d u); + spMat getDuMatrix(Eigen::Matrix U); + + Eigen::VectorXd getDvVector(Eigen::Vector2d u); + spMat getDvMatrix(Eigen::Matrix U); +}; + +struct NurbsBase1D +{ + NurbsBase1D(){;}; + NurbsBase1D(Eigen::VectorXd u_knots, Eigen::VectorXd weights, int degree_u=3); + int degree_u = degree_u; + Eigen::VectorXd u_knots; + Eigen::VectorXd weights; + std::vector> u_functions; + std::vector> Du_functions; + std::vector> DDu_functions; + + void computeFirstDerivatives(); + void computeSecondDerivatives(); + + Eigen::VectorXd getInfluenceVector(double u); + spMat getInfluenceMatrix(Eigen::VectorXd u); + + Eigen::VectorXd getDuVector(double u); + spMat getDuMatrix(Eigen::VectorXd u); + +}; + +} + +#endif diff --git a/src/Mod/MeshPart/App/MeshFlatteningPy.cpp b/src/Mod/MeshPart/App/MeshFlatteningPy.cpp new file mode 100644 index 0000000000..b23d38d255 --- /dev/null +++ b/src/Mod/MeshPart/App/MeshFlatteningPy.cpp @@ -0,0 +1,125 @@ +/*************************************************************************** + * Copyright (c) 2017 Lorenz Lechner * + * * + * This file is part of the FreeCAD CAx development system. * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Library General Public * + * License as published by the Free Software Foundation; either * + * version 2 of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU Library General Public License for more details. * + * * + * You should have received a copy of the GNU Library General Public * + * License along with this library; see the file COPYING.LIB. If not, * + * write to the Free Software Foundation, Inc., 59 Temple Place, * + * Suite 330, Boston, MA 02111-1307, USA * + * * + ***************************************************************************/ + + +#include "PreCompiled.h" +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "MeshFlattening.h" +#include "MeshFlatteningLscmRelax.h" +#include "MeshFlatteningNurbs.h" + +#include +#include + + + +namespace py = pybind11; + +// void FaceUnwrapper_constructor(FaceUnwrapper& instance, Part::TopoShapeFacePy Face) +FaceUnwrapper* FaceUnwrapper_constructor(py::object face) +{ + if (PyObject_TypeCheck(face.ptr(), &(Part::TopoShapeFacePy::Type))) + { + const Part::TopoShapeFacePy* f = static_cast(face.ptr()); + const TopoDS_Face& myFace = TopoDS::Face(f->getTopoShapePtr()->getShape()); + return new FaceUnwrapper(myFace); + } + else + throw std::invalid_argument("FaceUnwrapper should be initialized with Part.Face"); +} + +ColMat interpolateNurbsFacePy(FaceUnwrapper& instance, py::object face) +{ + std::cout << face.ptr()->ob_type->tp_name << std::endl; + std::cout << Part::TopoShapeFacePy::Type.tp_name << std::endl; + if (PyObject_TypeCheck(face.ptr(), &(Part::TopoShapeFacePy::Type))) + { + const Part::TopoShapeFacePy* f = static_cast(face.ptr()); + const TopoDS_Face& myFace = TopoDS::Face(f->getTopoShapePtr()->getShape()); + return instance.interpolateNurbsFace(myFace); + } + else + throw std::invalid_argument("FaceUnwrapper.interpolateNurbs should be initialized with Part.Face"); +} + +PYBIND11_MODULE(flatmesh, m) +{ + m.doc() = "functions to unwrapp faces/ meshes"; + + py::class_(m, "LscmRelax") + .def(py::init, ColMat, std::vector>()) + .def("lscm", &LscmRelax::lscm) + .def("relax", &LscmRelax::relax) + .def("rotate_by_min_bound_area", &LscmRelax::rotate_by_min_bound_area) + .def("transform", &LscmRelax::transform) + .def_readonly("rhs", &LscmRelax::rhs) + .def_readonly("MATRIX", &LscmRelax::MATRIX) + .def_property_readonly("area", &LscmRelax::get_area) + .def_property_readonly("flat_area", &LscmRelax::get_flat_area) + .def_property_readonly("flat_vertices", [](LscmRelax& L){return L.flat_vertices.transpose();}, py::return_value_policy::copy) + .def_property_readonly("flat_vertices_3D", &LscmRelax::get_flat_vertices_3D); + + py::class_(m, "NurbsBase2D") + .def(py::init()) + .def("computeFirstDerivatives", &nurbs::NurbsBase2D::computeFirstDerivatives) + .def("getInfluenceVector", &nurbs::NurbsBase2D::getInfluenceVector) + .def("getInfluenceMatrix", &nurbs::NurbsBase2D::getInfluenceMatrix) + .def("getDuVector", &nurbs::NurbsBase2D::getDuVector) + .def("getDuMatrix", &nurbs::NurbsBase2D::getDuMatrix) + .def("getDvVector", &nurbs::NurbsBase2D::getDvVector) + .def("getDvMatrix", &nurbs::NurbsBase2D::getDvMatrix); + + py::class_(m, "NurbsBase1D") + .def(py::init()) + .def("computeFirstDerivatives", &nurbs::NurbsBase1D::computeFirstDerivatives) + .def("getInfluenceVector", &nurbs::NurbsBase1D::getInfluenceVector) + .def("getInfluenceMatrix", &nurbs::NurbsBase1D::getInfluenceMatrix) + .def("getDuVector", &nurbs::NurbsBase1D::getDuVector) + .def("getDuMatrix", &nurbs::NurbsBase1D::getDuMatrix); + + py::class_(m, "FaceUnwrapper") + .def(py::init(&FaceUnwrapper_constructor)) + .def("findFlatNodes", &FaceUnwrapper::findFlatNodes) + .def("interpolateNurbsFace", &interpolateNurbsFacePy) + .def_readonly("tris", &FaceUnwrapper::tris) + .def_readonly("nodes", &FaceUnwrapper::xyz_nodes) + .def_readonly("uv_nodes", &FaceUnwrapper::uv_nodes) + .def_readonly("ze_nodes", &FaceUnwrapper::ze_nodes) + .def_readonly("ze_poles", &FaceUnwrapper::ze_poles) + .def_readonly("A", &FaceUnwrapper::A); +};