FreeCAD

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

25
#ifndef _PreComp_
26
#include <algorithm>
27
#include <map>
28
#include <queue>
29
#include <stdexcept>
30
#endif
31

32
#include <Base/Exception.h>
33
#include <Base/Stream.h>
34
#include <Base/Swap.h>
35

36
#include "Algorithm.h"
37
#include "Builder.h"
38
#include "Evaluation.h"
39
#include "Iterator.h"
40
#include "MeshIO.h"
41
#include "MeshKernel.h"
42
#include "Smoothing.h"
43

44

45
using namespace MeshCore;
46

47
MeshKernel::MeshKernel()
48
{
49
    _clBoundBox.SetVoid();
50
}
51

52
MeshKernel::MeshKernel(const MeshKernel& rclMesh)
53
{
54
    *this = rclMesh;
55
}
56

57
MeshKernel::MeshKernel(MeshKernel&& rclMesh)
58
{
59
    *this = rclMesh;
60
}
61

62
MeshKernel& MeshKernel::operator=(const MeshKernel& rclMesh)
63
{
64
    if (this != &rclMesh) {  // must be a different instance
65
        this->_aclPointArray = rclMesh._aclPointArray;
66
        this->_aclFacetArray = rclMesh._aclFacetArray;
67
        this->_clBoundBox = rclMesh._clBoundBox;
68
        this->_bValid = rclMesh._bValid;
69
    }
70
    return *this;
71
}
72

73
MeshKernel& MeshKernel::operator=(MeshKernel&& rclMesh)
74
{
75
    if (this != &rclMesh) {  // must be a different instance
76
        this->_aclPointArray = std::move(rclMesh._aclPointArray);
77
        this->_aclFacetArray = std::move(rclMesh._aclFacetArray);
78
        this->_clBoundBox = rclMesh._clBoundBox;
79
        this->_bValid = rclMesh._bValid;
80
    }
81
    return *this;
82
}
83

84
MeshKernel& MeshKernel::operator=(const std::vector<MeshGeomFacet>& rclFAry)
85
{
86
    MeshBuilder builder(*this);
87
    builder.Initialize(rclFAry.size());
88

89
    for (const auto& it : rclFAry) {
90
        builder.AddFacet(it);
91
    }
92

93
    builder.Finish();
94

95
    return *this;
96
}
97

98
void MeshKernel::Assign(const MeshPointArray& rPoints,
99
                        const MeshFacetArray& rFacets,
100
                        bool checkNeighbourHood)
101
{
102
    _aclPointArray = rPoints;
103
    _aclFacetArray = rFacets;
104
    RecalcBoundBox();
105
    if (checkNeighbourHood) {
106
        RebuildNeighbours();
107
    }
108
}
109

110
void MeshKernel::Adopt(MeshPointArray& rPoints, MeshFacetArray& rFacets, bool checkNeighbourHood)
111
{
112
    _aclPointArray.swap(rPoints);
113
    _aclFacetArray.swap(rFacets);
114
    RecalcBoundBox();
115
    if (checkNeighbourHood) {
116
        RebuildNeighbours();
117
    }
118
}
119

120
void MeshKernel::Swap(MeshKernel& mesh)
121
{
122
    this->_aclPointArray.swap(mesh._aclPointArray);
123
    this->_aclFacetArray.swap(mesh._aclFacetArray);
124
    this->_clBoundBox = mesh._clBoundBox;
125
}
126

127
MeshKernel& MeshKernel::operator+=(const MeshGeomFacet& rclSFacet)
128
{
129
    this->AddFacet(rclSFacet);
130
    return *this;
131
}
132

133
void MeshKernel::AddFacet(const MeshGeomFacet& rclSFacet)
134
{
135
    MeshFacet clFacet;
136

137
    // set corner points
138
    for (int i = 0; i < 3; i++) {
139
        _clBoundBox.Add(rclSFacet._aclPoints[i]);
140
        clFacet._aulPoints[i] = _aclPointArray.GetOrAddIndex(rclSFacet._aclPoints[i]);
141
    }
142

143
    // adjust orientation to normal
144
    AdjustNormal(clFacet, rclSFacet.GetNormal());
145

146
    FacetIndex ulCt = _aclFacetArray.size();
147

148
    // set neighbourhood
149
    PointIndex ulP0 = clFacet._aulPoints[0];
150
    PointIndex ulP1 = clFacet._aulPoints[1];
151
    PointIndex ulP2 = clFacet._aulPoints[2];
152
    FacetIndex ulCC = 0;
153
    for (TMeshFacetArray::iterator pF = _aclFacetArray.begin(); pF != _aclFacetArray.end();
154
         ++pF, ulCC++) {
155
        for (int i = 0; i < 3; i++) {
156
            PointIndex ulP = pF->_aulPoints[i];
157
            PointIndex ulQ = pF->_aulPoints[(i + 1) % 3];
158
            if (ulQ == ulP0 && ulP == ulP1) {
159
                clFacet._aulNeighbours[0] = ulCC;
160
                pF->_aulNeighbours[i] = ulCt;
161
            }
162
            else if (ulQ == ulP1 && ulP == ulP2) {
163
                clFacet._aulNeighbours[1] = ulCC;
164
                pF->_aulNeighbours[i] = ulCt;
165
            }
166
            else if (ulQ == ulP2 && ulP == ulP0) {
167
                clFacet._aulNeighbours[2] = ulCC;
168
                pF->_aulNeighbours[i] = ulCt;
169
            }
170
        }
171
    }
172

173
    // insert facet into array
174
    _aclFacetArray.push_back(clFacet);
175
}
176

177
MeshKernel& MeshKernel::operator+=(const std::vector<MeshGeomFacet>& rclFAry)
178
{
179
    this->AddFacets(rclFAry);
180
    return *this;
181
}
182

183
void MeshKernel::AddFacets(const std::vector<MeshGeomFacet>& rclFAry)
184
{
185
    // Create a temp. kernel to get the topology of the passed triangles
186
    // and merge them with this kernel. This keeps properties and flags
187
    // of this mesh.
188
    MeshKernel tmp;
189
    tmp = rclFAry;
190
    Merge(tmp);
191
}
192

