FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2012 Imetric 3D GmbH                                    *
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 <algorithm>
26
#endif
27

28
#include <Base/Sequencer.h>
29

30
#include "Grid.h"
31
#include "Iterator.h"
32
#include "Trim.h"
33

34

35
using namespace MeshCore;
36

37
MeshTrimming::MeshTrimming(MeshKernel& rclM,
38
                           const Base::ViewProjMethod* pclProj,
39
                           const Base::Polygon2d& rclPoly)
40
    : myMesh(rclM)
41
    , myProj(pclProj)
42
    , myPoly(rclPoly)
43
{}
44

45
void MeshTrimming::SetInnerOrOuter(TMode tMode)
46
{
47
    switch (tMode) {
48
        case INNER:
49
            myInner = true;
50
            break;
51
        case OUTER:
52
            myInner = false;
53
            break;
54
    }
55
}
56

57
void MeshTrimming::CheckFacets(const MeshFacetGrid& rclGrid,
58
                               std::vector<FacetIndex>& raulFacets) const
59
{
60
    std::vector<FacetIndex>::iterator it;
61
    MeshFacetIterator clIter(myMesh, 0);
62

63
    // cut inner: use grid to accelerate search
64
    if (myInner) {
65
        Base::BoundBox3f clBBox3d;
66
        Base::BoundBox2d clViewBBox, clPolyBBox;
67
        std::vector<FacetIndex> aulAllElements;
68

69
        // BBox of polygon
70
        clPolyBBox = myPoly.CalcBoundBox();
71
        MeshGridIterator clGridIter(rclGrid);
72
        // traverse all BBoxes
73
        for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
74
            clBBox3d = clGridIter.GetBoundBox();
75
            clViewBBox = clBBox3d.ProjectBox(myProj);
76
            if (clViewBBox.Intersect(clPolyBBox)) {
77
                // save all elements in AllElements
78
                clGridIter.GetElements(aulAllElements);
79
            }
80
        }
81

82
        // remove double elements
83
        std::sort(aulAllElements.begin(), aulAllElements.end());
84
        aulAllElements.erase(std::unique(aulAllElements.begin(), aulAllElements.end()),
85
                             aulAllElements.end());
86

87
        Base::SequencerLauncher seq("Check facets for intersection...", aulAllElements.size());
88

89
        for (it = aulAllElements.begin(); it != aulAllElements.end(); ++it) {
90
            MeshGeomFacet clFacet = myMesh.GetFacet(*it);
91
            if (HasIntersection(clFacet)) {
92
                raulFacets.push_back(*it);
93
            }
94
            seq.next();
95
        }
96
    }
97
    // cut outer
98
    else {
99
        Base::SequencerLauncher seq("Check facets for intersection...", myMesh.CountFacets());
100
        for (clIter.Init(); clIter.More(); clIter.Next()) {
101
            if (HasIntersection(*clIter)) {
102
                raulFacets.push_back(clIter.Position());
103
            }
104
            seq.next();
105
        }
106
    }
107
}
108

109
bool MeshTrimming::HasIntersection(const MeshGeomFacet& rclFacet) const
110
{
111
    Base::Polygon2d clPoly;
112
    Base::Line2d clFacLine, clPolyLine;
113
    Base::Vector2d S;
114
    // is corner of facet inside the polygon
115
    for (auto pnt : rclFacet._aclPoints) {
116
        Base::Vector3f clPt2d = myProj->operator()(pnt);
117
        if (myPoly.Contains(Base::Vector2d(clPt2d.x, clPt2d.y)) == myInner) {
118
            return true;
119
        }
120
        else {
121
            clPoly.Add(Base::Vector2d(clPt2d.x, clPt2d.y));
122
        }
123
    }
124

125
    // is corner of polygon inside the facet
126
    for (size_t j = 0; j < myPoly.GetCtVectors(); j++) {
127
        if (clPoly.Contains(myPoly[j])) {
128
            return true;
129
        }
130
    }
131

132
    // check for other intersections
133
    for (size_t j = 0; j < myPoly.GetCtVectors(); j++) {
134
        clPolyLine.clV1 = myPoly[j];
135
        clPolyLine.clV2 = myPoly[(j + 1) % myPoly.GetCtVectors()];
136

137
        for (size_t i = 0; i < 3; i++) {
138
            clFacLine.clV1 = clPoly[i];
139
            clFacLine.clV2 = clPoly[(i + 1) % 3];
140

141
            if (clPolyLine.IntersectAndContain(clFacLine, S)) {
142
                return true;
143
            }
144
        }
145
    }
146

147
    // no intersection
148
    return false;
149
}
150

