FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2005 Berthold Grupp                                     *
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

25
#ifndef _PreComp_
26
#include <fstream>
27
#include <ios>
28
#endif
29

30
#include <Base/Builder3D.h>
31
#include <Base/Sequencer.h>
32

33
#include "Algorithm.h"
34
#include "Builder.h"
35
#include "Definitions.h"
36
#include "Elements.h"
37
#include "Grid.h"
38
#include "Iterator.h"
39
#include "SetOperations.h"
40
#include "Triangulation.h"
41
#include "Visitor.h"
42

43

44
using namespace Base;
45
using namespace MeshCore;
46

47

48
SetOperations::SetOperations(const MeshKernel& cutMesh1,
49
                             const MeshKernel& cutMesh2,
50
                             MeshKernel& result,
51
                             OperationType opType,
52
                             float minDistanceToPoint)
53
    : _cutMesh0(cutMesh1)
54
    , _cutMesh1(cutMesh2)
55
    , _resultMesh(result)
56
    , _operationType(opType)
57
    , _minDistanceToPoint(minDistanceToPoint)
58
{}
59

60
void SetOperations::Do()
61
{
62
    _minDistanceToPoint = 0.000001f;
63
    float saveMinMeshDistance = MeshDefinitions::_fMinPointDistance;
64
    MeshDefinitions::SetMinPointDistance(0.000001f);
65

66
    //  Base::Sequencer().start("set operation", 5);
67

68
    // _builder.clear();
69

70
    // Base::Sequencer().next();
71
    std::set<FacetIndex> facetsCuttingEdge0, facetsCuttingEdge1;
72
    Cut(facetsCuttingEdge0, facetsCuttingEdge1);
73

74
    // no intersection curve of the meshes found
75
    if (facetsCuttingEdge0.empty() || facetsCuttingEdge1.empty()) {
76
        switch (_operationType) {
77
            case Union: {
78
                _resultMesh = _cutMesh0;
79
                _resultMesh.Merge(_cutMesh1);
80
            } break;
81
            case Intersect: {
82
                _resultMesh.Clear();
83
            } break;
84
            case Difference:
85
            case Inner:
86
            case Outer: {
87
                _resultMesh = _cutMesh0;
88
            } break;
89
            default: {
90
                _resultMesh.Clear();
91
                break;
92
            }
93
        }
94

95
        MeshDefinitions::SetMinPointDistance(saveMinMeshDistance);
96
        return;
97
    }
98

99
    for (auto i = 0UL; i < _cutMesh0.CountFacets(); i++) {
100
        if (facetsCuttingEdge0.find(i) == facetsCuttingEdge0.end()) {
101
            _newMeshFacets[0].push_back(_cutMesh0.GetFacet(i));
102
        }
103
    }
104

105
    for (auto i = 0UL; i < _cutMesh1.CountFacets(); i++) {
106
        if (facetsCuttingEdge1.find(i) == facetsCuttingEdge1.end()) {
107
            _newMeshFacets[1].push_back(_cutMesh1.GetFacet(i));
108
        }
109
    }
110

111
    // Base::Sequencer().next();
112
    TriangulateMesh(_cutMesh0, 0);
113

114
    // Base::Sequencer().next();
115
    TriangulateMesh(_cutMesh1, 1);
116

117
    float mult0 {}, mult1 {};
118
    switch (_operationType) {
119
        case Union:
120
            mult0 = -1.0f;
121
            mult1 = -1.0f;
122
            break;
123
        case Intersect:
124
            mult0 = 1.0f;
125
            mult1 = 1.0f;
126
            break;
127
        case Difference:
128
            mult0 = -1.0f;
129
            mult1 = 1.0f;
130
            break;
131
        case Inner:
132
            mult0 = 1.0f;
133
            mult1 = 0.0f;
134
            break;
135
        case Outer:
136
            mult0 = -1.0f;
137
            mult1 = 0.0f;
138
            break;
139
        default:
140
            mult0 = 0.0f;
141
            mult1 = 0.0f;
142
            break;
143
    }
144

145
    // Base::Sequencer().next();
146
    CollectFacets(0, mult0);
147
    // Base::Sequencer().next();
148
    CollectFacets(1, mult1);
149

150
    std::vector<MeshGeomFacet> facets;
151

152
    std::vector<MeshGeomFacet>::iterator itf;
153
    for (itf = _facetsOf[0].begin(); itf != _facetsOf[0].end(); ++itf) {
154
        if (_operationType == Difference) {  // toggle normal
155
            std::swap(itf->_aclPoints[0], itf->_aclPoints[1]);
156
            itf->CalcNormal();
157
        }
158

159
        facets.push_back(*itf);
160
    }
161

162
    for (itf = _facetsOf[1].begin(); itf != _facetsOf[1].end(); ++itf) {
163
        facets.push_back(*itf);
164
    }
165

166
    _resultMesh = facets;
167

168
    // Base::Sequencer().stop();
169
    // _builder.saveToFile("c:/temp/vdbg.iv");
170

171
    MeshDefinitions::SetMinPointDistance(saveMinMeshDistance);
172
}
173