193
unsigned long MeshKernel::AddFacets(const std::vector<MeshFacet>& rclFAry, bool checkManifolds)
194
{
195
    // Build map of edges of the referencing facets we want to append
196
#ifdef FC_DEBUG
197
    unsigned long countPoints = CountPoints();
198
#endif
199

200
    // if the manifold check shouldn't be done then just add all faces
201
    if (!checkManifolds) {
202
        FacetIndex countFacets = CountFacets();
203
        FacetIndex countValid = rclFAry.size();
204
        _aclFacetArray.reserve(countFacets + countValid);
205

206
        // just add all faces now
207
        for (const auto& pF : rclFAry) {
208
            _aclFacetArray.push_back(pF);
209
        }
210

211
        RebuildNeighbours(countFacets);
212
        return _aclFacetArray.size();
213
    }
214

215
    this->_aclPointArray.ResetInvalid();
216
    FacetIndex k = CountFacets();
217
    std::map<std::pair<PointIndex, PointIndex>, std::list<FacetIndex>> edgeMap;
218
    for (std::vector<MeshFacet>::const_iterator pF = rclFAry.begin(); pF != rclFAry.end();
219
         ++pF, k++) {
220
        // reset INVALID flag for all candidates
221
        pF->ResetFlag(MeshFacet::INVALID);
222
        for (int i = 0; i < 3; i++) {
223
#ifdef FC_DEBUG
224
            assert(pF->_aulPoints[i] < countPoints);
225
#endif
226
            this->_aclPointArray[pF->_aulPoints[i]].SetFlag(MeshPoint::INVALID);
227
            PointIndex ulT0 = pF->_aulPoints[i];
228
            PointIndex ulT1 = pF->_aulPoints[(i + 1) % 3];
229
            PointIndex ulP0 = std::min<PointIndex>(ulT0, ulT1);
230
            PointIndex ulP1 = std::max<PointIndex>(ulT0, ulT1);
231
            edgeMap[std::make_pair(ulP0, ulP1)].push_front(k);
232
        }
233
    }
234

235
    // Check for the above edges in the current facet array
236
    k = 0;
237
    for (MeshFacetArray::_TIterator pF = _aclFacetArray.begin(); pF != _aclFacetArray.end();
238
         ++pF, k++) {
239
        // if none of the points references one of the edges ignore the facet
240
        if (!this->_aclPointArray[pF->_aulPoints[0]].IsFlag(MeshPoint::INVALID)
241
            && !this->_aclPointArray[pF->_aulPoints[1]].IsFlag(MeshPoint::INVALID)
242
            && !this->_aclPointArray[pF->_aulPoints[2]].IsFlag(MeshPoint::INVALID)) {
243
            continue;
244
        }
245
        for (int i = 0; i < 3; i++) {
246
            PointIndex ulT0 = pF->_aulPoints[i];
247
            PointIndex ulT1 = pF->_aulPoints[(i + 1) % 3];
248
            PointIndex ulP0 = std::min<PointIndex>(ulT0, ulT1);
249
            PointIndex ulP1 = std::max<PointIndex>(ulT0, ulT1);
250
            std::pair<PointIndex, PointIndex> edge = std::make_pair(ulP0, ulP1);
251
            std::map<std::pair<PointIndex, PointIndex>, std::list<FacetIndex>>::iterator pI =
252
                edgeMap.find(edge);
253
            // Does the current facet share the same edge?
254
            if (pI != edgeMap.end()) {
255
                pI->second.push_front(k);
256
            }
257
        }
258
    }
259

260
    this->_aclPointArray.ResetInvalid();
261

262
    // Now let's see for which edges we might get manifolds, if so we don't add the corresponding
263
    // candidates
264
    FacetIndex countFacets = CountFacets();
265
    std::map<std::pair<PointIndex, PointIndex>, std::list<FacetIndex>>::iterator pE;
266
    for (pE = edgeMap.begin(); pE != edgeMap.end(); ++pE) {
267
        if (pE->second.size() > 2) {
268
            for (FacetIndex it : pE->second) {
269
                if (it >= countFacets) {
270
                    // this is a candidate
271
                    FacetIndex index = it - countFacets;
272
                    rclFAry[index].SetFlag(MeshFacet::INVALID);
273
                }
274
            }
275
        }
276
    }
277

278
    // Do not insert directly to the data structure because we should get the correct size of new
279
    // facets, otherwise std::vector reallocates too much memory which can't be freed so easily
280
    MeshIsNotFlag<MeshFacet> flag;
281
    FacetIndex countValid =
282
        std::count_if(rclFAry.begin(), rclFAry.end(), [flag](const MeshFacet& f) {
283
            return flag(f, MeshFacet::INVALID);
284
        });
285
    _aclFacetArray.reserve(_aclFacetArray.size() + countValid);
286
    // now start inserting the facets to the data structure and set the correct neighbourhood as
287
    // well
288
    FacetIndex startIndex = CountFacets();
289
    for (const auto& pF : rclFAry) {
290
        if (!pF.IsFlag(MeshFacet::INVALID)) {
291
            _aclFacetArray.push_back(pF);
292
            pF.SetProperty(startIndex++);
293
        }
294
    }
295

296
    // resolve neighbours
297
    for (pE = edgeMap.begin(); pE != edgeMap.end(); ++pE) {
298
        PointIndex ulP0 = pE->first.first;
299
        PointIndex ulP1 = pE->first.second;
300
        if (pE->second.size() == 1)  // border facet
301
        {
302
            FacetIndex ulF0 = pE->second.front();
303
            if (ulF0 >= countFacets) {
304
                ulF0 -= countFacets;
305
                std::vector<MeshFacet>::const_iterator pF = rclFAry.begin() + ulF0;
306
                if (!pF->IsFlag(MeshFacet::INVALID)) {
307
                    ulF0 = pF->_ulProp;
308
                }
309
                else {
310
                    ulF0 = FACET_INDEX_MAX;
311
                }
312
            }
313

314
            if (ulF0 != FACET_INDEX_MAX) {
315
                unsigned short usSide = _aclFacetArray[ulF0].Side(ulP0, ulP1);
316
                assert(usSide != USHRT_MAX);
317
                _aclFacetArray[ulF0]._aulNeighbours[usSide] = FACET_INDEX_MAX;
318
            }
319
        }
320
        else if (pE->second.size() == 2)  // normal facet with neighbour
321
        {
322
            // we must check if both facets are part of the mesh now
323
            FacetIndex ulF0 = pE->second.front();
324
            if (ulF0 >= countFacets) {
325
                ulF0 -= countFacets;
326
                std::vector<MeshFacet>::const_iterator pF = rclFAry.begin() + ulF0;
327
                if (!pF->IsFlag(MeshFacet::INVALID)) {
328
                    ulF0 = pF->_ulProp;
329
                }
330
                else {
331
                    ulF0 = FACET_INDEX_MAX;
332
                }
333
            }
334
            FacetIndex ulF1 = pE->second.back();
335
            if (ulF1 >= countFacets) {
336
                ulF1 -= countFacets;
337
                std::vector<MeshFacet>::const_iterator pF = rclFAry.begin() + ulF1;
338
                if (!pF->IsFlag(MeshFacet::INVALID)) {
339
                    ulF1 = pF->_ulProp;
340
                }
341
                else {
342
                    ulF1 = FACET_INDEX_MAX;
343
                }
344
            }
345

346
            if (ulF0 != FACET_INDEX_MAX) {
347
                unsigned short usSide = _aclFacetArray[ulF0].Side(ulP0, ulP1);
348
                assert(usSide != USHRT_MAX);
349
                _aclFacetArray[ulF0]._aulNeighbours[usSide] = ulF1;
350
            }
351

352
            if (ulF1 != FACET_INDEX_MAX) {
353
                unsigned short usSide = _aclFacetArray[ulF1].Side(ulP0, ulP1);
354
                assert(usSide != USHRT_MAX);
355
                _aclFacetArray[ulF1]._aulNeighbours[usSide] = ulF0;
356
            }
357
        }
358
    }
359

360
    return _aclFacetArray.size();
361
}
362