151
bool MeshTrimming::PolygonContainsCompleteFacet(bool bInner, FacetIndex ulIndex) const
152
{
153
    const MeshFacet& rclFacet = myMesh._aclFacetArray[ulIndex];
154
    for (PointIndex ptIndex : rclFacet._aulPoints) {
155
        const MeshPoint& rclFacPt = myMesh._aclPointArray[ptIndex];
156
        Base::Vector3f clPt = (*myProj)(rclFacPt);
157
        if (myPoly.Contains(Base::Vector2d(clPt.x, clPt.y)) != bInner) {
158
            return false;
159
        }
160
    }
161

162
    return true;
163
}
164

165
bool MeshTrimming::IsPolygonPointInFacet(FacetIndex ulIndex, Base::Vector3f& clPoint)
166
{
167
    Base::Vector2d A, B, C, P;
168
    float u {}, v {}, w {}, fDetPAC {}, fDetPBC {}, fDetPAB {}, fDetABC {};
169
    Base::Polygon2d clFacPoly;
170
    const MeshGeomFacet& rclFacet = myMesh.GetFacet(ulIndex);
171

172
    for (auto pnt : rclFacet._aclPoints) {
173
        Base::Vector3f clPt = (*myProj)(pnt);
174
        clFacPoly.Add(Base::Vector2d(clPt.x, clPt.y));
175
    }
176

177
    A = clFacPoly[0];
178
    B = clFacPoly[1];
179
    C = clFacPoly[2];
180
    fDetABC =
181
        static_cast<float>(A.x * B.y + A.y * C.x + B.x * C.y - (B.y * C.x + A.y * B.x + A.x * C.y));
182

183
    for (size_t j = 0; j < myPoly.GetCtVectors(); j++) {
184
        // facet contains a polygon point -> calculate the corresponding 3d-point
185
        if (clFacPoly.Contains(myPoly[j])) {
186
            P = myPoly[j];
187
            fDetPAC = static_cast<float>(A.x * P.y + A.y * C.x + P.x * C.y
188
                                         - (P.y * C.x + A.y * P.x + A.x * C.y));
189
            fDetPBC = static_cast<float>(P.x * B.y + P.y * C.x + B.x * C.y
190
                                         - (B.y * C.x + P.y * B.x + P.x * C.y));
191
            fDetPAB = static_cast<float>(A.x * B.y + A.y * P.x + B.x * P.y
192
                                         - (B.y * P.x + A.y * B.x + A.x * P.y));
193
            u = fDetPBC / fDetABC;
194
            v = fDetPAC / fDetABC;
195
            w = fDetPAB / fDetABC;
196

197
            // point is on edge or no valid convex combination
198
            if (u == 0.0f || v == 0.0f || w == 0.0f || fabs(u + v + w - 1.0f) >= 0.001f) {
199
                return false;
200
            }
201
            // 3d point
202
            clPoint = u * rclFacet._aclPoints[0] + v * rclFacet._aclPoints[1]
203
                + w * rclFacet._aclPoints[2];
204

205
            return true;
206
        }
207
    }
208

209
    return false;
210
}
211

212
bool MeshTrimming::GetIntersectionPointsOfPolygonAndFacet(
213
    FacetIndex ulIndex,
214
    int& iSide,
215
    std::vector<Base::Vector3f>& raclPoints) const
