FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2008 Juergen Riegel <juergen.riegel@web.de>             *
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
#ifdef FC_OS_LINUX
26
#include <unistd.h>
27
#endif
28
#endif
29

30
#include <Base/Builder3D.h>
31
#include <Base/Console.h>
32
#include <Mod/Mesh/App/Core/Evaluation.h>
33
#include <Mod/Mesh/App/Core/Iterator.h>
34
#include <Mod/Mesh/App/Core/MeshKernel.h>
35
#include <Mod/Mesh/App/Core/TopoAlgorithm.h>
36

37
#include "MeshAlgos.h"
38

39

40
using namespace MeshPart;
41
using namespace MeshCore;
42

43

44
void MeshAlgos::offset(MeshCore::MeshKernel* Mesh, float fSize)
45
{
46
    std::vector<Base::Vector3f> normals = Mesh->CalcVertexNormals();
47

48
    unsigned int i = 0;
49
    // go through all the Vertex normals
50
    for (std::vector<Base::Vector3f>::iterator It = normals.begin(); It != normals.end();
51
         ++It, i++) {
52
        // and move each mesh point in the normal direction
53
        Mesh->MovePoint(i, It->Normalize() * fSize);
54
    }
55
    Mesh->RecalcBoundBox();
56
}
57

58

59
void MeshAlgos::offsetSpecial2(MeshCore::MeshKernel* Mesh, float fSize)
60
{
61
    Base::Builder3D builder;
62
    std::vector<Base::Vector3f> PointNormals = Mesh->CalcVertexNormals();
63
    std::vector<Base::Vector3f> FaceNormals;
64
    std::set<MeshCore::FacetIndex> fliped;
65

66
    MeshFacetIterator it(*Mesh);
67
    for (it.Init(); it.More(); it.Next()) {
68
        FaceNormals.push_back(it->GetNormal().Normalize());
69
    }
70

71
    unsigned int i = 0;
72

73
    // go through all the Vertex normals
74
    for (std::vector<Base::Vector3f>::iterator It = PointNormals.begin(); It != PointNormals.end();
75
         ++It, i++) {
76
        Base::Line3f line {Mesh->GetPoint(i), Mesh->GetPoint(i) + It->Normalize() * fSize};
77
        Base::DrawStyle drawStyle;
78
        builder.addNode(Base::LineItem {line, drawStyle});
79
        // and move each mesh point in the normal direction
80
        Mesh->MovePoint(i, It->Normalize() * fSize);
81
    }
82
    Mesh->RecalcBoundBox();
83

84
    MeshTopoAlgorithm alg(*Mesh);
85

86
    for (int l = 0; l < 1; l++) {
87
        for (it.Init(), i = 0; it.More(); it.Next(), i++) {
88
            if (it->IsFlag(MeshFacet::INVALID)) {
89
                continue;
90
            }
91
            // calculate the angle between them
92
            float angle = acos((FaceNormals[i] * it->GetNormal())
93
                               / (it->GetNormal().Length() * FaceNormals[i].Length()));
94
            if (angle > 1.6) {
95
                Base::DrawStyle drawStyle;
96
                drawStyle.pointSize = 4.0F;
97
                Base::PointItem item {it->GetGravityPoint(),
98
                                      drawStyle,
99
                                      Base::ColorRGB {1.0F, 0.0F, 0.0F}};
100
                builder.addNode(item);
101
                fliped.insert(it.Position());
102
            }
103
        }
104

105
        // if there are no flipped triangles -> stop
106
        // int f =fliped.size();
107
        if (fliped.empty()) {
108
            break;
109
        }
110

111
        for (MeshCore::FacetIndex It : fliped) {
112
            alg.CollapseFacet(It);
113
        }
114
        fliped.clear();
115
    }
116

117
    alg.Cleanup();
118

119
    // search for intersected facets
120
    MeshCore::MeshEvalSelfIntersection eval(*Mesh);
121
    std::vector<std::pair<MeshCore::FacetIndex, MeshCore::FacetIndex>> faces;
122
    eval.GetIntersections(faces);
123

124

125
    builder.saveToLog();
126
}
127