174
void SetOperations::Cut(std::set<FacetIndex>& facetsCuttingEdge0,
175
                        std::set<FacetIndex>& facetsCuttingEdge1)
176
{
177
    MeshFacetGrid grid1(_cutMesh0, 20);
178
    MeshFacetGrid grid2(_cutMesh1, 20);
179

180
    unsigned long ctGx1 {}, ctGy1 {}, ctGz1 {};
181
    grid1.GetCtGrids(ctGx1, ctGy1, ctGz1);
182

183
    for (auto gx1 = 0UL; gx1 < ctGx1; gx1++) {
184
        for (auto gy1 = 0UL; gy1 < ctGy1; gy1++) {
185
            for (auto gz1 = 0UL; gz1 < ctGz1; gz1++) {
186
                if (grid1.GetCtElements(gx1, gy1, gz1) > 0) {
187
                    std::vector<FacetIndex> vecFacets2;
188
                    grid2.Inside(grid1.GetBoundBox(gx1, gy1, gz1), vecFacets2);
189

190
                    if (!vecFacets2.empty()) {
191
                        std::set<FacetIndex> vecFacets1;
192
                        grid1.GetElements(gx1, gy1, gz1, vecFacets1);
193

194
                        std::set<FacetIndex>::iterator it1;
195
                        for (it1 = vecFacets1.begin(); it1 != vecFacets1.end(); ++it1) {
196
                            FacetIndex fidx1 = *it1;
197
                            MeshGeomFacet f1 = _cutMesh0.GetFacet(*it1);
198

199
                            std::vector<FacetIndex>::iterator it2;
200
                            for (it2 = vecFacets2.begin(); it2 != vecFacets2.end(); ++it2) {
201
                                FacetIndex fidx2 = *it2;
202
                                MeshGeomFacet f2 = _cutMesh1.GetFacet(fidx2);
203

204
                                MeshPoint p0, p1;
205

206
                                int isect = f1.IntersectWithFacet(f2, p0, p1);
207
                                if (isect > 0) {
208
                                    // optimize cut line if distance to nearest point is too small
209
                                    float minDist1 = _minDistanceToPoint,
210
                                          minDist2 = _minDistanceToPoint;
211
                                    MeshPoint np0 = p0, np1 = p1;
212
                                    for (int i = 0; i < 3; i++)  // NOLINT
213
                                    {
214
                                        float d1 = (f1._aclPoints[i] - p0).Length();
215
                                        float d2 = (f1._aclPoints[i] - p1).Length();
216
                                        if (d1 < minDist1) {
217
                                            minDist1 = d1;
218
                                            np0 = f1._aclPoints[i];
219
                                        }
220
                                        if (d2 < minDist2) {
221
                                            minDist2 = d2;
222
                                            p1 = f1._aclPoints[i];
223
                                        }
224
                                    }  // for (int i = 0; i < 3; i++)
225

226
                                    // optimize cut line if distance to nearest point is too small
227
                                    for (int i = 0; i < 3; i++)  // NOLINT
228
                                    {
229
                                        float d1 = (f2._aclPoints[i] - p0).Length();
230
                                        float d2 = (f2._aclPoints[i] - p1).Length();
231
                                        if (d1 < minDist1) {
232
                                            minDist1 = d1;
233
                                            np0 = f2._aclPoints[i];
234
                                        }
235
                                        if (d2 < minDist2) {
236
                                            minDist2 = d2;
237
                                            np1 = f2._aclPoints[i];
238
                                        }
239
                                    }  // for (int i = 0; i < 3; i++)
240

241
                                    MeshPoint mp0 = np0;
242
                                    MeshPoint mp1 = np1;
243

244
                                    if (mp0 != mp1) {
245
                                        facetsCuttingEdge0.insert(fidx1);
246
                                        facetsCuttingEdge1.insert(fidx2);
247

248
                                        _cutPoints.insert(mp0);
249
                                        _cutPoints.insert(mp1);
250

251
                                        std::pair<std::set<MeshPoint>::iterator, bool> pit0 =
252
                                            _cutPoints.insert(mp0);
253
                                        std::pair<std::set<MeshPoint>::iterator, bool> pit1 =
254
                                            _cutPoints.insert(mp1);
255

256
                                        _edges[Edge(mp0, mp1)] = EdgeInfo();
257

258
                                        _facet2points[0][fidx1].push_back(pit0.first);
259
                                        _facet2points[0][fidx1].push_back(pit1.first);
260
                                        _facet2points[1][fidx2].push_back(pit0.first);
261
                                        _facet2points[1][fidx2].push_back(pit1.first);
262
                                    }
263
                                    else {
264
                                        std::pair<std::set<MeshPoint>::iterator, bool> pit =
265
                                            _cutPoints.insert(mp0);
266

267
                                        // do not insert a facet when only one corner point cuts the
268
                                        // edge if (!((mp0 == f1._aclPoints[0]) || (mp0 ==
269
                                        // f1._aclPoints[1]) || (mp0 == f1._aclPoints[2])))
270
                                        {
271
                                            facetsCuttingEdge0.insert(fidx1);
272
                                            _facet2points[0][fidx1].push_back(pit.first);
273
                                        }
274

275
                                        // if (!((mp0 == f2._aclPoints[0]) || (mp0 ==
276
                                        // f2._aclPoints[1]) || (mp0 == f2._aclPoints[2])))
277
                                        {
278
                                            facetsCuttingEdge1.insert(fidx2);
279
                                            _facet2points[1][fidx2].push_back(pit.first);
280
                                        }
281
                                    }
282

283
                                }  // if (f1.IntersectWithFacet(f2, p0, p1))
284
                            }      // for (it2 = vecFacets2.begin(); it2 != vecFacets2.end(); ++it2)
285
                        }          // for (it1 = vecFacets1.begin(); it1 != vecFacets1.end(); ++it1)
286
                    }              // if (vecFacets2.size() > 0)
287
                }                  // if (grid1.GetCtElements(gx1, gy1, gz1) > 0)
288
            }                      // for (gz1 = 0; gz1 < ctGz1; gz1++)
289
        }                          // for (gy1 = 0; gy1 < ctGy1; gy1++)
290
    }                              // for (gx1 = 0; gx1 < ctGx1; gx1++)
291
}
292