216
{
217
    MeshGeomFacet clFac(myMesh.GetFacet(ulIndex));
218
    Base::Vector2d S;
219
    Base::Line2d clFacLine, clPolyLine;
220
    int iIntersections = 0;
221
    int iIntsctWithEdge0 = 0, iIntsctWithEdge1 = 0, iIntsctWithEdge2 = 0;
222

223
    // Edge with no intersection
224
    iSide = -1;
225

226
    for (size_t i = 0; i < myPoly.GetCtVectors(); i++) {
227
        // totally only four intersections allowed
228
        if (iIntersections == 4) {
229
            break;
230
        }
231

232
        Base::Vector2d P3(myPoly[i]), P4(myPoly[(i + 1) % myPoly.GetCtVectors()]);
233
        clPolyLine.clV1 = P3;
234
        clPolyLine.clV2 = P4;
235

236
        for (int j = 0; j < 3; j++) {
237
            Base::Vector3f clP1((*myProj)(clFac._aclPoints[j]));
238
            Base::Vector3f clP2((*myProj)(clFac._aclPoints[(j + 1) % 3]));
239
            Base::Vector2d P1(clP1.x, clP1.y);
240
            Base::Vector2d P2(clP2.x, clP2.y);
241
            clFacLine.clV1 = P1;
242
            clFacLine.clV2 = P2;
243

244
            if (clPolyLine.Intersect(P1, double(MESH_MIN_PT_DIST))) {
245
                // do not pick up corner points
246
                iIntersections++;
247
            }
248
            else if (clPolyLine.Intersect(P2, double(MESH_MIN_PT_DIST))) {
249
                // do not pick up corner points
250
                iIntersections++;
251
            }
252
            else if (clPolyLine.Intersect(clFacLine, S)) {
253
                bool bPushBack = true;
254
                float fP1P2 = static_cast<float>((P2 - P1).Length());
255
                float fSP1 = static_cast<float>((P1 - S).Length());
256
                float fSP2 = static_cast<float>((P2 - S).Length());
257

258
                float fP3P4 = static_cast<float>((P4 - P3).Length());
259
                float fSP3 = static_cast<float>((P3 - S).Length());
260
                float fSP4 = static_cast<float>((P4 - S).Length());
261
                // compute proportion of length
262
                float l = fSP1 / fP1P2;
263
                float m = fSP2 / fP1P2;
264

265
                float r = fSP3 / fP3P4;
266
                float s = fSP4 / fP3P4;
267

268
                // is intersection point convex combination?
269
                if ((fabs(l + m - 1.0f) < 0.001f) && (fabs(r + s - 1.0f) < 0.001f)) {
270
                    Base::Vector3f clIntersection(m * clFac._aclPoints[j]
271
                                                  + l * clFac._aclPoints[(j + 1) % 3]);
272

273
                    iIntersections++;
274

275
                    // only two intersections points per edge allowed
276
                    if (j == 0) {
277
                        if (iIntsctWithEdge0 == 2) {
278
                            bPushBack = false;
279
                        }
280
                        else {
281
                            iIntsctWithEdge0++;
282
                        }
283
                    }
284
                    else if (j == 1) {
285
                        if (iIntsctWithEdge1 == 2) {
286
                            bPushBack = false;
287
                        }
288
                        else {
289
                            iIntsctWithEdge1++;
290
                        }
291
                    }
292
                    else {
293
                        if (iIntsctWithEdge2 == 2) {
294
                            bPushBack = false;
295
                        }
296
                        else {
297
                            iIntsctWithEdge2++;
298
                        }
299
                    }
300

301
                    if (bPushBack) {
302
                        raclPoints.push_back(clIntersection);
303
                    }
304
                }
305
            }
306
        }
307
    }
308

309
    // check for rotating facet later
310
    if (iIntsctWithEdge0 == 0) {
311
        iSide = 0;
312
    }
313
    else if (iIntsctWithEdge1 == 0) {
314
        iSide = 1;
315
    }
316
    else if (iIntsctWithEdge2 == 0) {
317
        iSide = 2;
318
    }
319

320
    // further check (for rotating the facet)
321
    if (iIntsctWithEdge0 == 0 && iIntsctWithEdge1 == 0) {
322
        iSide = 1;
323
    }
324
    else if (iIntsctWithEdge0 == 0 && iIntsctWithEdge2 == 0) {
325
        iSide = 0;
326
    }
327
    else if (iIntsctWithEdge1 == 0 && iIntsctWithEdge2 == 0) {
328
        iSide = 2;
329
    }
330

331
    // and last another check
332
    if (iIntsctWithEdge0 * iIntsctWithEdge1 * iIntsctWithEdge2 > 0) {
333
        if (iIntsctWithEdge0 == 2) {
334
            iSide = 2;
335
        }
336
        else if (iIntsctWithEdge1 == 2) {
337
            iSide = 0;
338
        }
339
        else if (iIntsctWithEdge2 == 2) {
340
            iSide = 1;
341
        }
342
    }
343

344
    return iIntersections > 0;
345
}
346