128
void MeshAlgos::offsetSpecial(MeshCore::MeshKernel* Mesh, float fSize, float zmax, float zmin)
129
{
130
    std::vector<Base::Vector3f> normals = Mesh->CalcVertexNormals();
131

132
    unsigned int i = 0;
133
    // go through all the Vertex normals
134
    for (std::vector<Base::Vector3f>::iterator It = normals.begin(); It != normals.end();
135
         ++It, i++) {
136
        Base::Vector3f Pnt = Mesh->GetPoint(i);
137

138
        if (Pnt.z < zmax && Pnt.z > zmin) {
139
            Pnt.z = 0;
140
            Mesh->MovePoint(i, Pnt.Normalize() * fSize);
141
        }
142
        else {
143
            // and move each mesh point in the normal direction
144
            Mesh->MovePoint(i, It->Normalize() * fSize);
145
        }
146
    }
147
}
148

149

150
void MeshAlgos::coarsen(MeshCore::MeshKernel* /*Mesh*/, float /*f*/)
151
{
152
#ifdef FC_USE_GTS
153
    GtsSurface* surface;
154

155
    // create a GTS surface
156
    surface = MeshAlgos::createGTSSurface(Mesh);
157

158
    Mesh->Clear();
159

160
    guint stop_number = 100000;
161
    gdouble fold = 3.1415 / 180.;
162

163
    gts_surface_coarsen(surface,
164
                        NULL,
165
                        NULL,
166
                        NULL,
167
                        NULL,
168
                        (GtsStopFunc)gts_coarsen_stop_number,
169
                        &stop_number,
170
                        fold);
171

172
    // get the standard mesh
173
    fillMeshFromGTSSurface(Mesh, surface);
174
#endif
175
}
176

177

178
MeshCore::MeshKernel* MeshAlgos::boolean(MeshCore::MeshKernel* pMesh1,
179
                                         MeshCore::MeshKernel* /*pMesh2*/,
180
                                         MeshCore::MeshKernel* /*pResult*/,
181
                                         int /*Type*/)