363
unsigned long MeshKernel::AddFacets(const std::vector<MeshFacet>& rclFAry,
364
                                    const std::vector<Base::Vector3f>& rclPAry,
365
                                    bool checkManifolds)
366
{
367
    for (auto it : rclPAry) {
368
        _clBoundBox.Add(it);
369
    }
370
    this->_aclPointArray.insert(this->_aclPointArray.end(), rclPAry.begin(), rclPAry.end());
371
    return this->AddFacets(rclFAry, checkManifolds);
372
}
373

374
void MeshKernel::Merge(const MeshKernel& rKernel)
375
{
376
    if (this != &rKernel) {
377
        const MeshPointArray& rPoints = rKernel._aclPointArray;
378
        const MeshFacetArray& rFacets = rKernel._aclFacetArray;
379
        Merge(rPoints, rFacets);
380
    }
381
}
382

383
void MeshKernel::Merge(const MeshPointArray& rPoints, const MeshFacetArray& rFaces)
384
{
385
    if (rPoints.empty() || rFaces.empty()) {
386
        return;  // nothing to do
387
    }
388
    std::vector<PointIndex> increments(rPoints.size());
389

390
    FacetIndex countFacets = this->_aclFacetArray.size();
391
    // Reserve the additional memory to append the new facets
392
    this->_aclFacetArray.reserve(this->_aclFacetArray.size() + rFaces.size());
393

394
    // Copy the new faces immediately to the facet array
395
    MeshFacet face;
396
    for (const auto& it : rFaces) {
397
        face = it;
398
        for (PointIndex point : it._aulPoints) {
399
            increments[point]++;
400
        }
401

402
        // append to the facet array
403
        this->_aclFacetArray.push_back(face);
404
    }
405

406
    std::size_t countNewPoints =
407
        std::count_if(increments.begin(), increments.end(), [](PointIndex v) {
408
            return v > 0;
409
        });
410
    // Reserve the additional memory to append the new points
411
    PointIndex index = this->_aclPointArray.size();
412
    this->_aclPointArray.reserve(this->_aclPointArray.size() + countNewPoints);
413

414
    // Now we can start inserting the points and adjust the point indices of the faces
415
    for (std::vector<PointIndex>::iterator it = increments.begin(); it != increments.end(); ++it) {
416
        if (*it > 0) {
417
            // set the index of the point array
418
            *it = index++;
419
            const MeshPoint& rPt = rPoints[it - increments.begin()];
420
            this->_aclPointArray.push_back(rPt);
421
            _clBoundBox.Add(rPt);
422
        }
423
    }
424

425
    for (MeshFacetArray::_TIterator pF = this->_aclFacetArray.begin() + countFacets;
426
         pF != this->_aclFacetArray.end();
427
         ++pF) {
428
        for (PointIndex& index : pF->_aulPoints) {
429
            index = increments[index];
430
        }
431
    }
432

433
    // Since rFaces could consist of a subset of the actual facet array the
434
    // neighbour indices could be totally wrong so they must be rebuilt from
435
    // scratch. Fortunately, this needs only to be done for the newly inserted
436
    // facets -- not for all
437
    RebuildNeighbours(countFacets);
438
}
439

440
void MeshKernel::Cleanup()
441
{
442
    MeshCleanup meshCleanup(_aclPointArray, _aclFacetArray);
443
    meshCleanup.RemoveInvalids();
444
}
445

446
void MeshKernel::Clear()
447
{
448
    _aclPointArray.clear();
449
    _aclFacetArray.clear();
450

451
    // release memory
452
    MeshPointArray().swap(_aclPointArray);
453
    MeshFacetArray().swap(_aclFacetArray);
454

455
    _clBoundBox.SetVoid();
456
}
457

458
bool MeshKernel::DeleteFacet(const MeshFacetIterator& rclIter)
459
{
460
    FacetIndex ulNFacet {}, ulInd {};
461

462
    if (rclIter._clIter >= _aclFacetArray.end()) {
463
        return false;
464
    }
465

466
    // index of the facet to delete
467
    ulInd = rclIter._clIter - _aclFacetArray.begin();
468

469
    // invalidate neighbour indices of the neighbour facet to this facet
470
    for (FacetIndex nbIndex : rclIter._clIter->_aulNeighbours) {
471
        ulNFacet = nbIndex;
472
        if (ulNFacet != FACET_INDEX_MAX) {
473
            for (FacetIndex& nbOfNb : _aclFacetArray[ulNFacet]._aulNeighbours) {
474
                if (nbOfNb == ulInd) {
475
                    nbOfNb = FACET_INDEX_MAX;
476
                    break;
477
                }
478
            }
479
        }
480
    }
481

482
    // erase corner point if needed
483
    for (int i = 0; i < 3; i++) {
484
        if ((rclIter._clIter->_aulNeighbours[i] == FACET_INDEX_MAX)
485
            && (rclIter._clIter->_aulNeighbours[(i + 1) % 3] == FACET_INDEX_MAX)) {
486
            // no neighbours, possibly delete point
487
            ErasePoint(rclIter._clIter->_aulPoints[(i + 1) % 3], ulInd);
488
        }
489
    }
490

491
    // remove facet from array
492
    _aclFacetArray.Erase(_aclFacetArray.begin() + rclIter.Position());
493

494
    return true;
495
}
496