347
void MeshTrimming::AdjustFacet(MeshFacet& facet, int iInd)
348
{
349
    unsigned long tmp {};
350

351
    if (iInd == 1) {
352
        tmp = facet._aulPoints[0];
353
        facet._aulPoints[0] = facet._aulPoints[1];
354
        facet._aulPoints[1] = facet._aulPoints[2];
355
        facet._aulPoints[2] = tmp;
356
        tmp = facet._aulNeighbours[0];
357
        facet._aulNeighbours[0] = facet._aulNeighbours[1];
358
        facet._aulNeighbours[1] = facet._aulNeighbours[2];
359
        facet._aulNeighbours[2] = tmp;
360
    }
361
    else if (iInd == 2) {
362
        tmp = facet._aulPoints[0];
363
        facet._aulPoints[0] = facet._aulPoints[2];
364
        facet._aulPoints[2] = facet._aulPoints[1];
365
        facet._aulPoints[1] = tmp;
366
        tmp = facet._aulNeighbours[0];
367
        facet._aulNeighbours[0] = facet._aulNeighbours[2];
368
        facet._aulNeighbours[2] = facet._aulNeighbours[1];
369
        facet._aulNeighbours[1] = tmp;
370
    }
371
}
372

373
bool MeshTrimming::CreateFacets(FacetIndex ulFacetPos,
374
                                int iSide,
375
                                const std::vector<Base::Vector3f>& raclPoints,
376
                                std::vector<MeshGeomFacet>& aclNewFacets)
