FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2017 Lorenz Lechner                                     *
3
 *                                                                         *
4
 *   This file is part of the FreeCAD CAx development system.              *
5
 *                                                                         *
6
 *   This library is free software; you can redistribute it and/or         *
7
 *   modify it under the terms of the GNU Library General Public           *
8
 *   License as published by the Free Software Foundation; either          *
9
 *   version 2 of the License, or (at your option) any later version.      *
10
 *                                                                         *
11
 *   This library  is distributed in the hope that it will be useful,      *
12
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
14
 *   GNU Library General Public License for more details.                  *
15
 *                                                                         *
16
 *   You should have received a copy of the GNU Library General Public     *
17
 *   License along with this library; see the file COPYING.LIB. If not,    *
18
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
19
 *   Suite 330, Boston, MA  02111-1307, USA                                *
20
 *                                                                         *
21
 ***************************************************************************/
22

23
#include "PreCompiled.h"
24
#ifndef _PreComp_
25
#include <map>
26
#include <set>
27
#include <vector>
28

29
#include <BRep_Tool.hxx>
30
#include <Geom_BSplineSurface.hxx>
31
#include <Geom_Surface.hxx>
32
#include <Poly_Triangulation.hxx>
33
#include <Standard_Version.hxx>
34
#include <TColStd_Array1OfReal.hxx>
35
#include <TopLoc_Location.hxx>
36
#endif
37

38
#include "MeshFlattening.h"
39
#include "MeshFlatteningLscmRelax.h"
40

41

42
std::vector<ColMat<double, 3>> getBoundaries(ColMat<double, 3> vertices, ColMat<long, 3> tris)
43
{
44
    // get a hashtable for all edges
45
    // e: v1, v2, num
46
    std::map<std::set<long>, std::vector<long>> hash_map;
47
    std::vector<std::set<long>> hash_list;
48
    std::map<long, std::vector<long>> neighbour_map;
49
    std::vector<long> edge_vector_0;
50
    std::vector<std::vector<long>> edge_vector;
51

52

53
    for (long i = 0; i < tris.rows(); i++) {
54
        for (long j = 0; j < 3; j++) {
55
            long k = j + 1;
56
            if (k == 3) {
57
                k = 0;
58
            }
59
            long v1 = tris(i, j);
60
            long v2 = tris(i, k);
61
            std::set<long> hash {v1, v2};
62
            hash_list.push_back(hash);
63
            if (v1 < v2) {
64
                hash_map[hash] = std::vector<long> {v1, v2, 0};
65
            }
66
            else {
67
                hash_map[hash] = std::vector<long> {v2, v1, 0};
68
            }
69
        }
70
    }
71
    for (auto& hash : hash_list) {
72
        hash_map[hash][2] += 1;
73
    }
74

75
    for (auto& hash : hash_map) {
76
        if (hash.second[2] == 1) {
77
            long v0 = hash.second[0];
78
            long v1 = hash.second[1];
79

80
            neighbour_map[v0].push_back(v1);
81
            neighbour_map[v1].push_back(v0);
82
        }
83
    }
84

85

86
    while (neighbour_map.size() != 0) {
87
        long start_index = neighbour_map.begin()->first;
88
        long close_index = start_index;
89
        long next_index = neighbour_map[start_index][1];
90
        long temporary_next;
91
        edge_vector_0.clear();
92
        edge_vector_0.push_back(close_index);
93
        edge_vector_0.push_back(start_index);
94
        neighbour_map.erase(start_index);
95
        edge_vector_0.push_back(next_index);
96
        while (next_index != close_index) {
97
            temporary_next = neighbour_map[next_index][0];
98
            if (temporary_next != start_index) {
99
                start_index = next_index;
100
                next_index = temporary_next;
101
            }
102
            else {
103
                start_index = next_index;
104
                next_index = neighbour_map[start_index][1];
105
            }
106
            neighbour_map.erase(start_index);
107
            edge_vector_0.push_back(next_index);
108
        }
109
        edge_vector.push_back(edge_vector_0);
110
    }
111
    std::vector<ColMat<double, 3>> edges;
112
    for (auto& edge : edge_vector) {
113
        ColMat<double, 3> edge_vertices;
114
        edge_vertices.resize(edge.size(), 3);
115
        int i = 0;
116
        for (auto index : edge) {
117
            edge_vertices.row(i) = vertices.row(index);
118
            i++;
119
        }
120
        edges.push_back(edge_vertices);
121
    }
122
    return edges;
123
}
124