182
{
183
#ifdef FC_USE_GTS
184
    GtsSurface *s1, *s2, *s3;
185
    GtsSurfaceInter* si;
186
    GNode *tree1, *tree2;
187
    gboolean check_self_intersection = false;
188
    gboolean closed = true, is_open1, is_open2;
189

190

191
    // create a GTS surface
192
    s1 = MeshAlgos::createGTSSurface(pMesh1);
193
    s2 = MeshAlgos::createGTSSurface(pMesh2);
194

195
    // clear the mesh (memory)
196
    // Mesh1.clear();
197
    // Mesh2.clear();
198

199
    /* check that the surfaces are orientable manifolds */
200
    if (!gts_surface_is_orientable(s1)) {
201
        gts_object_destroy(GTS_OBJECT(s1));
202
        gts_object_destroy(GTS_OBJECT(s2));
203
        throw std::runtime_error("surface 1 is not an orientable manifold\n");
204
    }
205
    if (!gts_surface_is_orientable(s2)) {
206
        gts_object_destroy(GTS_OBJECT(s1));
207
        gts_object_destroy(GTS_OBJECT(s2));
208
        throw std::runtime_error("surface 2 is not an orientable manifold\n");
209
    }
210

211
    /* check that the surfaces are not self-intersecting */
212
    if (check_self_intersection) {
213
        GtsSurface* self_intersects;
214

215
        self_intersects = gts_surface_is_self_intersecting(s1);
216
        if (self_intersects != NULL) {
217
            //      if (verbose)
218
            //        gts_surface_print_stats (self_intersects, stderr);
219
            //      gts_surface_write (self_intersects, stdout);
220
            gts_object_destroy(GTS_OBJECT(self_intersects));
221
            gts_object_destroy(GTS_OBJECT(s1));
222
            gts_object_destroy(GTS_OBJECT(s2));
223
            throw std::runtime_error("surface is self-intersecting\n");
224
        }
225
        self_intersects = gts_surface_is_self_intersecting(s2);
226
        if (self_intersects != NULL) {
227
            //      if (verbose)
228
            //        gts_surface_print_stats (self_intersects, stderr);
229
            //      gts_surface_write (self_intersects, stdout);
230
            gts_object_destroy(GTS_OBJECT(self_intersects));
231
            gts_object_destroy(GTS_OBJECT(s1));
232
            gts_object_destroy(GTS_OBJECT(s2));
233
            throw std::runtime_error("surface is self-intersecting\n");
234
        }
235
    }
236

237
    /* build bounding box tree for first surface */
238
    tree1 = gts_bb_tree_surface(s1);
239
    is_open1 = gts_surface_volume(s1) < 0. ? true : false;
240

241
    /* build bounding box tree for second surface */
242
    tree2 = gts_bb_tree_surface(s2);
243
    is_open2 = gts_surface_volume(s2) < 0. ? true : false;
244

245
    si = gts_surface_inter_new(gts_surface_inter_class(), s1, s2, tree1, tree2, is_open1, is_open2);
246
    g_assert(gts_surface_inter_check(si, &closed));
247
    if (!closed) {
248
        gts_object_destroy(GTS_OBJECT(s1));
249
        gts_object_destroy(GTS_OBJECT(s2));
250
        gts_bb_tree_destroy(tree1, true);
251
        gts_bb_tree_destroy(tree2, true);
252
        throw "the intersection of 1 and  2 is not a closed curve\n";
253
    }
254

255
    s3 = gts_surface_new(gts_surface_class(),
256
                         gts_face_class(),
257
                         gts_edge_class(),
258
                         gts_vertex_class());
259
    if (Type == 0) {  // union
260
        gts_surface_inter_boolean(si, s3, GTS_1_OUT_2);
261
        gts_surface_inter_boolean(si, s3, GTS_2_OUT_1);
262
    }
263
    else if (Type == 1) {  // inter
264
        gts_surface_inter_boolean(si, s3, GTS_1_IN_2);
265
        gts_surface_inter_boolean(si, s3, GTS_2_IN_1);
266
    }
267
    else if (Type == 2) {  // diff
268
        gts_surface_inter_boolean(si, s3, GTS_1_OUT_2);
269
        gts_surface_inter_boolean(si, s3, GTS_2_IN_1);
270
        gts_surface_foreach_face(si->s2, (GtsFunc)gts_triangle_revert, NULL);
271
        gts_surface_foreach_face(s2, (GtsFunc)gts_triangle_revert, NULL);
272
    }
273
    else if (Type == 3) {  // cut inner
274
        gts_surface_inter_boolean(si, s3, GTS_1_IN_2);
275
    }
276
    else if (Type == 4) {  // cut outer
277
        gts_surface_inter_boolean(si, s3, GTS_1_OUT_2);
278
    }
279

280
    // check that the resulting surface is not self-intersecting
281
    if (check_self_intersection) {
282
        GtsSurface* self_intersects;
283

284
        self_intersects = gts_surface_is_self_intersecting(s3);
285
        if (self_intersects != NULL) {
286
            //      if (verbose)
287
            //        gts_surface_print_stats (self_intersects, stderr);
288
            //     gts_surface_write (self_intersects, stdout);
289
            gts_object_destroy(GTS_OBJECT(self_intersects));
290
            gts_object_destroy(GTS_OBJECT(s1));
291
            gts_object_destroy(GTS_OBJECT(s2));
292
            gts_object_destroy(GTS_OBJECT(s3));
293
            gts_object_destroy(GTS_OBJECT(si));
294
            gts_bb_tree_destroy(tree1, true);
295
            gts_bb_tree_destroy(tree2, true);
296
            throw std::runtime_error("the resulting surface is self-intersecting\n");
297
        }
298
    }
299
    // display summary information about the resulting surface
300
    //  if (verbose)
301
    //    gts_surface_print_stats (s3, stderr);
302
    // write resulting surface to standard output
303

304
    // get the standard mesh
305
    fillMeshFromGTSSurface(pResult, s3);
306

307

308
    // destroy surfaces
309
    gts_object_destroy(GTS_OBJECT(s1));
310
    gts_object_destroy(GTS_OBJECT(s2));
311
    //  gts_object_destroy (GTS_OBJECT (s3));
312
    //  gts_object_destroy (GTS_OBJECT (si));
313

314
    // destroy bounding box trees (including bounding boxes)
315
    //  gts_bb_tree_destroy (tree1, true);
316
    //  gts_bb_tree_destroy (tree2, true);
317

318
#endif
319
    return pMesh1;
320
}
321

