MeshFlattening: initial commit

This commit is contained in:
looooo
2017-09-07 18:39:13 +02:00
committed by wmayer
parent 05b8004001
commit ce33e2dbe2
8 changed files with 1533 additions and 0 deletions

View File

@@ -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})

View File

@@ -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 <Poly_Triangulation.hxx>
#include <BRep_Tool.hxx>
#include <Geom_Surface.hxx>
#include <Geom_BSplineSurface.hxx>
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<long> 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<double, 3> 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<spMat > solver;
solver.compute(A);
ColMat<double, 2> ze_poles;
ColMat<double, 3> 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;
}

View File

@@ -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 <BRepTools.hxx>
#include <TopoDS_Face.hxx>
#include <vector>
#include <Eigen/Geometry>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/SparseCore>
#include <Eigen/QR>
typedef Eigen::Vector3d Vector3;
typedef Eigen::Vector2d Vector2;
template <typename type, unsigned int size>
using ColMat = Eigen::Matrix<type, Eigen::Dynamic, size>;
template <typename type, unsigned int size>
using RowMat = Eigen::Matrix<type, size, Eigen::Dynamic>;
typedef Eigen::Triplet<double> trip;
typedef Eigen::SparseMatrix<double> spMat;
class FaceUnwrapper{
nurbs::NurbsBase2D nu;
public:
FaceUnwrapper(const TopoDS_Face & face);
void findFlatNodes();
ColMat<double, 3> interpolateNurbsFace(const TopoDS_Face& face);
bool use_nurbs = true;
// the mesh
ColMat<long, 3> tris; // input
ColMat<long, 1> fixed_nodes; // input
ColMat<double, 3> xyz_nodes; // compute from uv_mesh (xyz = A * poles)
ColMat<double, 2> uv_nodes; // input
ColMat<double, 2> ze_nodes; // copute
// nurbs
ColMat<double, 2> ze_poles; // compute
spMat A; // mapping between nurbs(poles) and mesh(vertices) computed with nurbs-basis-functions and uv_mesh
//
};
#endif // MESHFLATTENING

View File