293
void SetOperations::TriangulateMesh(const MeshKernel& cutMesh, int side)
294
{
295
    // Triangulate Mesh
296
    std::map<FacetIndex, std::list<std::set<MeshPoint>::iterator>>::iterator it1;
297
    for (it1 = _facet2points[side].begin(); it1 != _facet2points[side].end(); ++it1) {
298
        std::vector<Vector3f> points;
299
        std::set<MeshPoint> pointsSet;
300

301
        FacetIndex fidx = it1->first;
302
        MeshGeomFacet f = cutMesh.GetFacet(fidx);
303

304
        // if (side == 1)
305
        //     _builder.addSingleTriangle(f._aclPoints[0], f._aclPoints[1], f._aclPoints[2], 3, 0,
306
        //     1, 1);
307

308
        // facet corner points
309
        // const MeshFacet& mf = cutMesh._aclFacetArray[fidx];
310
        for (int i = 0; i < 3; i++)  // NOLINT
311
        {
312
            pointsSet.insert(f._aclPoints[i]);
313
            points.push_back(f._aclPoints[i]);
314
        }
315

316
        // triangulated facets
317
        std::list<std::set<MeshPoint>::iterator>::iterator it2;
318
        for (it2 = it1->second.begin(); it2 != it1->second.end(); ++it2) {
319
            if (pointsSet.find(*(*it2)) == pointsSet.end()) {
320
                pointsSet.insert(*(*it2));
321
                points.push_back(*(*it2));
322
            }
323
        }
324

325
        Vector3f normal = f.GetNormal();
326
        Vector3f base = points[0];
327
        Vector3f dirX = points[1] - points[0];
328
        dirX.Normalize();
329
        Vector3f dirY = dirX % normal;
330

331
        // project points to 2D plane
332
        std::vector<Vector3f>::iterator it;
333
        std::vector<Vector3f> vertices;
334
        for (it = points.begin(); it != points.end(); ++it) {
335
            Vector3f pv = *it;
336
            pv.TransformToCoordinateSystem(base, dirX, dirY);
337
            vertices.push_back(pv);
338
        }
339

340
        DelaunayTriangulator tria;
341
        tria.SetPolygon(vertices);
342
        tria.TriangulatePolygon();
343

344
        std::vector<MeshFacet> facets = tria.GetFacets();
345
        for (auto& it : facets) {
346
            if ((it._aulPoints[0] == it._aulPoints[1]) || (it._aulPoints[1] == it._aulPoints[2])
347
                || (it._aulPoints[2] == it._aulPoints[0])) {  // two same triangle corner points
348
                continue;
349
            }
350

351
            MeshGeomFacet facet(points[it._aulPoints[0]],
352
                                points[it._aulPoints[1]],
353
                                points[it._aulPoints[2]]);
354

355
            // if (side == 1)
356
            //  _builder.addSingleTriangle(facet._aclPoints[0], facet._aclPoints[1],
357
            //  facet._aclPoints[2], true, 3, 0, 1, 1);
358

359
            // if (facet.Area() < 0.0001f)
360
            //{ // too small facet
361
            //   continue;
362
            // }
363

364
            float dist0 =
365
                facet._aclPoints[0].DistanceToLine(facet._aclPoints[1],
366
                                                   facet._aclPoints[1] - facet._aclPoints[2]);
367
            float dist1 =
368
                facet._aclPoints[1].DistanceToLine(facet._aclPoints[0],
369
                                                   facet._aclPoints[0] - facet._aclPoints[2]);
370
            float dist2 =
371
                facet._aclPoints[2].DistanceToLine(facet._aclPoints[0],
372
                                                   facet._aclPoints[0] - facet._aclPoints[1]);
373

374
            if ((dist0 < _minDistanceToPoint) || (dist1 < _minDistanceToPoint)
375
                || (dist2 < _minDistanceToPoint)) {
376
                continue;
377
            }
378

379
            // dist0 = (facet._aclPoints[0] - facet._aclPoints[1]).Length();
380
            // dist1 = (facet._aclPoints[1] - facet._aclPoints[2]).Length();
381
            // dist2 = (facet._aclPoints[2] - facet._aclPoints[3]).Length();
382

383
            // if ((dist0 < _minDistanceToPoint) || (dist1 < _minDistanceToPoint) || (dist2 <
384
            // _minDistanceToPoint))
385
            //{
386
            //   continue;
387
            // }
388

389
            facet.CalcNormal();
390
            if ((facet.GetNormal() * f.GetNormal()) < 0.0f) {  // adjust normal
391
                std::swap(facet._aclPoints[0], facet._aclPoints[1]);
392
                facet.CalcNormal();
393
            }
394

395

396
            for (int j = 0; j < 3; j++) {
397
                std::map<Edge, EdgeInfo>::iterator eit =
398
                    _edges.find(Edge(facet._aclPoints[j], facet._aclPoints[(j + 1) % 3]));
399

400
                if (eit != _edges.end()) {
401

402
                    if (eit->second.fcounter[side] < 2) {
403
                        // if (side == 0)
404
                        //    _builder.addSingleTriangle(facet._aclPoints[0], facet._aclPoints[1],
405
                        //    facet._aclPoints[2], true, 3, 0, 1, 1);
406

407
                        eit->second.facet[side] = fidx;
408
                        eit->second.facets[side][eit->second.fcounter[side]] = facet;
409
                        eit->second.fcounter[side]++;
410
                        facet.SetFlag(
411
                            MeshFacet::MARKED);  // set all facets connected to an edge: MARKED
412
                    }
413
                }
414
            }
415

416
            _newMeshFacets[side].push_back(facet);
417

418
        }  // for (i = 0; i < (out->numberoftriangles * 3); i += 3)
419
    }      // for (it1 = _facet2points[side].begin(); it1 != _facet2points[side].end(); ++it1)
420
}
421