322

323
#ifdef FC_USE_GTS
324

325

326
/// helper function - construct a Edge out of two Vertexes if not already there
327
static GtsEdge* new_edge(GtsVertex* v1, GtsVertex* v2)
328
{
329
    GtsSegment* s = gts_vertices_are_connected(v1, v2);
330
    if (s == NULL) {
331
        return gts_edge_new(gts_edge_class(), v1, v2);
332
    }
333
    else {
334
        return GTS_EDGE(s);
335
    }
336
}
337

338

339
GtsSurface* MeshAlgos::createGTSSurface(MeshCore::MeshKernel* Mesh)
340
{
341
    GtsSurface* Surf = gts_surface_new(gts_surface_class(),
342
                                       gts_face_class(),
343
                                       gts_edge_class(),
344
                                       gts_vertex_class());
345

346
    unsigned long p1, p2, p3;
347
    Base::Vector3f Vertex;
348

349

350
    // Getting all the points
351
    GtsVertex** aVertex = (GtsVertex**)malloc(Mesh->CountPoints() * sizeof(GtsVertex*));
352
    for (unsigned int PIter = 0; PIter < Mesh->CountPoints(); PIter++) {
353
        Vertex = Mesh->GetPoint(PIter);
354
        aVertex[PIter] = gts_vertex_new(gts_vertex_class(), Vertex.x, Vertex.y, Vertex.z);
355
    }
356

357
    // cycling through the facets
358
    for (unsigned int pFIter = 0; pFIter < Mesh->CountFacets(); pFIter++) {
359
        // getting the three points of the facet
360
        Mesh->GetFacetPoints(pFIter, p1, p2, p3);
361

362
        // creating the edges and add the face to the surface
363
        gts_surface_add_face(Surf,
364
                             gts_face_new(Surf->face_class,
365
                                          new_edge(aVertex[p1], aVertex[p2]),
366
                                          new_edge(aVertex[p2], aVertex[p3]),
367
                                          new_edge(aVertex[p3], aVertex[p1])));
368
    }
369

370
    Base::Console().Log("GTS [%d faces, %d Points, %d Edges,%s ,%s]\n",
371
                        gts_surface_face_number(Surf),
372
                        gts_surface_vertex_number(Surf),
373
                        gts_surface_edge_number(Surf),
374
                        gts_surface_is_orientable(Surf) ? "orientable" : "not orientable",
375
                        gts_surface_is_self_intersecting(Surf) ? "self-intersections"
376
                                                               : "no self-intersection");
377

378
    return Surf;
379
}
380

381
/// helper function for the face (triangle iteration
382
static void onFaces(GtsTriangle* t, std::vector<MeshGeomFacet>* VAry)
383
{
384
    GtsVertex *mv0, *mv1, *mv2;
385

386
    gts_triangle_vertices(t, &mv0, &mv1, &mv2);
387

388
    VAry->push_back(MeshGeomFacet(Base::Vector3f(mv0->p.x, mv0->p.y, mv0->p.z),
389
                                  Base::Vector3f(mv1->p.x, mv1->p.y, mv1->p.z),
390
                                  Base::Vector3f(mv2->p.x, mv2->p.y, mv2->p.z)));
391
}
392

393
/*
394
static void onVertices(GtsVertex *v, MeshKernel *pKernel )
395
{
396
  Base::Vector3f Point(GTS_POINT(v)->x,GTS_POINT(v)->y,GTS_POINT(v)->z);
397
}*/
398