497
bool MeshKernel::DeleteFacet(FacetIndex ulInd)
498
{
499
    if (ulInd >= _aclFacetArray.size()) {
500
        return false;
501
    }
502

503
    MeshFacetIterator clIter(*this);
504
    clIter.Set(ulInd);
505

506
    return DeleteFacet(clIter);
507
}
508

509
void MeshKernel::DeleteFacets(const std::vector<FacetIndex>& raulFacets)
510
{
511
    _aclPointArray.SetProperty(0);
512

513
    // number of referencing facets per point
514
    for (const auto& pF : _aclFacetArray) {
515
        _aclPointArray[pF._aulPoints[0]]._ulProp++;
516
        _aclPointArray[pF._aulPoints[1]]._ulProp++;
517
        _aclPointArray[pF._aulPoints[2]]._ulProp++;
518
    }
519

520
    // invalidate facet and adjust number of point references
521
    _aclFacetArray.ResetInvalid();
522
    for (FacetIndex index : raulFacets) {
523
        MeshFacet& rclFacet = _aclFacetArray[index];
524
        rclFacet.SetInvalid();
525
        _aclPointArray[rclFacet._aulPoints[0]]._ulProp--;
526
        _aclPointArray[rclFacet._aulPoints[1]]._ulProp--;
527
        _aclPointArray[rclFacet._aulPoints[2]]._ulProp--;
528
    }
529

530
    // invalidate all unreferenced points
531
    _aclPointArray.ResetInvalid();
532
    for (auto& pP : _aclPointArray) {
533
        if (pP._ulProp == 0) {
534
            pP.SetInvalid();
535
        }
536
    }
537

538
    RemoveInvalids();
539
    RecalcBoundBox();
540
}
541

542
bool MeshKernel::DeletePoint(PointIndex ulInd)
543
{
544
    if (ulInd >= _aclPointArray.size()) {
545
        return false;
546
    }
547

548
    MeshPointIterator clIter(*this);
549
    clIter.Set(ulInd);
550

551
    return DeletePoint(clIter);
552
}
553

554
bool MeshKernel::DeletePoint(const MeshPointIterator& rclIter)
555
{
556
    MeshFacetIterator pFIter(*this), pFEnd(*this);
557
    std::vector<MeshFacetIterator> clToDel;
558
    PointIndex ulInd {};
559

560
    // index of the point to delete
561
    ulInd = rclIter._clIter - _aclPointArray.begin();
562

563
    pFIter.Begin();
564
    pFEnd.End();
565

566
    // check corner points of all facets
567
    while (pFIter < pFEnd) {
568
        for (PointIndex ptIndex : pFIter._clIter->_aulPoints) {
569
            if (ulInd == ptIndex) {
570
                clToDel.push_back(pFIter);
571
            }
572
        }
573
        ++pFIter;
574
    }
575

576
    // iterators (facets) sort by index
577
    std::sort(clToDel.begin(), clToDel.end());
578

579
    // delete each facet separately (from back to front to avoid to
580
    // invalidate the iterators)
581
    for (size_t i = clToDel.size(); i > 0; i--) {
582
        DeleteFacet(clToDel[i - 1]);
583
    }
584
    return true;
585
}
586

587
void MeshKernel::DeletePoints(const std::vector<PointIndex>& raulPoints)
588
{
589
    _aclPointArray.ResetInvalid();
590
    for (PointIndex ptIndex : raulPoints) {
591
        _aclPointArray[ptIndex].SetInvalid();
592
    }
593

594
    // delete facets if at least one corner point is invalid
595
    _aclPointArray.SetProperty(0);
596
    for (auto& pF : _aclFacetArray) {
597
        MeshPoint& rclP0 = _aclPointArray[pF._aulPoints[0]];
598
        MeshPoint& rclP1 = _aclPointArray[pF._aulPoints[1]];
599
        MeshPoint& rclP2 = _aclPointArray[pF._aulPoints[2]];
600

601
        if (!rclP0.IsValid() || !rclP1.IsValid() || !rclP2.IsValid()) {
602
            pF.SetInvalid();
603
        }
604
        else {
605
            pF.ResetInvalid();
606
            rclP0._ulProp++;
607
            rclP1._ulProp++;
608
            rclP2._ulProp++;
609
        }
610
    }
611

612
    // invalidate all unreferenced points to delete them
613
    for (auto& pP : _aclPointArray) {
614
        if (pP._ulProp == 0) {
615
            pP.SetInvalid();
616
        }
617
    }
618

619
    RemoveInvalids();
620
    RecalcBoundBox();
621
}
622

623
void MeshKernel::ErasePoint(PointIndex ulIndex, FacetIndex ulFacetIndex, bool bOnlySetInvalid)
624
{
625
    std::vector<MeshFacet>::iterator pFIter, pFEnd, pFNot;
626

627
    pFIter = _aclFacetArray.begin();
628
    pFNot = _aclFacetArray.begin() + ulFacetIndex;
629
    pFEnd = _aclFacetArray.end();
630

631
    // check all facets
632
    while (pFIter < pFNot) {
633
        for (PointIndex ptIndex : pFIter->_aulPoints) {
634
            if (ptIndex == ulIndex) {
635
                return;  // point still referenced ==> do not delete
636
            }
637
        }
638
        ++pFIter;
639
    }
640

641
    ++pFIter;
642
    while (pFIter < pFEnd) {
643
        for (PointIndex ptIndex : pFIter->_aulPoints) {
644
            if (ptIndex == ulIndex) {
645
                return;  // point still referenced ==> do not delete
646
            }
647
        }
648
        ++pFIter;
649
    }
650

651

652
    if (!bOnlySetInvalid) {
653
        // completely remove point
654
        _aclPointArray.erase(_aclPointArray.begin() + ulIndex);
655

656
        // correct point indices of the facets
657
        pFIter = _aclFacetArray.begin();
658
        while (pFIter < pFEnd) {
659
            for (PointIndex& ptIndex : pFIter->_aulPoints) {
660
                if (ptIndex > ulIndex) {
661
                    ptIndex--;
662
                }
663
            }
664
            ++pFIter;
665
        }
666
    }
667
    else {  // only invalidate
668
        _aclPointArray[ulIndex].SetInvalid();
669
    }
670
}
671