422
void SetOperations::CollectFacets(int side, float mult)
423
{
424
    // float distSave = MeshDefinitions::_fMinPointDistance;
425
    // MeshDefinitions::SetMinPointDistance(1.0e-4f);
426

427
    MeshKernel mesh;
428
    MeshBuilder mb(mesh);
429
    mb.Initialize(_newMeshFacets[side].size());
430
    std::vector<MeshGeomFacet>::iterator it;
431
    for (it = _newMeshFacets[side].begin(); it != _newMeshFacets[side].end(); ++it) {
432
        // if (it->IsFlag(MeshFacet::MARKED))
433
        //{
434
        //   _builder.addSingleTriangle(it->_aclPoints[0], it->_aclPoints[1], it->_aclPoints[2],
435
        //   true, 3.0, 0.0, 1.0, 1.0);
436
        // }
437
        mb.AddFacet(*it, true);
438
    }
439
    mb.Finish();
440

441
    MeshAlgorithm algo(mesh);
442
    algo.ResetFacetFlag(static_cast<MeshFacet::TFlagType>(MeshFacet::VISIT | MeshFacet::TMP0));
443

444
    // bool hasFacetsNotVisited = true; // until facets not visited
445
    // search for facet not visited
446
    MeshFacetArray::_TConstIterator itf;
447
    const MeshFacetArray& rFacets = mesh.GetFacets();
448
    for (itf = rFacets.begin(); itf != rFacets.end(); ++itf) {
449
        if (!itf->IsFlag(MeshFacet::VISIT)) {  // Facet found, visit neighbours
450
            std::vector<FacetIndex> facets;
451
            facets.push_back(itf - rFacets.begin());  // add seed facet
452
            CollectFacetVisitor visitor(mesh, facets, _edges, side, mult, _builder);
453
            mesh.VisitNeighbourFacets(visitor, itf - rFacets.begin());
454

455
            if (visitor._addFacets == 0) {  // mark all facets to add it to the result
456
                algo.SetFacetsFlag(facets, MeshFacet::TMP0);
457
            }
458
        }
459
    }
460

461
    // add all facets to the result vector
462
    for (itf = rFacets.begin(); itf != rFacets.end(); ++itf) {
463
        if (itf->IsFlag(MeshFacet::TMP0)) {
464
            _facetsOf[side].push_back(mesh.GetFacet(*itf));
465
        }
466
    }
467

468
    // MeshDefinitions::SetMinPointDistance(distSave);
469
}
470