377
{
378
    MeshGeomFacet clFac;
379

380
    // no valid triangulation possible
381
    if (iSide == -1) {
382
        return false;
383
    }
384

385
    // no intersection point found => triangle is only touched at a corner point
386
    if (raclPoints.empty()) {
387
        MeshFacet& facet = myMesh._aclFacetArray[ulFacetPos];
388
        int iCtPtsIn = 0;
389
        int iCtPtsOn = 0;
390
        Base::Vector3f clFacPnt;
391
        Base::Vector2d clProjPnt;
392
        for (PointIndex ptIndex : facet._aulPoints) {
393
            clFacPnt = (*myProj)(myMesh._aclPointArray[ptIndex]);
394
            clProjPnt = Base::Vector2d(clFacPnt.x, clFacPnt.y);
395
            if (myPoly.Intersect(clProjPnt, double(MESH_MIN_PT_DIST))) {
396
                ++iCtPtsOn;
397
            }
398
            else if (myPoly.Contains(clProjPnt) == myInner) {
399
                ++iCtPtsIn;
400
            }
401
        }
402

403
        // in this case we can use the original triangle
404
        if (iCtPtsIn != (3 - iCtPtsOn)) {
405
            aclNewFacets.push_back(myMesh.GetFacet(ulFacetPos));
406
        }
407
    }
408
    // one intersection point found => triangle is also touched at a corner point
409
    else if (raclPoints.size() == 1) {
410
        Base::Vector3f clP(raclPoints[0]);
411
        clP = ((*myProj)(clP));
412
        Base::Vector2d P(clP.x, clP.y);
413
        MeshGeomFacet clFac(myMesh.GetFacet(ulFacetPos));
414

415
        // determine the edge containing the intersection point
416
        Base::Line2d clFacLine;
417
        for (int j = 0; j < 3; j++) {
418
            Base::Vector3f clP1((*myProj)(clFac._aclPoints[j]));
419
            Base::Vector3f clP2((*myProj)(clFac._aclPoints[(j + 1) % 3]));
420
            Base::Vector2d P1(clP1.x, clP1.y);
421
            Base::Vector2d P2(clP2.x, clP2.y);
422
            clFacLine.clV1 = P1;
423
            clFacLine.clV2 = P2;
424

425
            if (clFacLine.Intersect(P, double(MESH_MIN_PT_DIST))) {
426
                if (myPoly.Contains(P1) == myInner) {
427
                    MeshGeomFacet clNew;
428
                    clNew._aclPoints[0] = raclPoints[0];
429
                    clNew._aclPoints[1] = clFac._aclPoints[(j + 1) % 3];
430
                    clNew._aclPoints[2] = clFac._aclPoints[(j + 2) % 3];
431
                    aclNewFacets.push_back(clNew);
432
                    break;
433
                }
434
                else if (myPoly.Contains(P2) == myInner) {
435
                    MeshGeomFacet clNew;
436
                    clNew._aclPoints[0] = raclPoints[0];
437
                    clNew._aclPoints[1] = clFac._aclPoints[(j + 2) % 3];
438
                    clNew._aclPoints[2] = clFac._aclPoints[j];
439
                    aclNewFacets.push_back(clNew);
440
                    break;
441
                }
442
            }
443
        }
444
    }
445
    // two intersection points found
446
    else if (raclPoints.size() == 2) {
447
        MeshFacet& facet = myMesh._aclFacetArray[ulFacetPos];
448
        AdjustFacet(facet, iSide);
449
        Base::Vector3f clP1(raclPoints[0]), clP2(raclPoints[1]);
450

451
        if (iSide == 1) {
452
            // swap P1 and P2
453
            clP1 = raclPoints[1];
454
            clP2 = raclPoints[0];
455
        }
456

457
        // check which facets can be inserted
458
        int iCtPts = 0;
459
        Base::Vector3f clFacPnt;
460
        Base::Vector2d clProjPnt;
461
        for (PointIndex ptIndex : facet._aulPoints) {
462
            clFacPnt = (*myProj)(myMesh._aclPointArray[ptIndex]);
463
            clProjPnt = Base::Vector2d(clFacPnt.x, clFacPnt.y);
464
            if (myPoly.Contains(clProjPnt) == myInner) {
465
                ++iCtPts;
466
            }
467
        }
468

469
        if (iCtPts == 2) {
470
            // erstes Dreieck
471
            clFac._aclPoints[0] = clP1;
472
            clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[2]];
473
            clFac._aclPoints[2] = clP2;
474
            aclNewFacets.push_back(clFac);
475
        }
476
        else if (iCtPts == 1) {
477
            // erstes Dreieck
478
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
479
            clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[1]];
480
            clFac._aclPoints[2] = clP2;
481
            aclNewFacets.push_back(clFac);
482
            // zweites Dreieck
483
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[1]];
484
            clFac._aclPoints[1] = clP1;
485
            clFac._aclPoints[2] = clP2;
486
            aclNewFacets.push_back(clFac);
487
        }
488
        else {
489
            return false;
490
        }
491
    }
492
    // four intersection points found