672
void MeshKernel::RemoveInvalids()
673
{
674
    std::vector<unsigned long> aulDecrements;
675
    std::vector<unsigned long>::iterator pDIter;
676
    unsigned long ulDec {};
677
    MeshPointArray::_TIterator pPIter, pPEnd;
678
    MeshFacetArray::_TIterator pFIter, pFEnd;
679

680
    // generate array of decrements
681
    aulDecrements.resize(_aclPointArray.size());
682
    pDIter = aulDecrements.begin();
683
    ulDec = 0;
684
    pPEnd = _aclPointArray.end();
685
    for (pPIter = _aclPointArray.begin(); pPIter != pPEnd; ++pPIter) {
686
        *pDIter++ = ulDec;
687
        if (!pPIter->IsValid()) {
688
            ulDec++;
689
        }
690
    }
691

692
    // correct point indices of the facets
693
    pFEnd = _aclFacetArray.end();
694
    for (pFIter = _aclFacetArray.begin(); pFIter != pFEnd; ++pFIter) {
695
        if (pFIter->IsValid()) {
696
            pFIter->_aulPoints[0] -= aulDecrements[pFIter->_aulPoints[0]];
697
            pFIter->_aulPoints[1] -= aulDecrements[pFIter->_aulPoints[1]];
698
            pFIter->_aulPoints[2] -= aulDecrements[pFIter->_aulPoints[2]];
699
        }
700
    }
701

702
    // delete point, number of valid points
703
    unsigned long ulNewPts =
704
        std::count_if(_aclPointArray.begin(), _aclPointArray.end(), [](const MeshPoint& p) {
705
            return p.IsValid();
706
        });
707
    // tmp. point array
708
    MeshPointArray aclTempPt(ulNewPts);
709
    MeshPointArray::_TIterator pPTemp = aclTempPt.begin();
710
    pPEnd = _aclPointArray.end();
711
    for (pPIter = _aclPointArray.begin(); pPIter != pPEnd; ++pPIter) {
712
        if (pPIter->IsValid()) {
713
            *pPTemp++ = *pPIter;
714
        }
715
    }
716

717
    // free memory
718
    //_aclPointArray = aclTempPt;
719
    // aclTempPt.clear();
720
    _aclPointArray.swap(aclTempPt);
721
    MeshPointArray().swap(aclTempPt);
722

723
    // generate array of facet decrements
724
    aulDecrements.resize(_aclFacetArray.size());
725
    pDIter = aulDecrements.begin();
726
    ulDec = 0;
727
    pFEnd = _aclFacetArray.end();
728
    for (pFIter = _aclFacetArray.begin(); pFIter != pFEnd; ++pFIter, ++pDIter) {
729
        *pDIter = ulDec;
730
        if (!pFIter->IsValid()) {
731
            ulDec++;
732
        }
733
    }
734

735
    // correct neighbour indices of the facets
736
    pFEnd = _aclFacetArray.end();
737
    for (pFIter = _aclFacetArray.begin(); pFIter != pFEnd; ++pFIter) {
738
        if (pFIter->IsValid()) {
739
            for (FacetIndex& nbIndex : pFIter->_aulNeighbours) {
740
                FacetIndex k = nbIndex;
741
                if (k != FACET_INDEX_MAX) {
742
                    if (_aclFacetArray[k].IsValid()) {
743
                        nbIndex -= aulDecrements[k];
744
                    }
745
                    else {
746
                        nbIndex = FACET_INDEX_MAX;
747
                    }
748
                }
749
            }
750
        }
751
    }
752

753
    // delete facets, number of valid facets
754
    unsigned long ulDelFacets =
755
        std::count_if(_aclFacetArray.begin(), _aclFacetArray.end(), [](const MeshFacet& f) {
756
            return f.IsValid();
757
        });
758
    MeshFacetArray aclFArray(ulDelFacets);
759
    MeshFacetArray::_TIterator pFTemp = aclFArray.begin();
760
    pFEnd = _aclFacetArray.end();
761
    for (pFIter = _aclFacetArray.begin(); pFIter != pFEnd; ++pFIter) {
762
        if (pFIter->IsValid()) {
763
            *pFTemp++ = *pFIter;
764
        }
765
    }
766

767
    // free memory
768
    //_aclFacetArray = aclFArray;
769
    _aclFacetArray.swap(aclFArray);
770
}
771

772
void MeshKernel::CutFacets(const MeshFacetGrid& rclGrid,
773
                           const Base::ViewProjMethod* pclProj,
774
                           const Base::Polygon2d& rclPoly,
775
                           bool bCutInner,
776
                           std::vector<MeshGeomFacet>& raclFacets)
777
{
778
    std::vector<FacetIndex> aulFacets;
779

780
    MeshAlgorithm(*this).CheckFacets(rclGrid, pclProj, rclPoly, bCutInner, aulFacets);
781

782
    for (FacetIndex it : aulFacets) {
783
        raclFacets.push_back(GetFacet(it));
784
    }
785

786
    DeleteFacets(aulFacets);
787
}
788

789
void MeshKernel::CutFacets(const MeshFacetGrid& rclGrid,
790
                           const Base::ViewProjMethod* pclProj,
791
                           const Base::Polygon2d& rclPoly,
792
                           bool bInner,
793
                           std::vector<FacetIndex>& raclCutted)
794
{
795
    MeshAlgorithm(*this).CheckFacets(rclGrid, pclProj, rclPoly, bInner, raclCutted);
796
    DeleteFacets(raclCutted);
797
}
798

799
std::vector<PointIndex> MeshKernel::GetFacetPoints(const std::vector<FacetIndex>& facets) const
800
{
801
    std::vector<PointIndex> points;
802
    for (FacetIndex it : facets) {
803
        PointIndex p0 {}, p1 {}, p2 {};
804
        GetFacetPoints(it, p0, p1, p2);
805
        points.push_back(p0);
806
        points.push_back(p1);
807
        points.push_back(p2);
808
    }
809

810
    std::sort(points.begin(), points.end());
811
    points.erase(std::unique(points.begin(), points.end()), points.end());
812
    return points;
813
}
814