471
SetOperations::CollectFacetVisitor::CollectFacetVisitor(const MeshKernel& mesh,
472
                                                        std::vector<FacetIndex>& facets,
473
                                                        std::map<Edge, EdgeInfo>& edges,
474
                                                        int side,
475
                                                        float mult,
476
                                                        Base::Builder3D& builder)
477
    : _facets(facets)
478
    , _mesh(mesh)
479
    , _edges(edges)
480
    , _side(side)
481
    , _mult(mult)
482
    , _builder(builder)
483
{}
484

485
bool SetOperations::CollectFacetVisitor::Visit(const MeshFacet& rclFacet,
486
                                               const MeshFacet& rclFrom,
487
                                               FacetIndex ulFInd,
488
                                               unsigned long ulLevel)
489
{
490
    (void)rclFacet;
491
    (void)rclFrom;
492
    (void)ulLevel;
493
    _facets.push_back(ulFInd);
494
    return true;
495
}
496

497
// static int matchCounter = 0;
498
bool SetOperations::CollectFacetVisitor::AllowVisit(const MeshFacet& rclFacet,
499
                                                    const MeshFacet& rclFrom,
500
                                                    FacetIndex ulFInd,
501
                                                    unsigned long ulLevel,
502
                                                    unsigned short neighbourIndex)