125
FaceUnwrapper::FaceUnwrapper(const TopoDS_Face& face)
126
{
127
    long i = 0;
128
    //  transform to nurbs:
129
    TopLoc_Location location;
130

131
    //  triangulate:
132
    const Handle(Poly_Triangulation)& triangulation = BRep_Tool::Triangulation(face, location);
133

134
    if (triangulation.IsNull()) {
135
        throw std::runtime_error("null triangulation in face construction");
136
    }
137

138
    Standard_Integer numNodes = triangulation->NbNodes();
139
    Standard_Integer numTriangles = triangulation->NbTriangles();
140

141
    //  compute uv coordinates
142
    if (triangulation->HasUVNodes()) {
143
        this->uv_nodes.resize(numNodes, 2);
144
        i = 0;
145
        for (Standard_Integer index = 1; index <= numNodes; ++index) {
146
            const gp_Pnt2d& _uv_node = triangulation->UVNode(index);
147
            this->uv_nodes.row(i) << _uv_node.X(), _uv_node.Y();
148
            i++;
149
        }
150
    }
151
    //
152
    this->xyz_nodes.resize(numNodes, 3);
153
    i = 0;
154
    for (Standard_Integer index = 1; index <= numNodes; ++index) {
155
        gp_Pnt _node = triangulation->Node(index);
156
        this->xyz_nodes.row(i) << _node.X(), _node.Y(), _node.Z();
157
        i++;
158
    }
159

160
    this->tris.resize(numTriangles, 3);
161
    i = 0;
162
    for (Standard_Integer index = 1; index <= numTriangles; ++index) {
163
        int n1, n2, n3;
164
        const Poly_Triangle& _tri = triangulation->Triangle(index);
165
        _tri.Get(n1, n2, n3);
166
        this->tris.row(i) << n1 - 1, n2 - 1, n3 - 1;
167
        i++;
168
    }
169
}
170

171
void FaceUnwrapper::findFlatNodes(int steps, double val)
172
{
173
    std::vector<long> fixed_pins;  // TODO: INPUT
174
    lscmrelax::LscmRelax mesh_flattener(this->xyz_nodes.transpose(),
175
                                        this->tris.transpose(),
176
                                        fixed_pins);
177
    mesh_flattener.lscm();
178
    for (int j = 0; j < steps; j++) {
179
        mesh_flattener.relax(val);
180
    }
181
    this->ze_nodes = mesh_flattener.flat_vertices.transpose();
182
}
183

184
ColMat<double, 3> FaceUnwrapper::interpolateFlatFace(const TopoDS_Face& face)
185
{
186
    if (this->uv_nodes.size() == 0) {
187
        throw(std::runtime_error("no uv-coordinates found, interpolating with nurbs is only "
188
                                 "possible if the Flattener was constructed with a nurbs."));
189
    }
190

191
    // extract xyz poles, knots, weights, degree
192
    const Handle(Geom_Surface)& _surface = BRep_Tool::Surface(face);
193
    const Handle(Geom_BSplineSurface)& _bspline = Handle(Geom_BSplineSurface)::DownCast(_surface);
194

195
    const TColStd_Array1OfReal& _uknots = _bspline->UKnotSequence();
196
    const TColStd_Array1OfReal& _vknots = _bspline->VKnotSequence();
197

198
    Eigen::VectorXd weights;
199
    weights.resize(_bspline->NbUPoles() * _bspline->NbVPoles());
200
    long i = 0;
201
    for (long u = 1; u <= _bspline->NbUPoles(); u++) {
202
        for (long v = 1; v <= _bspline->NbVPoles(); v++) {
203
            weights[i] = _bspline->Weight(u, v);
204
            i++;
205
        }
206
    }
207

208
    Eigen::VectorXd u_knots;
209
    Eigen::VectorXd v_knots;
210
    u_knots.resize(_uknots.Length());
211
    v_knots.resize(_vknots.Length());
212
    for (long u = 1; u <= _uknots.Length(); u++) {
213
        u_knots[u - 1] = _uknots.Value(u);
214
    }
215
    for (long v = 1; v <= _vknots.Length(); v++) {
216
        v_knots[v - 1] = _vknots.Value(v);
217
    }
218

219

220
    nu = nurbs::NurbsBase2D(u_knots, v_knots, weights, _bspline->UDegree(), _bspline->VDegree());
221
    A = nu.getInfluenceMatrix(this->uv_nodes);
222

223
    Eigen::LeastSquaresConjugateGradient<spMat> solver;
224
    solver.compute(A);
225
    ColMat<double, 2> ze_poles;
226
    ColMat<double, 3> flat_poles;
227
    ze_poles.resize(weights.rows(), 2);
228
    flat_poles.resize(weights.rows(), 3);
229
    flat_poles.setZero();
230
    ze_poles = solver.solve(ze_nodes);
231
    flat_poles.col(0) << ze_poles.col(0);
232
    flat_poles.col(1) << ze_poles.col(1);
233
    return flat_poles;
234
}
235

236

237
FaceUnwrapper::FaceUnwrapper(ColMat<double, int(3)> xyz_nodes, ColMat<long int, int(3)> tris)
238
{
239
    this->tris = tris;
240
    this->xyz_nodes = xyz_nodes;
241
}
242

243
std::vector<ColMat<double, 3>> FaceUnwrapper::getFlatBoundaryNodes()
244
{
245
    if (this->ze_nodes.size() == 0) {
246
        throw(std::runtime_error("flat vertices not xet computed"));
247
    }
248

249
    ColMat<double, 3> flat_vertices;
250
    flat_vertices.resize(this->ze_nodes.rows(), 3);
251
    flat_vertices.setZero();
252
    flat_vertices.col(0) << this->ze_nodes.col(0);
253
    flat_vertices.col(1) << this->ze_nodes.col(1);
254
    return getBoundaries(flat_vertices, this->tris);
255
}
256