815
std::vector<FacetIndex> MeshKernel::GetPointFacets(const std::vector<PointIndex>& points) const
816
{
817
    _aclPointArray.ResetFlag(MeshPoint::TMP0);
818
    _aclFacetArray.ResetFlag(MeshFacet::TMP0);
819
    for (PointIndex point : points) {
820
        _aclPointArray[point].SetFlag(MeshPoint::TMP0);
821
    }
822

823
    // mark facets if at least one corner point is marked
824
    for (const auto& pF : _aclFacetArray) {
825
        const MeshPoint& rclP0 = _aclPointArray[pF._aulPoints[0]];
826
        const MeshPoint& rclP1 = _aclPointArray[pF._aulPoints[1]];
827
        const MeshPoint& rclP2 = _aclPointArray[pF._aulPoints[2]];
828

829
        if (rclP0.IsFlag(MeshPoint::TMP0) || rclP1.IsFlag(MeshPoint::TMP0)
830
            || rclP2.IsFlag(MeshPoint::TMP0)) {
831
            pF.SetFlag(MeshFacet::TMP0);
832
        }
833
    }
834

835
    std::vector<FacetIndex> facets;
836
    MeshAlgorithm(*this).GetFacetsFlag(facets, MeshFacet::TMP0);
837
    return facets;
838
}
839

840
std::vector<FacetIndex> MeshKernel::HasFacets(const MeshPointIterator& rclIter) const
841
{
842
    PointIndex ulPtInd = rclIter.Position();
843
    std::vector<MeshFacet>::const_iterator pFIter = _aclFacetArray.begin();
844
    std::vector<MeshFacet>::const_iterator pFBegin = _aclFacetArray.begin();
845
    std::vector<MeshFacet>::const_iterator pFEnd = _aclFacetArray.end();
846
    std::vector<FacetIndex> aulBelongs;
847

848
    while (pFIter < pFEnd) {
849
        for (PointIndex point : pFIter->_aulPoints) {
850
            if (point == ulPtInd) {
851
                aulBelongs.push_back(pFIter - pFBegin);
852
                break;
853
            }
854
        }
855
        ++pFIter;
856
    }
857

858
    return aulBelongs;
859
}
860

861
MeshPointArray MeshKernel::GetPoints(const std::vector<PointIndex>& indices) const
862
{
863
    MeshPointArray ary;
864
    ary.reserve(indices.size());
865
    for (PointIndex it : indices) {
866
        ary.push_back(this->_aclPointArray[it]);
867
    }
868
    return ary;
869
}
870

871
MeshFacetArray MeshKernel::GetFacets(const std::vector<FacetIndex>& indices) const
872
{
873
    MeshFacetArray ary;
874
    ary.reserve(indices.size());
875
    for (FacetIndex it : indices) {
876
        ary.push_back(this->_aclFacetArray[it]);
877
    }
878
    return ary;
879
}
880

881
void MeshKernel::Write(std::ostream& rclOut) const
882
{
883
    if (!rclOut || rclOut.bad()) {
884
        return;
885
    }
886

887
    Base::OutputStream str(rclOut);
888

889
    // Write a header with a "magic number" and a version
890
    str << static_cast<uint32_t>(0xA0B0C0D0);
891
    str << static_cast<uint32_t>(0x010000);
892

893
    char szInfo[257];  // needs an additional byte for zero-termination
894
    strcpy(szInfo,
895
           "MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-"
896
           "MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-"
897
           "MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-MESH-"
898
           "MESH-MESH-MESH-\n");
899
    rclOut.write(szInfo, 256);
900

901
    // write the number of points and facets
902
    str << static_cast<uint32_t>(CountPoints()) << static_cast<uint32_t>(CountFacets());
903

904
    // write the data
905
    for (const auto& it : _aclPointArray) {
906
        str << it.x << it.y << it.z;
907
    }
908

909
    for (const auto& it : _aclFacetArray) {
910
        str << static_cast<uint32_t>(it._aulPoints[0]) << static_cast<uint32_t>(it._aulPoints[1])
911
            << static_cast<uint32_t>(it._aulPoints[2]);
912
        str << static_cast<uint32_t>(it._aulNeighbours[0])
913
            << static_cast<uint32_t>(it._aulNeighbours[1])
914
            << static_cast<uint32_t>(it._aulNeighbours[2]);
915
    }
916

917
    str << _clBoundBox.MinX << _clBoundBox.MaxX;
918
    str << _clBoundBox.MinY << _clBoundBox.MaxY;
919
    str << _clBoundBox.MinZ << _clBoundBox.MaxZ;
920
}
921