503
{
504
    (void)ulFInd;
505
    (void)ulLevel;
506
    if (rclFacet.IsFlag(MeshFacet::MARKED) && rclFrom.IsFlag(MeshFacet::MARKED)) {
507
        // facet connected to an edge
508
        PointIndex pt0 = rclFrom._aulPoints[neighbourIndex],
509
                   pt1 = rclFrom._aulPoints[(neighbourIndex + 1) % 3];
510
        Edge edge(_mesh.GetPoint(pt0), _mesh.GetPoint(pt1));
511

512
        std::map<Edge, EdgeInfo>::iterator it = _edges.find(edge);
513

514
        if (it != _edges.end()) {
515
            if (_addFacets == -1) {
516
                // determine if the facets should add or not only once
517
                MeshGeomFacet facet = _mesh.GetFacet(rclFrom);  // triangulated facet
518
                MeshGeomFacet facetOther =
519
                    it->second
520
                        .facets[1 - _side][0];  // triangulated facet from same edge and other mesh
521
                Vector3f normalOther = facetOther.GetNormal();
522
                // Vector3f normal = facet.GetNormal();
523

524
                Vector3f edgeDir = it->first.pt1 - it->first.pt2;
525
                Vector3f ocDir = (edgeDir % (facet.GetGravityPoint() - it->first.pt1)) % edgeDir;
526
                ocDir.Normalize();
527
                Vector3f ocDirOther =
528
                    (edgeDir % (facetOther.GetGravityPoint() - it->first.pt1)) % edgeDir;
529
                ocDirOther.Normalize();
530

531
                // Vector3f dir = ocDir % normal;
532
                // Vector3f dirOther = ocDirOther % normalOther;
533

534
                bool match = ((ocDir * normalOther) * _mult) < 0.0f;
535

536
                // if (matchCounter == 1)
537
                //{
538
                //   // _builder.addSingleArrow(it->second.pt1, it->second.pt1 + edgeDir, 3,
539
                //   0.0, 1.0, 0.0);
540

541
                //  _builder.addSingleTriangle(facet._aclPoints[0], facet._aclPoints[1],
542
                //  facet._aclPoints[2], true, 3.0, 1.0, 0.0, 0.0);
543
                //  // _builder.addSingleArrow(facet.GetGravityPoint(), facet.GetGravityPoint() +
544
                //  ocDir, 3, 1.0, 0.0, 0.0); _builder.addSingleArrow(facet.GetGravityPoint(),
545
                //  facet.GetGravityPoint() + normal, 3, 1.0, 0.5, 0.0);
546
                //  // _builder.addSingleArrow(facet.GetGravityPoint(), facet.GetGravityPoint() +
547
                //  dir, 3, 1.0, 1.0, 0.0);
548

549
                //  _builder.addSingleTriangle(facetOther._aclPoints[0], facetOther._aclPoints[1],
550
                //  facetOther._aclPoints[2], true, 3.0, 0.0, 0.0, 1.0);
551
                //  // _builder.addSingleArrow(facetOther.GetGravityPoint(),
552
                //  facetOther.GetGravityPoint() + ocDirOther, 3, 0.0, 0.0, 1.0);
553
                //  _builder.addSingleArrow(facetOther.GetGravityPoint(),
554
                //  facetOther.GetGravityPoint() + normalOther, 3, 0.0, 0.5, 1.0);
555
                //  // _builder.addSingleArrow(facetOther.GetGravityPoint(),
556
                //  facetOther.GetGravityPoint() + dirOther, 3, 0.0, 1.0, 1.0);
557

558
                //}
559

560
                // float scalar = dir * dirOther * _mult;
561
                // bool match = scalar > 0.0f;
562

563

564
                // MeshPoint pt0 = it->first.pt1;
565
                // MeshPoint pt1 = it->first.pt2;
566

567
                // int i, n0 = -1, n1 = -1, m0 = -1, m1 = -1;
568
                // for (i = 0; i < 3; i++)
569
                //{
570
                //   if ((n0 == -1) && (facet._aclPoints[i] == pt0))
571
                //     n0 = i;
572
                //   if ((n1 == -1) && (facet._aclPoints[i] == pt1))
573
                //     n1 = i;
574
                //   if ((m0 == -1) && (facetOther._aclPoints[i] == pt0))
575
                //     m0 = i;
576
                //   if ((m1 == -1) && (facetOther._aclPoints[i] == pt1))
577
                //     m1 = i;
578
                // }
579

580
                // if ((n0 != -1) && (n1 != -1) && (m0 != -1) && (m1 != -1))
581
                //{
582
                //   bool orient_n = n1 > n0;
583
                //   bool orient_m = m1 > m0;
584

585
                //  Vector3f dirN = facet._aclPoints[n1] - facet._aclPoints[n0];
586
                //  Vector3f dirM = facetOther._aclPoints[m1] - facetOther._aclPoints[m0];
587

588
                //  if (matchCounter == 1)
589
                //  {
590
                //    _builder.addSingleArrow(facet.GetGravityPoint(), facet.GetGravityPoint() +
591
                //    dirN, 3, 1.0, 1.0, 0.0); _builder.addSingleArrow(facetOther.GetGravityPoint(),
592
                //    facetOther.GetGravityPoint() + dirM, 3, 0.0, 1.0, 1.0);
593
                //  }
594

595
                //  if (_mult > 0.0)
596
                //    match = orient_n == orient_m;
597
                //  else
598
                //    match = orient_n != orient_m;
599
                //}
600

601
                if (match) {
602
                    _addFacets = 0;
603
                }
604
                else {
605
                    _addFacets = 1;
606
                }
607

608
                // matchCounter++;
609
            }
610

611
            return false;
612
        }
613
    }
614

615
    return true;
616
}
617

618
// ----------------------------------------------------------------------------
619