493
    else if (raclPoints.size() == 4) {
494
        MeshFacet& facet = myMesh._aclFacetArray[ulFacetPos];
495
        AdjustFacet(facet, iSide);
496

497
        clFac = myMesh.GetFacet(ulFacetPos);
498
        // intersection points
499
        Base::Vector3f clP1(raclPoints[0]), clP2(raclPoints[1]), clP3(raclPoints[2]),
500
            clP4(raclPoints[3]);
501

502
        // check which facets can be inserted
503
        int iCtPts = 0;
504
        Base::Vector3f clFacPnt;
505
        Base::Vector2d clProjPnt;
506
        for (PointIndex ptIndex : facet._aulPoints) {
507
            clFacPnt = (*myProj)(myMesh._aclPointArray[ptIndex]);
508
            clProjPnt = Base::Vector2d(clFacPnt.x, clFacPnt.y);
509
            if (myPoly.Contains(clProjPnt) == myInner) {
510
                ++iCtPts;
511
            }
512
        }
513

514
        // sort the intersection points in a certain order
515
        if (iCtPts == 0 || iCtPts == 3) {
516
            if (iSide == 1) {
517
                // swap the points
518
                clP1 = clP2;
519
                clP2 = raclPoints[0];
520
                clP3 = clP4;
521
                clP4 = raclPoints[2];
522
            }
523

524
            if ((clP1 - clFac._aclPoints[1]).Length() > (clP3 - clFac._aclPoints[1]).Length()) {
525
                // swap P1 and P3
526
                Base::Vector3f tmp(clP1);
527
                clP1 = clP3;
528
                clP3 = tmp;
529
            }
530
            if ((clP2 - clFac._aclPoints[0]).Length() > (clP4 - clFac._aclPoints[0]).Length()) {
531
                // swap P2 and P4
532
                Base::Vector3f tmp(clP2);
533
                clP2 = clP4;
534
                clP4 = tmp;
535
            }
536
        }
537
        else {
538
            if (iSide == 0) {
539
                Base::Vector3f clNormal(clFac.GetNormal());
540
                MeshGeomFacet clTmpFac;
541
                clTmpFac._aclPoints[0] = clFac._aclPoints[1];
542
                clTmpFac._aclPoints[1] = clP2;
543
                clTmpFac._aclPoints[2] = clP1;
544
                if (clTmpFac.GetNormal() * clNormal > 0) {
545
                    Base::Vector3f tmp(clP1);
546
                    clP1 = clP2;
547
                    clP2 = tmp;
548
                }
549
                else {
550
                    Base::Vector3f tmp(clP1);
551
                    clP1 = clP4;
552
                    clP4 = clP2;
553
                    clP2 = clP3;
554
                    clP3 = tmp;
555
                }
556
            }
557
            else if (iSide == 1) {
558
                if ((clP2 - clFac._aclPoints[1]).Length() > (clP4 - clFac._aclPoints[1]).Length()) {
559
                    Base::Vector3f tmp(clP1);
560
                    clP1 = clP4;
561
                    clP4 = tmp;
562
                    tmp = clP2;
563
                    clP2 = clP3;
564
                    clP3 = tmp;
565
                }
566
                else {
567
                    Base::Vector3f tmp(clP1);
568
                    clP1 = clP2;
569
                    clP2 = tmp;
570
                    tmp = clP3;
571
                    clP3 = clP4;
572
                    clP4 = tmp;
573
                }
574
            }
575
            else {
576
                if ((clP1 - clFac._aclPoints[1]).Length() > (clP3 - clFac._aclPoints[1]).Length()) {
577
                    Base::Vector3f tmp(clP1);
578
                    clP1 = clP3;
579
                    clP3 = tmp;
580
                    tmp = clP2;
581
                    clP2 = clP4;
582
                    clP4 = tmp;
583
                }
584
            }
585
        }
586

587
        // now create the new facets
588
        if (iCtPts == 0) {
589
            // insert first facet
590
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
591
            clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[1]];
592
            clFac._aclPoints[2] = clP1;
593
            aclNewFacets.push_back(clFac);
594
            // insert second facet
595
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
596
            clFac._aclPoints[1] = clP1;
597
            clFac._aclPoints[2] = clP2;
598
            aclNewFacets.push_back(clFac);
599
            // finally insert third facet
600
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[2]];
601
            clFac._aclPoints[1] = clP4;
602
            clFac._aclPoints[2] = clP3;
603
            aclNewFacets.push_back(clFac);
604
        }