922
void MeshKernel::Read(std::istream& rclIn)
923
{
924
    if (!rclIn || rclIn.bad()) {
925
        return;
926
    }
927

928
    // get header
929
    Base::InputStream str(rclIn);
930

931
    // Read the header with a "magic number" and a version
932
    uint32_t magic {}, version {}, swap_magic {}, swap_version {};
933
    str >> magic >> version;
934
    swap_magic = magic;
935
    Base::SwapEndian(swap_magic);
936
    swap_version = version;
937
    Base::SwapEndian(swap_version);
938
    uint32_t open_edge = 0xffffffff;  // value to mark an open edge
939

940
    // is it the new or old format?
941
    bool new_format = false;
942
    if (magic == 0xA0B0C0D0 && version == 0x010000) {
943
        new_format = true;
944
    }
945
    else if (swap_magic == 0xA0B0C0D0 && swap_version == 0x010000) {
946
        new_format = true;
947
        str.setByteOrder(Base::Stream::BigEndian);
948
    }
949

950
    if (new_format) {
951
        char szInfo[256];
952
        rclIn.read(szInfo, 256);
953

954
        // read the number of points and facets
955
        uint32_t uCtPts = 0, uCtFts = 0;
956
        str >> uCtPts >> uCtFts;
957

958
        try {
959
            // read the data
960
            MeshPointArray pointArray;
961
            pointArray.resize(uCtPts);
962
            for (auto& it : pointArray) {
963
                str >> it.x >> it.y >> it.z;
964
            }
965

966
            MeshFacetArray facetArray;
967
            facetArray.resize(uCtFts);
968

969
            uint32_t v1 {}, v2 {}, v3 {};
970
            for (auto& it : facetArray) {
971
                str >> v1 >> v2 >> v3;
972

973
                // make sure to have valid indices
974
                if (v1 >= uCtPts || v2 >= uCtPts || v3 >= uCtPts) {
975
                    throw Base::BadFormatError("Invalid data structure");
976
                }
977

978
                it._aulPoints[0] = v1;
979
                it._aulPoints[1] = v2;
980
                it._aulPoints[2] = v3;
981

982
                // On systems where an 'unsigned long' is a 64-bit value
983
                // the empty neighbour must be explicitly set to 'FACET_INDEX_MAX'
984
                // because in algorithms this value is always used to check
985
                // for open edges.
986
                str >> v1 >> v2 >> v3;
987

988
                // make sure to have valid indices
989
                if (v1 >= uCtFts && v1 < open_edge) {
990
                    throw Base::BadFormatError("Invalid data structure");
991
                }
992
                if (v2 >= uCtFts && v2 < open_edge) {
993
                    throw Base::BadFormatError("Invalid data structure");
994
                }
995
                if (v3 >= uCtFts && v3 < open_edge) {
996
                    throw Base::BadFormatError("Invalid data structure");
997
                }
998

999
                if (v1 < open_edge) {
1000
                    it._aulNeighbours[0] = v1;
1001
                }
1002
                else {
1003
                    it._aulNeighbours[0] = FACET_INDEX_MAX;
1004
                }
1005

1006
                if (v2 < open_edge) {
1007
                    it._aulNeighbours[1] = v2;
1008
                }
1009
                else {
1010
                    it._aulNeighbours[1] = FACET_INDEX_MAX;
1011
                }
1012

1013
                if (v3 < open_edge) {
1014
                    it._aulNeighbours[2] = v3;
1015
                }
1016
                else {
1017
                    it._aulNeighbours[2] = FACET_INDEX_MAX;
1018
                }
1019
            }
1020

1021
            str >> _clBoundBox.MinX >> _clBoundBox.MaxX;
1022
            str >> _clBoundBox.MinY >> _clBoundBox.MaxY;
1023
            str >> _clBoundBox.MinZ >> _clBoundBox.MaxZ;
1024

1025
            // If we reach this block no exception occurred and we can safely assign the mesh
1026
            _aclPointArray.swap(pointArray);
1027
            _aclFacetArray.swap(facetArray);
1028
        }
1029
        catch (std::exception&) {
1030
            // Special handling of std::length_error
1031
            throw Base::BadFormatError("Reading from stream failed");
1032
        }
1033
    }
1034
    else {
1035
        // The old formats
1036
        unsigned long uCtPts = magic, uCtFts = version;
1037
        MeshPointArray pointArray;
1038
        MeshFacetArray facetArray;
1039

1040
        float ratio = 0;
1041
        if (uCtPts > 0) {
1042
            ratio = static_cast<float>(uCtFts) / static_cast<float>(uCtPts);
1043
        }
1044

1045
        // without edge array
1046
        if (ratio < 2.5f) {
1047
            // the stored mesh kernel might be empty
1048
            if (uCtPts > 0) {
1049
                pointArray.resize(uCtPts);
1050
                rclIn.read((char*)&(pointArray[0]), uCtPts * sizeof(MeshPoint));
1051
            }
1052
            if (uCtFts > 0) {
1053
                facetArray.resize(uCtFts);
1054
                rclIn.read((char*)&(facetArray[0]), uCtFts * sizeof(MeshFacet));
1055
            }
1056
            rclIn.read((char*)&_clBoundBox, sizeof(Base::BoundBox3f));
1057
        }
1058
        else {
1059
            // with edge array
1060
            unsigned long uCtEdges = uCtFts;
1061
            str >> magic;
1062
            uCtFts = magic;
1063
            pointArray.resize(uCtPts);
1064
            for (auto& it : pointArray) {
1065
                str >> it.x >> it.y >> it.z;
1066
            }
1067
            uint32_t dummy {};
1068
            for (unsigned long i = 0; i < uCtEdges; i++) {
1069
                str >> dummy;
1070
            }
1071
            uint32_t v1 {}, v2 {}, v3 {};
1072
            facetArray.resize(uCtFts);
1073
            for (auto& it : facetArray) {
1074
                str >> v1 >> v2 >> v3;
1075
                it._aulNeighbours[0] = v1;
1076
                it._aulNeighbours[1] = v2;
1077
                it._aulNeighbours[2] = v3;
1078
                str >> v1 >> v2 >> v3;
1079
                it._aulPoints[0] = v1;
1080
                it._aulPoints[1] = v2;
1081
                it._aulPoints[2] = v3;
1082
                str >> it._ucFlag;
1083
            }
1084

1085
            str >> _clBoundBox.MinX >> _clBoundBox.MinY >> _clBoundBox.MinZ >> _clBoundBox.MaxX
1086
                >> _clBoundBox.MaxY >> _clBoundBox.MaxZ;
1087
        }
1088

1089
        for (auto& it : facetArray) {
1090
            for (int i = 0; i < 3; i++) {
1091
                if (it._aulPoints[i] >= uCtPts) {
1092
                    throw Base::BadFormatError("Invalid data structure");
1093
                }
1094
                if (it._aulNeighbours[i] < FACET_INDEX_MAX && it._aulNeighbours[i] >= uCtFts) {
1095
                    throw Base::BadFormatError("Invalid data structure");
1096
                }
1097
            }
1098
        }
1099

1100
        _aclPointArray.swap(pointArray);
1101
        _aclFacetArray.swap(facetArray);
1102
    }
1103
}
1104

1105
void MeshKernel::operator*=(const Base::Matrix4D& rclMat)
1106
{
1107
    this->Transform(rclMat);
1108
}
1109

1110
void MeshKernel::Transform(const Base::Matrix4D& rclMat)
1111
{
1112
    MeshPointArray::_TIterator clPIter = _aclPointArray.begin(), clPEIter = _aclPointArray.end();
1113
    Base::Matrix4D clMatrix(rclMat);
1114

1115
    _clBoundBox.SetVoid();
1116
    while (clPIter < clPEIter) {
1117
        *clPIter *= clMatrix;
1118
        _clBoundBox.Add(*clPIter);
1119
        clPIter++;
1120
    }
1121
}
1122

1123
void MeshKernel::Smooth(int iterations, float stepsize)
1124
{
1125
    (void)stepsize;
1126
    LaplaceSmoothing(*this).Smooth(iterations);
1127
}
1128