620
bool MeshIntersection::hasIntersection() const
621
{
622
    Base::BoundBox3f bbox1 = kernel1.GetBoundBox();
623
    Base::BoundBox3f bbox2 = kernel2.GetBoundBox();
624
    if (!(bbox1 && bbox2)) {
625
        return false;
626
    }
627

628
    if (testIntersection(kernel1, kernel2)) {
629
        return true;
630
    }
631

632
    return false;
633
}
634

635
void MeshIntersection::getIntersection(std::list<MeshIntersection::Tuple>& intsct) const
636
{
637
    const MeshKernel& k1 = kernel1;
638
    const MeshKernel& k2 = kernel2;
639

640
    // Contains bounding boxes for every facet of 'k1'
641
    std::vector<Base::BoundBox3f> boxes1;
642
    MeshFacetIterator cMFI1(k1);
643
    for (cMFI1.Begin(); cMFI1.More(); cMFI1.Next()) {
644
        boxes1.push_back((*cMFI1).GetBoundBox());
645
    }
646

647
    // Contains bounding boxes for every facet of 'k2'
648
    std::vector<Base::BoundBox3f> boxes2;
649
    MeshFacetIterator cMFI2(k2);
650
    for (cMFI2.Begin(); cMFI2.More(); cMFI2.Next()) {
651
        boxes2.push_back((*cMFI2).GetBoundBox());
652
    }
653

654
    // Splits the mesh using grid for speeding up the calculation
655
    MeshFacetGrid cMeshFacetGrid(k1);
656

657
    const MeshFacetArray& rFaces2 = k2.GetFacets();
658
    Base::SequencerLauncher seq("Checking for intersections...", rFaces2.size());
659
    int index = 0;
660
    MeshGeomFacet facet1, facet2;
661
    Base::Vector3f pt1, pt2;
662

663
    // Iterate over the facets of the 2nd mesh and find the grid elements of the 1st mesh
664
    for (MeshFacetArray::_TConstIterator it = rFaces2.begin(); it != rFaces2.end(); ++it, index++) {
665
        seq.next();
666
        std::vector<FacetIndex> elements;
667
        cMeshFacetGrid.Inside(boxes2[index], elements, true);
668

669
        cMFI2.Set(index);
670
        facet2 = *cMFI2;
671

672
        for (FacetIndex element : elements) {
673
            if (boxes2[index] && boxes1[element]) {
674
                cMFI1.Set(element);
675
                facet1 = *cMFI1;
676
                int ret = facet1.IntersectWithFacet(facet2, pt1, pt2);
677
                if (ret == 2) {
678
                    Tuple d;
679
                    d.p1 = pt1;
680
                    d.p2 = pt2;
681
                    d.f1 = element;
682
                    d.f2 = index;
683
                    intsct.push_back(d);
684
                }
685
            }
686
        }
687
    }
688
}
689

690
bool MeshIntersection::testIntersection(const MeshKernel& k1, const MeshKernel& k2)
691
{
692
    // Contains bounding boxes for every facet of 'k1'
693
    std::vector<Base::BoundBox3f> boxes1;
694
    MeshFacetIterator cMFI1(k1);
695
    for (cMFI1.Begin(); cMFI1.More(); cMFI1.Next()) {
696
        boxes1.push_back((*cMFI1).GetBoundBox());
697
    }
698

699
    // Contains bounding boxes for every facet of 'k2'
700
    std::vector<Base::BoundBox3f> boxes2;
701
    MeshFacetIterator cMFI2(k2);
702
    for (cMFI2.Begin(); cMFI2.More(); cMFI2.Next()) {
703
        boxes2.push_back((*cMFI2).GetBoundBox());
704
    }
705

706
    // Splits the mesh using grid for speeding up the calculation
707
    MeshFacetGrid cMeshFacetGrid(k1);
708

709
    const MeshFacetArray& rFaces2 = k2.GetFacets();
710
    Base::SequencerLauncher seq("Checking for intersections...", rFaces2.size());
711
    int index = 0;
712
    MeshGeomFacet facet1, facet2;
713
    Base::Vector3f pt1, pt2;
714

715
    // Iterate over the facets of the 2nd mesh and find the grid elements of the 1st mesh
716
    for (MeshFacetArray::_TConstIterator it = rFaces2.begin(); it != rFaces2.end(); ++it, index++) {
717
        seq.next();
718
        std::vector<FacetIndex> elements;
719
        cMeshFacetGrid.Inside(boxes2[index], elements, true);
720

721
        cMFI2.Set(index);
722
        facet2 = *cMFI2;
723

724
        for (FacetIndex element : elements) {
725
            if (boxes2[index] && boxes1[element]) {
726
                cMFI1.Set(element);
727
                facet1 = *cMFI1;
728
                int ret = facet1.IntersectWithFacet(facet2, pt1, pt2);
729
                if (ret == 2) {
730
                    // abort after the first detected self-intersection
731
                    return true;
732
                }
733
            }
734
        }
735
    }
736

737
    return false;
738
}
739