399
void MeshAlgos::fillMeshFromGTSSurface(MeshCore::MeshKernel* pMesh, GtsSurface* pSurface)
400
{
401
    std::vector<MeshGeomFacet> VAry;
402

403
    // remove old mesh
404
    pMesh->Clear();
405

406
    //  gts_surface_foreach_vertex(pSurface,(GtsFunc) onVertices,&MeshK);
407
    gts_surface_foreach_face(pSurface, (GtsFunc)onFaces, &VAry);
408

409
    // destroy surfaces
410
    gts_object_destroy(GTS_OBJECT(pSurface));
411

412
    // put the facets the simple way in the mesh, totp is recalculated!
413
    (*pMesh) = VAry;
414
}
415

416
#endif
417

418
#include <BRep_Tool.hxx>
419
#include <GeomAPI_IntCS.hxx>
420
#include <GeomLProp_CLProps.hxx>
421
#include <Geom_Curve.hxx>
422
#include <Geom_Plane.hxx>
423
#include <TopExp.hxx>
424
#include <TopExp_Explorer.hxx>
425
#include <TopoDS.hxx>
426
#include <TopoDS_Edge.hxx>
427
#include <TopoDS_Vertex.hxx>
428
#include <TopoDS_Wire.hxx>
429

430
void MeshAlgos::cutByShape(const TopoDS_Shape& aShape,
431
                           const MeshCore::MeshKernel* pMesh,
432
                           MeshCore::MeshKernel* pToolMesh)
433
{
434

435
    // calculate the projection for each Edge
436
    //  CurveProjectorShape Project(aShape,*pMesh);
437
    CurveProjectorWithToolMesh Project(aShape, *pMesh, *pToolMesh);
438

439
    // IntersectionLine Lines;
440
    //  MeshWithProperty *ResultMesh = new MeshWithProperty();
441

442

443
    // boolean(pMesh,ToolMesh,ResultMesh,1);
444
}
445

446
/*
447
void MeshAlgos::doIntersection(const MeshWithProperty &pMesh,const MeshWithProperty
448
ToolMesh,IntersectionLine &Lines)
449
  {
450

451

452
  }
453

454
*/
455

456
void MeshAlgos::cutByCurve(MeshCore::MeshKernel* pMesh,
457
                           const std::vector<CurveProjector::FaceSplitEdge>& vSplitEdges)
458
{
459
    MeshTopoAlgorithm cTopAlg(*pMesh);
460

461
    for (const auto& it : vSplitEdges) {
462
        cTopAlg.SplitFacet(it.ulFaceIndex, it.p1, it.p2);
463
    }
464
}
465

466
class _VertexCompare
467
{
468
public:
469
    bool operator()(const TopoDS_Vertex& rclV1, const TopoDS_Vertex& rclV2) const
470
    {
471
        if (rclV1.IsSame(rclV2) == Standard_True) {
472
            return false;
473
        }
474

475
        gp_XYZ clP1 = BRep_Tool::Pnt(rclV1).XYZ();
476
        gp_XYZ clP2 = BRep_Tool::Pnt(rclV2).XYZ();
477

478
        if (fabs(clP1.X() - clP2.X()) < dE) {
479
            if (fabs(clP1.Y() - clP2.Y()) < dE) {
480
                return clP1.Z() < clP2.Z();
481
            }
482
            else {
483
                return clP1.Y() < clP2.Y();
484
            }
485
        }
486
        else {
487
            return clP1.X() < clP2.X();
488
        }
489
    }
490

491
    double dE = 1.0e-5;
492
};
493

494

495
void MeshAlgos::LoftOnCurve(MeshCore::MeshKernel& ResultMesh,
496
                            const TopoDS_Shape& Shape,
497
                            const std::vector<Base::Vector3f>& poly,
498
                            const Base::Vector3f& up,
499
                            float MaxSize)