@@ -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 <Eigen/IterativeLinearSolvers>
#include <Eigen/SVD>
#include <iostream>
#include <algorithm>
#include <cmath>
#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<double> trip;
typedef Eigen::SparseMatrix<double> spMat;
ColMat<double, 2> map_to_2D(ColMat<double, 3> points)
{
ColMat<double, 4> mat(points.size(), 4);
mat << points, ColMat<double, 1>::Ones(points.size());
Eigen::JacobiSVD<Eigen::MatrixXd> 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<double, 3, 2> transform;
transform.col(0) = x;
transform.col(1) = y;
return points * transform;
}
unsigned int get_max_distance(Vector3 point, RowMat<double, 3> 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<double, 3> vertices,
RowMat<long, 3> triangles,
std::vector<long> 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<double, 3> d_q_l_g = this->q_l_m - this->q_l_g;
// for every triangle
Eigen::Matrix<double, 3, 6> B;
Eigen::Matrix<double, 2, 2> T;
Eigen::Matrix<double, 6, 6> K_m;
Eigen::Matrix<double, 6, 1> 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<trip> K_g_triplets;
Vector2 v1, v2, v3, v12, v23, v31;
long row_pos, col_pos;
double A;
rhs.setZero();
for (long i=0; i<this->triangles.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<long> 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<double, 2, 2> 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<spMat, Eigen::Lower> 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<trip> triple_list;
long i;
double x21, x31, y31, x32;
// 1. create the triplet list (t * 2, v * 2)
for(i=0; i<this->triangles.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<trip> rhs_triplets;
std::vector<trip> 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<spMat, Eigen::COLAMDOrdering<int> > solver;
Eigen::LeastSquaresConjugateGradient<spMat > 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<x2, x3, y3>
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<x2, x3, y3>
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<Vector3> fixed_3d_points;
// for (unsigned int index: this->fixed_pins)
// fixed_3d_points.push_back(this->vertices[index]);
// std::vector<Vector2> 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<double, 3> LscmRelax::get_flat_vertices_3D()
{
ColMat<double, 2> mat = this->flat_vertices.transpose();
ColMat<double, 3> mat3d(mat.rows(), 3);
mat3d << mat, ColMat<double, 1>::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<double, 2, 2> 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<long> LscmRelax::get_fem_fixed_pins()
{
long min_x_index = 0;
double min_x = this->vertices(0, 0);
for (long i=0; i<this->flat_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; i<this->flat_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<long>{min_x_index * 2, min_x_index * 2 + 1, max_x_index * 2 + 1};
}

View File

@@ -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 <vector>
#include <memory>
#include <tuple>
#include <Eigen/Geometry>
#include <Eigen/IterativeLinearSolvers>
#include "MeshFlatteningLscmRelax.h"
typedef Eigen::SparseMatrix<double> spMat;
typedef Eigen::Vector3d Vector3;
typedef Eigen::Vector2d Vector2;
template <typename type, unsigned int size>
using ColMat = Eigen::Matrix<type, Eigen::Dynamic, size>;
template <typename type, unsigned int size>
using RowMat = Eigen::Matrix<type, size, Eigen::Dynamic>;
class LscmRelax{
private:
ColMat<double, 3> q_l_g; // the position of the 3d triangles at there locale coord sys
ColMat<double, 3> 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<long> new_order;
std::vector<long> old_order;
Eigen::Matrix<double, 3, 3> C;
Eigen::VectorXd sol;
std::vector<long> get_fem_fixed_pins();
public:
LscmRelax(
RowMat<double, 3> vertices,
RowMat<long, 3> triangles,
std::vector<long> fixed_pins);
std::vector<long> fixed_pins;
RowMat<double, 3> vertices;
RowMat<long, 3> triangles;
RowMat<double, 2> flat_vertices;
ColMat<double, 1> rhs;
Eigen::MatrixXd MATRIX;
double nue=0.0;
double elasticity=1.;
void lscm();
void relax(double);
ColMat<double, 3> get_flat_vertices_3D();
void rotate_by_min_bound_area();
void transform(bool scale=false);
double get_area();
double get_flat_area();
};
#endif

View File

@@ -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 <iostream>
#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<double(double)> 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<double(double)> 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<trip> &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<double, Eigen::Dynamic, 2> U)
{
std::vector<trip> 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<double, Eigen::Dynamic, 2> U)
{
std::vector<trip> 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<double, Eigen::Dynamic, 2> U)
{
std::vector<trip> 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<trip> 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<trip> 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;
}
}

View File

@@ -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 <Eigen/Geometry>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/SparseCore>
namespace nurbs{
typedef Eigen::Triplet<double> trip;
typedef Eigen::SparseMatrix<double> 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<std::function<double(double)>> u_functions;
std::vector<std::function<double(double)>> v_functions;
std::vector<std::function<double(double)>> Du_functions;
std::vector<std::function<double(double)>> Dv_functions;
std::vector<std::function<double(double)>> DDu_functions;
std::vector<std::function<double(double)>> DDv_functions;
void computeFirstDerivatives();
void computeSecondDerivatives();
Eigen::VectorXd getInfluenceVector(Eigen::Vector2d u);
spMat getInfluenceMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U);
Eigen::VectorXd getDuVector(Eigen::Vector2d u);
spMat getDuMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> U);
Eigen::VectorXd getDvVector(Eigen::Vector2d u);
spMat getDvMatrix(Eigen::Matrix<double, Eigen::Dynamic, 2> 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<std::function<double(double)>> u_functions;
std::vector<std::function<double(double)>> Du_functions;
std::vector<std::function<double(double)>> 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

View File

@@ -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 <Mod/Part/App/TopoShapeFacePy.h>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Sparse>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/operators.h>
#include <pybind11/numpy.h>
#include <pybind11/eigen.h>
#include <memory>
#include <vector>
#include <tuple>
#include <map>
#include <stdexcept>
#include "MeshFlattening.h"
#include "MeshFlatteningLscmRelax.h"
#include "MeshFlatteningNurbs.h"
#include <TopoDS_Face.hxx>
#include <TopoDS.hxx>
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<Part::TopoShapeFacePy*>(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<double, 3> 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<Part::TopoShapeFacePy*>(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_<LscmRelax>(m, "LscmRelax")
.def(py::init<ColMat<double, 3>, ColMat<long, 3>, std::vector<long>>())
.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_<nurbs::NurbsBase2D>(m, "NurbsBase2D")
.def(py::init<Eigen::VectorXd, Eigen::VectorXd, Eigen::VectorXd, int, int>())
.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_<nurbs::NurbsBase1D>(m, "NurbsBase1D")
.def(py::init<Eigen::VectorXd, Eigen::VectorXd, int>())
.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_<FaceUnwrapper>(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);
};