605
        else if (iCtPts == 1) {
606
            // insert first facet
607
            clFac._aclPoints[0] = clP1;
608
            clFac._aclPoints[1] = clP2;
609
            clFac._aclPoints[2] = myMesh._aclPointArray[facet._aulPoints[1]];
610
            aclNewFacets.push_back(clFac);
611
            // finally insert second facet
612
            clFac._aclPoints[0] = clP4;
613
            clFac._aclPoints[1] = clP3;
614
            clFac._aclPoints[2] = myMesh._aclPointArray[facet._aulPoints[2]];
615
            aclNewFacets.push_back(clFac);
616
        }
617
        else if (iCtPts == 2) {
618
            // insert first facet
619
            clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
620
            clFac._aclPoints[1] = clP2;
621
            clFac._aclPoints[2] = clP4;
622
            aclNewFacets.push_back(clFac);
623
            // insert second facet
624
            clFac._aclPoints[0] = clP1;
625
            clFac._aclPoints[1] = clP4;
626
            clFac._aclPoints[2] = clP2;
627
            aclNewFacets.push_back(clFac);
628
            // finally insert third facet
629
            clFac._aclPoints[0] = clP1;
630
            clFac._aclPoints[1] = clP3;
631
            clFac._aclPoints[2] = clP4;
632
            aclNewFacets.push_back(clFac);
633
        }
634
        else {
635
            // insert first facet
636
            clFac._aclPoints[0] = clP1;
637
            clFac._aclPoints[1] = clP3;
638
            clFac._aclPoints[2] = clP4;
639
            aclNewFacets.push_back(clFac);
640
            // finally insert second facet
641
            clFac._aclPoints[0] = clP1;
642
            clFac._aclPoints[1] = clP4;
643
            clFac._aclPoints[2] = clP2;
644
            aclNewFacets.push_back(clFac);
645
        }
646
    }
647
    else {
648
        return false;
649
    }
650

651
    return true;
652
}
653

654
bool MeshTrimming::CreateFacets(FacetIndex ulFacetPos,
655
                                int iSide,
656
                                const std::vector<Base::Vector3f>& raclPoints,
657
                                Base::Vector3f& clP3,
658
                                std::vector<MeshGeomFacet>& aclNewFacets)
659
{
660
    // no valid triangulation possible
661
    if (iSide == -1 || raclPoints.size() < 2) {
662
        return false;
663
    }
664

665
    Base::Vector3f clP1(raclPoints[0]);
666
    Base::Vector3f clP2(raclPoints[1]);
667

668
    MeshFacet& facet = myMesh._aclFacetArray[ulFacetPos];
669
    AdjustFacet(facet, iSide);
670

671
    MeshGeomFacet clFac;
672

673
    float fDistEdgeP1 = clP1.DistanceToLineSegment(myMesh._aclPointArray[facet._aulPoints[1]],
674
                                                   myMesh._aclPointArray[facet._aulPoints[2]])
675
                            .Length();
676
    float fDistEdgeP2 = clP2.DistanceToLineSegment(myMesh._aclPointArray[facet._aulPoints[1]],
677
                                                   myMesh._aclPointArray[facet._aulPoints[2]])
678
                            .Length();
679

680
    // swap P1 and P2
681
    if (fDistEdgeP2 < fDistEdgeP1) {
682
        Base::Vector3f tmp(clP1);
683
        clP1 = clP2;
684
        clP2 = tmp;
685
    }
686

687
    // check which facets should be inserted
688
    int iCtPts = 0;
689
    Base::Vector3f clFacPnt;
690
    Base::Vector2d clProjPnt;
691
    for (PointIndex ptIndex : facet._aulPoints) {
692
        clFacPnt = (*myProj)(myMesh._aclPointArray[ptIndex]);
693
        clProjPnt = Base::Vector2d(clFacPnt.x, clFacPnt.y);
694
        if (myPoly.Contains(clProjPnt) == myInner) {
695
            ++iCtPts;
696
        }
697
    }
698
    if (iCtPts == 3) {
699
        clFac = myMesh.GetFacet(ulFacetPos);
700
        if ((clP1 - clFac._aclPoints[1]).Length() > (clP2 - clFac._aclPoints[1]).Length()) {
701
            Base::Vector3f tmp(clP1);
702
            clP1 = clP2;
703
            clP2 = tmp;
704
        }
705
        // only one facet
706
        clFac._aclPoints[0] = clP1;
707
        clFac._aclPoints[1] = clP2;
708
        clFac._aclPoints[2] = clP3;
709
        aclNewFacets.push_back(clFac);
710
    }
711
    else if (iCtPts == 2) {
712
        // first facet
713
        clFac._aclPoints[0] = clP1;
714
        clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[2]];
715
        clFac._aclPoints[2] = clP3;
716
        aclNewFacets.push_back(clFac);
717
        // second facet
718
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[2]];
719
        clFac._aclPoints[1] = clP2;