500
{
501
    TopExp_Explorer Ex;
502
    Standard_Real fBegin, fEnd;
503
    std::vector<MeshGeomFacet> cVAry;
504
    std::map<TopoDS_Vertex, std::vector<Base::Vector3f>, _VertexCompare> ConnectMap;
505

506
    for (Ex.Init(Shape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
507
        // get the edge and the belonging Vertexes
508
        TopoDS_Edge Edge = (TopoDS_Edge&)Ex.Current();
509
        TopoDS_Vertex V1, V2;
510
        TopExp::Vertices(Edge, V1, V2);
511
        bool bBegin = false, bEnd = false;
512
        // getting the geometric curve and the interval
513
        GeomLProp_CLProps prop(BRep_Tool::Curve(Edge, fBegin, fEnd), 1, 0.0000000001);
514
        int res = int((fEnd - fBegin) / MaxSize);
515
        // do at least 2 segments
516
        if (res < 2) {
517
            res = 2;
518
        }
519
        gp_Dir Tangent;
520

521
        std::vector<Base::Vector3f> prePoint(poly.size());
522
        std::vector<Base::Vector3f> actPoint(poly.size());
523

524
        // checking if there is already a end to connect
525
        if (ConnectMap.find(V1) != ConnectMap.end()) {
526
            bBegin = true;
527
            prePoint = ConnectMap[V1];
528
        }
529

530
        if (ConnectMap.find(V2) != ConnectMap.end()) {
531
            bEnd = true;
532
        }
533

534
        for (long i = 0; i < res; i++) {
535

536
            // get point and tangent at the position, up is fix for the moment
537
            prop.SetParameter(fBegin + ((fEnd - fBegin) * float(i)) / float(res - 1));
538
            prop.Tangent(Tangent);
539
            Base::Vector3f Tng((float)Tangent.X(), (float)Tangent.Y(), (float)Tangent.Z());
540
            Base::Vector3f Ptn((float)prop.Value().X(),
541
                               (float)prop.Value().Y(),
542
                               (float)prop.Value().Z());
543
            Base::Vector3f Up(up);
544
            // normalize and calc the third vector of the plane coordinatesystem
545
            Tng.Normalize();
546
            Up.Normalize();
547
            Base::Vector3f Third(Tng % Up);
548

549
            //      Base::Console().Log("Pos: %f %f %f \n",Ptn.x,Ptn.y,Ptn.z);
550

551
            unsigned int l = 0;
552
            std::vector<Base::Vector3f>::const_iterator It;
553

554
            // got through the profile
555
            for (It = poly.begin(); It != poly.end(); ++It, l++) {
556
                actPoint[l] = ((Third * It->x) + (Up * It->y) + (Tng * It->z) + Ptn);
557
            }
558

559
            if (i == res - 1 && !bEnd) {
560
                // remember the last row to connect to a otger edge with the same vertex
561
                ConnectMap[V2] = actPoint;
562
            }
563

564
            if (i == 1 && bBegin) {
565
                // using the end of an other edge as start
566
                prePoint = ConnectMap[V1];
567
            }
568

569
            if (i == 0 && !bBegin) {
570
                // remember the first row for connection to a edge with the same vertex
571
                ConnectMap[V1] = actPoint;
572
            }
573

574
            if (i)  // not the first row or something to connect to
575
            {
576
                for (l = 0; l < actPoint.size(); l++) {
577
                    if (l)  // not first point in row
578
                    {
579
                        if (i == res - 1 && bEnd) {  // if last row and a end to connect
580
                            actPoint = ConnectMap[V2];
581
                        }
582

583
                        Base::Vector3f p1 = prePoint[l - 1], p2 = actPoint[l - 1], p3 = prePoint[l],
584
                                       p4 = actPoint[l];
585

586
                        cVAry.emplace_back(p1, p2, p3);
587
                        cVAry.emplace_back(p3, p2, p4);
588
                    }
589
                }
590
            }
591

592
            prePoint = actPoint;
593
        }
594
    }
595

596
    ResultMesh.AddFacets(cVAry);
597
}
598