1129
void MeshKernel::RecalcBoundBox() const
1130
{
1131
    _clBoundBox.SetVoid();
1132
    for (const auto& pI : _aclPointArray) {
1133
        _clBoundBox.Add(pI);
1134
    }
1135
}
1136

1137
std::vector<Base::Vector3f> MeshKernel::CalcVertexNormals() const
1138
{
1139
    std::vector<Base::Vector3f> normals;
1140

1141
    normals.resize(CountPoints());
1142

1143
    PointIndex p1 {}, p2 {}, p3 {};
1144
    unsigned int ct = CountFacets();
1145
    for (unsigned int pFIter = 0; pFIter < ct; pFIter++) {
1146
        GetFacetPoints(pFIter, p1, p2, p3);
1147
        Base::Vector3f Norm = (GetPoint(p2) - GetPoint(p1)) % (GetPoint(p3) - GetPoint(p1));
1148

1149
        normals[p1] += Norm;
1150
        normals[p2] += Norm;
1151
        normals[p3] += Norm;
1152
    }
1153

1154
    return normals;
1155
}
1156

1157
std::vector<Base::Vector3f> MeshKernel::GetFacetNormals(const std::vector<FacetIndex>& facets) const
1158
{
1159
    std::vector<Base::Vector3f> normals;
1160
    normals.reserve(facets.size());
1161

1162
    for (FacetIndex it : facets) {
1163
        const MeshFacet& face = _aclFacetArray[it];
1164

1165
        const Base::Vector3f& p1 = _aclPointArray[face._aulPoints[0]];
1166
        const Base::Vector3f& p2 = _aclPointArray[face._aulPoints[1]];
1167
        const Base::Vector3f& p3 = _aclPointArray[face._aulPoints[2]];
1168

1169
        Base::Vector3f n = (p2 - p1) % (p3 - p1);
1170
        n.Normalize();
1171
        normals.emplace_back(n);
1172
    }
1173

1174
    return normals;
1175
}
1176

1177
// Evaluation
1178
float MeshKernel::GetSurface() const
1179
{
1180
    float fSurface = 0.0;
1181
    MeshFacetIterator cIter(*this);
1182
    for (cIter.Init(); cIter.More(); cIter.Next()) {
1183
        fSurface += cIter->Area();
1184
    }
1185

1186
    return fSurface;
1187
}
1188

1189
float MeshKernel::GetSurface(const std::vector<FacetIndex>& aSegment) const
1190
{
1191
    float fSurface = 0.0;
1192
    MeshFacetIterator cIter(*this);
1193

1194
    for (FacetIndex it : aSegment) {
1195
        cIter.Set(it);
1196
        fSurface += cIter->Area();
1197
    }
1198

1199
    return fSurface;
1200
}
1201

1202
float MeshKernel::GetVolume() const
1203
{
1204
    // MeshEvalSolid cSolid(*this);
1205
    // if ( !cSolid.Evaluate() )
1206
    //     return 0.0f; // no solid
1207

1208
    float fVolume = 0.0;
1209
    MeshFacetIterator cIter(*this);
1210
    Base::Vector3f p1, p2, p3;
1211
    for (cIter.Init(); cIter.More(); cIter.Next()) {
1212
        const MeshGeomFacet& rclF = *cIter;
1213
        p1 = rclF._aclPoints[0];
1214
        p2 = rclF._aclPoints[1];
1215
        p3 = rclF._aclPoints[2];
1216

1217
        fVolume += (-p3.x * p2.y * p1.z + p2.x * p3.y * p1.z + p3.x * p1.y * p2.z
1218
                    - p1.x * p3.y * p2.z - p2.x * p1.y * p3.z + p1.x * p2.y * p3.z);
1219
    }
1220

1221
    fVolume /= 6.0f;
1222
    fVolume = fabs(fVolume);
1223

1224
    return fVolume;
1225
}
1226

1227
bool MeshKernel::HasOpenEdges() const
1228
{
1229
    MeshEvalSolid eval(*this);
1230
    return !eval.Evaluate();
1231
}
1232

1233
bool MeshKernel::HasNonManifolds() const
1234
{
1235
    MeshEvalTopology eval(*this);
1236
    return !eval.Evaluate();
1237
}
1238

1239
bool MeshKernel::HasSelfIntersections() const
1240
{
1241
    MeshEvalSelfIntersection eval(*this);
1242
    return !eval.Evaluate();
1243
}
1244

1245
// Iterators
1246
MeshFacetIterator MeshKernel::FacetIterator() const
1247
{
1248
    MeshFacetIterator it(*this);
1249
    it.Begin();
1250
    return it;
1251
}
1252

1253
MeshPointIterator MeshKernel::PointIterator() const
1254
{
1255
    MeshPointIterator it(*this);
1256
    it.Begin();
1257
    return it;
1258
}
1259

1260
void MeshKernel::GetEdges(std::vector<MeshGeomEdge>& edges) const
1261
{
1262
    std::set<MeshBuilder::Edge> tmp;
1263

1264
    for (const auto& it : _aclFacetArray) {
1265
        for (int i = 0; i < 3; i++) {
1266
            tmp.insert(MeshBuilder::Edge(it._aulPoints[i],
1267
                                         it._aulPoints[(i + 1) % 3],
1268
                                         it._aulNeighbours[i]));
1269
        }
1270
    }
1271

1272
    edges.reserve(tmp.size());
1273
    for (const auto& it2 : tmp) {
1274
        MeshGeomEdge edge;
1275
        edge._aclPoints[0] = this->_aclPointArray[it2.pt1];
1276
        edge._aclPoints[1] = this->_aclPointArray[it2.pt2];
1277
        edge._bBorder = it2.facetIdx == FACET_INDEX_MAX;
1278

1279
        edges.push_back(edge);
1280
    }
1281
}
1282

1283
unsigned long MeshKernel::CountEdges() const
1284
{
1285
    unsigned long openEdges = 0, closedEdges = 0;
1286

1287
    for (const auto& it : _aclFacetArray) {
1288
        for (FacetIndex nbFacet : it._aulNeighbours) {
1289
            if (nbFacet == FACET_INDEX_MAX) {
1290
                openEdges++;
1291
            }
1292
            else {
1293
                closedEdges++;
1294
            }
1295
        }
1296
    }
1297

1298
    return (openEdges + (closedEdges / 2));
1299
}
1300