720
        clFac._aclPoints[2] = clP3;
721
        aclNewFacets.push_back(clFac);
722
    }
723
    else if (iCtPts == 1) {
724
        // first facet
725
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
726
        clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[1]];
727
        clFac._aclPoints[2] = clP3;
728
        aclNewFacets.push_back(clFac);
729
        // second facet
730
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[1]];
731
        clFac._aclPoints[1] = clP1;
732
        clFac._aclPoints[2] = clP3;
733
        aclNewFacets.push_back(clFac);
734
        // third facet
735
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
736
        clFac._aclPoints[1] = clP3;
737
        clFac._aclPoints[2] = clP2;
738
        aclNewFacets.push_back(clFac);
739
    }
740
    else if (iCtPts == 0) {
741
        clFac = myMesh.GetFacet(ulFacetPos);
742
        if ((clP1 - clFac._aclPoints[1]).Length() > (clP2 - clFac._aclPoints[1]).Length()) {
743
            Base::Vector3f tmp(clP1);
744
            clP1 = clP2;
745
            clP2 = tmp;
746
        }
747
        // first facet
748
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[2]];
749
        clFac._aclPoints[1] = clP3;
750
        clFac._aclPoints[2] = clP2;
751
        aclNewFacets.push_back(clFac);
752
        // second facet
753
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[2]];
754
        clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[0]];
755
        clFac._aclPoints[2] = clP3;
756
        aclNewFacets.push_back(clFac);
757
        // third facet
758
        clFac._aclPoints[0] = myMesh._aclPointArray[facet._aulPoints[0]];
759
        clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[1]];
760
        clFac._aclPoints[2] = clP3;
761
        aclNewFacets.push_back(clFac);
762
        // and finally fourth facet
763
        clFac._aclPoints[0] = clP3;
764
        clFac._aclPoints[1] = myMesh._aclPointArray[facet._aulPoints[1]];
765
        clFac._aclPoints[2] = clP1;
766
        aclNewFacets.push_back(clFac);
767
    }
768

769
    return true;
770
}
771

772
void MeshTrimming::TrimFacets(const std::vector<FacetIndex>& raulFacets,
773
                              std::vector<MeshGeomFacet>& aclNewFacets)
774
{
775
    Base::Vector3f clP;
776
    std::vector<Base::Vector3f> clIntsct;
777
    int iSide {};
778

779
    Base::SequencerLauncher seq("trimming facets...", raulFacets.size());
780
    for (FacetIndex index : raulFacets) {
781
        clIntsct.clear();
782
        if (!IsPolygonPointInFacet(index, clP)) {
783
            // facet must be trimmed
784
            if (!PolygonContainsCompleteFacet(myInner, index)) {
785
                // generate new facets
786
                if (GetIntersectionPointsOfPolygonAndFacet(index, iSide, clIntsct)) {
787
                    CreateFacets(index, iSide, clIntsct, myTriangles);
788
                }
789
            }
790
        }
791
        // facet contains a polygon point
792
        else {
793
            // generate new facets
794
            if (GetIntersectionPointsOfPolygonAndFacet(index, iSide, clIntsct)) {
795
                CreateFacets(index, iSide, clIntsct, clP, myTriangles);
796
            }
797
        }
798
        seq.next();
799
    }
800

801
    aclNewFacets = myTriangles;
802
}
803