740
void MeshIntersection::connectLines(bool onlyclosed,
741
                                    const std::list<MeshIntersection::Tuple>& rdata,
742
                                    std::list<std::list<MeshIntersection::Triple>>& lines)
743
{
744
    float fMinEps = minDistance * minDistance;
745

746
    std::list<Tuple> data = rdata;
747
    while (!data.empty()) {
748
        std::list<Tuple>::iterator pF;
749
        std::list<Triple> newPoly;
750

751
        // add first line and delete from the list
752
        Triple front, back;
753
        front.f1 = data.begin()->f1;
754
        front.f2 = data.begin()->f2;
755
        front.p = data.begin()->p1;  // current start point of the polyline
756
        back.f1 = data.begin()->f1;
757
        back.f2 = data.begin()->f2;
758
        back.p = data.begin()->p2;  // current end point of the polyline
759
        newPoly.push_back(front);
760
        newPoly.push_back(back);
761
        data.erase(data.begin());
762

763
        // search for the next line on the begin/end of the polyline and add it
764
        std::list<Tuple>::iterator pFront, pEnd;
765
        bool bFoundLine {};
766
        do {
767
            float fFrontMin = fMinEps, fEndMin = fMinEps;
768
            bool bFrontFirst = false, bEndFirst = false;
769

770
            pFront = data.end();
771
            pEnd = data.end();
772
            bFoundLine = false;
773

774
            for (pF = data.begin(); pF != data.end(); ++pF) {
775
                if (Base::DistanceP2(front.p, pF->p1) < fFrontMin) {
776
                    fFrontMin = Base::DistanceP2(front.p, pF->p1);
777
                    pFront = pF;
778
                    bFrontFirst = true;
779
                }
780
                else if (Base::DistanceP2(back.p, pF->p1) < fEndMin) {
781
                    fEndMin = Base::DistanceP2(back.p, pF->p1);
782
                    pEnd = pF;
783
                    bEndFirst = true;
784
                }
785
                else if (Base::DistanceP2(front.p, pF->p2) < fFrontMin) {
786
                    fFrontMin = Base::DistanceP2(front.p, pF->p2);
787
                    pFront = pF;
788
                    bFrontFirst = false;
789
                }
790
                else if (Base::DistanceP2(back.p, pF->p2) < fEndMin) {
791
                    fEndMin = Base::DistanceP2(back.p, pF->p2);
792
                    pEnd = pF;
793
                    bEndFirst = false;
794
                }
795

796
                if (fFrontMin == 0.0f || fEndMin == 0.0f) {
797
                    break;
798
                }
799
            }
800

801
            if (pFront != data.end()) {
802
                bFoundLine = true;
803
                if (bFrontFirst) {
804
                    front.f1 = pFront->f1;
805
                    front.f2 = pFront->f2;
806
                    front.p = pFront->p2;
807
                    newPoly.push_front(front);
808
                }
809
                else {
810
                    front.f1 = pFront->f1;
811
                    front.f2 = pFront->f2;
812
                    front.p = pFront->p1;
813
                    newPoly.push_front(front);
814
                }
815

816
                data.erase(pFront);
817
            }
818

819
            if (pEnd != data.end()) {
820
                bFoundLine = true;
821
                if (bEndFirst) {
822
                    back.f1 = pEnd->f1;
823
                    back.f2 = pEnd->f2;
824
                    back.p = pEnd->p2;
825
                    newPoly.push_back(back);
826
                }
827
                else {
828
                    back.f1 = pEnd->f1;
829
                    back.f2 = pEnd->f2;
830
                    back.p = pEnd->p1;
831
                    newPoly.push_back(back);
832
                }
833

834
                data.erase(pEnd);
835
            }
836
        } while (bFoundLine);
837

838
        if (onlyclosed) {
839
            if (newPoly.size() > 2
840
                && Base::DistanceP2(newPoly.front().p, newPoly.back().p) < fMinEps) {
841
                lines.push_back(newPoly);
842
            }
843
        }
844
        else {
845
            lines.push_back(newPoly);
846
        }
847
    }
848
}
849