FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2010 Werner Mayer <wmayer[at]users.sourceforge.net>     *
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

27
#include <BRepMesh_IncrementalMesh.hxx>
28
#include <BRepTools.hxx>
29
#include <Standard_Version.hxx>
30
#include <TopoDS_Shape.hxx>
31
#endif
32

33
#include <Base/Console.h>
34
#include <Base/Tools.h>
35
#include <Mod/Mesh/App/Mesh.h>
36
#include <Mod/Part/App/BRepMesh.h>
37
#include <Mod/Part/App/TopoShape.h>
38

39
#include "Mesher.h"
40

41
#ifdef HAVE_SMESH
42
#if defined(__clang__)
43
#pragma clang diagnostic push
44
#pragma clang diagnostic ignored "-Woverloaded-virtual"
45
#pragma clang diagnostic ignored "-Wextra-semi"
46
#elif defined(__GNUC__)
47
#pragma GCC diagnostic push
48
#pragma GCC diagnostic ignored "-Wpedantic"
49
#endif
50

51
#include <SMESHDS_Mesh.hxx>
52
#include <SMESH_Gen.hxx>
53
#include <SMESH_Mesh.hxx>
54
#include <StdMeshers_MaxLength.hxx>
55

56
#include <StdMeshers_Arithmetic1D.hxx>
57
#include <StdMeshers_AutomaticLength.hxx>
58
#include <StdMeshers_Deflection1D.hxx>
59
#include <StdMeshers_LocalLength.hxx>
60
#if SMESH_VERSION_MAJOR <= 9 && SMESH_VERSION_MINOR < 10
61
#include <StdMeshers_MEFISTO_2D.hxx>
62
#endif
63
#include <StdMeshers_MaxElementArea.hxx>
64
#include <StdMeshers_NumberOfSegments.hxx>
65
#include <StdMeshers_QuadranglePreference.hxx>
66
#include <StdMeshers_Quadrangle_2D.hxx>
67
#include <StdMeshers_Regular_1D.hxx>
68

69
#include <StdMeshers_LengthFromEdges.hxx>
70
#include <StdMeshers_NotConformAllowed.hxx>
71
#if defined(HAVE_NETGEN)
72
#include <NETGENPlugin_Hypothesis_2D.hxx>
73
#include <NETGENPlugin_NETGEN_2D.hxx>
74
#include <NETGENPlugin_SimpleHypothesis_2D.hxx>
75
#endif  // HAVE_NETGEN
76
#if defined(__clang__)
77
#pragma clang diagnostic pop
78
#elif defined(__GNUC__)
79
#pragma GCC diagnostic pop
80
#endif
81
#endif  // HAVE_SMESH
82

83
using namespace MeshPart;
84

85
SMESH_Gen* Mesher::_mesh_gen = nullptr;
86

87

88
MeshingOutput::MeshingOutput()
89
{
90
    buffer.reserve(80);
91
}
92

93
int MeshingOutput::overflow(int c)
94
{
95
    if (c != EOF) {
96
        buffer.push_back((char)c);
97
    }
98
    return c;
99
}
100

101
int MeshingOutput::sync()
102
{
103
    // Print as log as this might be verbose
104
    if (!buffer.empty()) {
105
        if (buffer.find("failed") != std::string::npos) {
106
            std::string::size_type pos = buffer.find(" : ");
107
            std::string sub;
108
            if (pos != std::string::npos) {
109
                // chop the last newline
110
                sub = buffer.substr(pos + 3, buffer.size() - pos - 4);
111
            }
112
            else {
113
                sub = buffer;
114
            }
115
            Base::Console().Error("%s", sub.c_str());
116
        }
117
        buffer.clear();
118
    }
119
    return 0;
120
}
121

122
// ----------------------------------------------------------------------------
123

124
namespace MeshPart
125
{
126

127
class BrepMesh
128
{
129
    bool segments;
130
    std::vector<uint32_t> colors;
131

132
public:
133
    BrepMesh(bool s, const std::vector<uint32_t>& c)
134
        : segments(s)
135
        , colors(c)
136
    {}
137

138
    Mesh::MeshObject* create(const std::vector<Part::TopoShape::Domain>& domains) const
139
    {
140
        std::vector<Base::Vector3d> points;
141
        std::vector<Part::TopoShape::Facet> facets;
142
        Part::BRepMesh mesh;
143
        mesh.getFacesFromDomains(domains, points, facets);
144

145
        MeshCore::MeshFacetArray faces;
146
        faces.reserve(facets.size());
147
        std::transform(facets.cbegin(),
148
                       facets.cend(),
149
                       std::back_inserter(faces),
150
                       [](const Part::TopoShape::Facet& face) {
151
                           return MeshCore::MeshFacet(face.I1, face.I2, face.I3);
152
                       });
153

154
        MeshCore::MeshPointArray verts;
155
        verts.reserve(points.size());
156
        for (const auto& it : points) {
157
            verts.emplace_back(float(it.x), float(it.y), float(it.z));
158
        }
159

160
        MeshCore::MeshKernel kernel;
161
        kernel.Adopt(verts, faces, true);
162

163
        // mesh segments
164
        std::vector<std::vector<MeshCore::FacetIndex>> meshSegments;
165

166
        std::map<uint32_t, std::vector<std::size_t>> colorMap;
167
        for (std::size_t i = 0; i < colors.size(); i++) {
168
            colorMap[colors[i]].push_back(i);
169
        }
170

171
        bool createSegm = (colors.size() == domains.size());
172

173
        // add a segment for the face
174
        if (createSegm || this->segments) {
175
            auto segments = mesh.createSegments();
176
            meshSegments.reserve(segments.size());
177
            std::transform(segments.cbegin(),
178
                           segments.cend(),
179
                           std::back_inserter(meshSegments),
180
                           [](const Part::BRepMesh::Segment& segm) {
181
                               std::vector<MeshCore::FacetIndex> faces;
182
                               faces.insert(faces.end(), segm.cbegin(), segm.cend());
183
                               return faces;
184
                           });
185
        }
186

187
        Mesh::MeshObject* meshdata = new Mesh::MeshObject();
188
        meshdata->swap(kernel);
189
        if (createSegm) {
190
            int index = 0;
191
            for (const auto& it : colorMap) {
192
                Mesh::Segment segm(meshdata, false);
193
                for (auto jt : it.second) {
194
                    segm.addIndices(meshSegments[jt]);
195
                }
196
                segm.save(true);
197
                std::stringstream str;
198
                str << "patch" << index++;
199
                segm.setName(str.str());
200
                App::Color col;
201
                col.setPackedValue(it.first);
202
                segm.setColor(col.asHexString());
203
                meshdata->addSegment(segm);
204
            }
205
        }
206
        else {
207
            for (const auto& it : meshSegments) {
208
                meshdata->addSegment(it);
209
            }
210
        }
211
        return meshdata;
212
    }
213
};
214
}  // namespace MeshPart
215

216
// ----------------------------------------------------------------------------
217

218
Mesher::Mesher(const TopoDS_Shape& s)
219
    : shape(s)
220
{}
221

222
Mesher::~Mesher() = default;
223

224
Mesh::MeshObject* Mesher::createStandard() const
225
{
226
    if (!shape.IsNull()) {
227
        BRepTools::Clean(shape);
228
        BRepMesh_IncrementalMesh aMesh(shape, deflection, relative, angularDeflection);
229
    }
230

231
    std::vector<Part::TopoShape::Domain> domains;
232
    Part::TopoShape(shape).getDomains(domains);
233

234
    BrepMesh brepmesh(this->segments, this->colors);
235
    return brepmesh.create(domains);
236
}
237

238
Mesh::MeshObject* Mesher::createMesh() const
239
{
240
    // OCC standard mesher
241
    if (method == Standard) {
242
        return createStandard();
243
    }
244

245
#ifndef HAVE_SMESH
246
    throw Base::RuntimeError("SMESH is not available on this platform");
247
#else
248
    std::list<SMESH_Hypothesis*> hypoth;
249

250
    if (!Mesher::_mesh_gen) {
251
        Mesher::_mesh_gen = new SMESH_Gen();
252
    }
253
    SMESH_Gen* meshgen = Mesher::_mesh_gen;
254

255
#if SMESH_VERSION_MAJOR >= 9
256
    SMESH_Mesh* mesh = meshgen->CreateMesh(true);
257
#else
258
    SMESH_Mesh* mesh = meshgen->CreateMesh(0, true);
259
#endif
260

261
    int hyp = 0;
262

263
    switch (method) {
264
#if defined(HAVE_NETGEN)
265
        case Netgen: {
266
#if SMESH_VERSION_MAJOR >= 9
267
            NETGENPlugin_Hypothesis_2D* hyp2d = new NETGENPlugin_Hypothesis_2D(hyp++, meshgen);
268
#else
269
            NETGENPlugin_Hypothesis_2D* hyp2d = new NETGENPlugin_Hypothesis_2D(hyp++, 0, meshgen);
270
#endif
271

272
            if (fineness >= 0 && fineness < 5) {
273
                hyp2d->SetFineness(NETGENPlugin_Hypothesis_2D::Fineness(fineness));
274
            }
275
            // user defined values
276
            else {
277
                if (growthRate > 0) {
278
                    hyp2d->SetGrowthRate(growthRate);
279
                }
280
                if (nbSegPerEdge > 0) {
281
                    hyp2d->SetNbSegPerEdge(nbSegPerEdge);
282
                }
283
                if (nbSegPerRadius > 0) {
284
                    hyp2d->SetNbSegPerRadius(nbSegPerRadius);
285
                }
286
            }
287

288
            if (maxLen > 0) {
289
                hyp2d->SetMaxSize(maxLen);
290
            }
291
            if (minLen > 0) {
292
                hyp2d->SetMinSize(minLen);
293
            }
294

295
            hyp2d->SetQuadAllowed(allowquad);
296
            hyp2d->SetOptimize(optimize);
297
            hyp2d->SetSecondOrder(
298
                secondOrder);  // apply bisecting to create four triangles out of one
299
            hypoth.push_back(hyp2d);
300

301
#if SMESH_VERSION_MAJOR >= 9
302
            NETGENPlugin_NETGEN_2D* alg2d = new NETGENPlugin_NETGEN_2D(hyp++, meshgen);
303
#else
304
            NETGENPlugin_NETGEN_2D* alg2d = new NETGENPlugin_NETGEN_2D(hyp++, 0, meshgen);
305
#endif
306
            hypoth.push_back(alg2d);
307
        } break;
308
#endif
309
#if SMESH_VERSION_MAJOR <= 9 && SMESH_VERSION_MINOR < 10
310
#if defined(HAVE_MEFISTO)
311
        case Mefisto: {
312
            if (maxLength > 0) {
313
#if SMESH_VERSION_MAJOR >= 9
314
                StdMeshers_MaxLength* hyp1d = new StdMeshers_MaxLength(hyp++, meshgen);
315
#else
316
                StdMeshers_MaxLength* hyp1d = new StdMeshers_MaxLength(hyp++, 0, meshgen);
317
#endif
318
                hyp1d->SetLength(maxLength);
319
                hypoth.push_back(hyp1d);
320
            }
321
            else if (localLength > 0) {
322
#if SMESH_VERSION_MAJOR >= 9
323
                StdMeshers_LocalLength* hyp1d = new StdMeshers_LocalLength(hyp++, meshgen);
324
#else
325
                StdMeshers_LocalLength* hyp1d = new StdMeshers_LocalLength(hyp++, 0, meshgen);
326
#endif
327
                hyp1d->SetLength(localLength);
328
                hypoth.push_back(hyp1d);
329
            }
330
            else if (maxArea > 0) {
331
#if SMESH_VERSION_MAJOR >= 9
332
                StdMeshers_MaxElementArea* hyp2d = new StdMeshers_MaxElementArea(hyp++, meshgen);
333
#else
334
                StdMeshers_MaxElementArea* hyp2d = new StdMeshers_MaxElementArea(hyp++, 0, meshgen);
335
#endif
336
                hyp2d->SetMaxArea(maxArea);
337
                hypoth.push_back(hyp2d);
338
            }
339
            else if (deflection > 0) {
340
#if SMESH_VERSION_MAJOR >= 9
341
                StdMeshers_Deflection1D* hyp1d = new StdMeshers_Deflection1D(hyp++, meshgen);
342
#else
343
                StdMeshers_Deflection1D* hyp1d = new StdMeshers_Deflection1D(hyp++, 0, meshgen);
344
#endif
345
                hyp1d->SetDeflection(deflection);
346
                hypoth.push_back(hyp1d);
347
            }
348
            else if (minLen > 0 && maxLen > 0) {
349
#if SMESH_VERSION_MAJOR >= 9
350
                StdMeshers_Arithmetic1D* hyp1d = new StdMeshers_Arithmetic1D(hyp++, meshgen);
351
#else
352
                StdMeshers_Arithmetic1D* hyp1d = new StdMeshers_Arithmetic1D(hyp++, 0, meshgen);
353
#endif
354
                hyp1d->SetLength(minLen, false);
355
                hyp1d->SetLength(maxLen, true);
356
                hypoth.push_back(hyp1d);
357
            }
358
            else {
359
#if SMESH_VERSION_MAJOR >= 9
360
                StdMeshers_AutomaticLength* hyp1d = new StdMeshers_AutomaticLength(hyp++, meshgen);
361
#else
362
                StdMeshers_AutomaticLength* hyp1d =
363
                    new StdMeshers_AutomaticLength(hyp++, 0, meshgen);
364
#endif
365
                hypoth.push_back(hyp1d);
366
            }
367

368
            {
369
#if SMESH_VERSION_MAJOR >= 9
370
                StdMeshers_NumberOfSegments* hyp1d =
371
                    new StdMeshers_NumberOfSegments(hyp++, meshgen);
372
#else
373
                StdMeshers_NumberOfSegments* hyp1d =
374
                    new StdMeshers_NumberOfSegments(hyp++, 0, meshgen);
375
#endif
376
                hyp1d->SetNumberOfSegments(1);
377
                hypoth.push_back(hyp1d);
378
            }
379

380
            if (regular) {
381
#if SMESH_VERSION_MAJOR >= 9
382
                StdMeshers_Regular_1D* hyp1d = new StdMeshers_Regular_1D(hyp++, meshgen);
383
#else
384
                StdMeshers_Regular_1D* hyp1d = new StdMeshers_Regular_1D(hyp++, 0, meshgen);
385
#endif
386
                hypoth.push_back(hyp1d);
387
            }
388

389
#if SMESH_VERSION_MAJOR >= 9
390
            StdMeshers_MEFISTO_2D* alg2d = new StdMeshers_MEFISTO_2D(hyp++, meshgen);
391
#else
392
            StdMeshers_MEFISTO_2D* alg2d = new StdMeshers_MEFISTO_2D(hyp++, 0, meshgen);
393
#endif
394
            hypoth.push_back(alg2d);
395
        } break;
396
#endif
397
#endif
398
        default:
399
            break;
400
    }
401

402
    // Set new cout
403
    MeshingOutput stdcout;
404
    std::streambuf* oldcout = std::cout.rdbuf(&stdcout);
405

406
    // Apply the hypothesis and create the mesh
407
    mesh->ShapeToMesh(shape);
408
    for (int i = 0; i < hyp; i++) {
409
        mesh->AddHypothesis(shape, i);
410
    }
411
    meshgen->Compute(*mesh, mesh->GetShapeToMesh());
412

413
    // Restore old cout
414
    std::cout.rdbuf(oldcout);
415

416
    // build up the mesh structure
417
    Mesh::MeshObject* meshdata = createFrom(mesh);
418

419
    // clean up
420
    TopoDS_Shape aNull;
421
    mesh->ShapeToMesh(aNull);
422
    mesh->Clear();
423
    delete mesh;
424
    for (auto it : hypoth) {
425
        delete it;
426
    }
427

428
    return meshdata;
429
#endif  // HAVE_SMESH
430
}
431

432
Mesh::MeshObject* Mesher::createFrom(SMESH_Mesh* mesh) const
433
{
434
    // build up the mesh structure
435
    SMDS_FaceIteratorPtr aFaceIter = mesh->GetMeshDS()->facesIterator();
436
    SMDS_NodeIteratorPtr aNodeIter = mesh->GetMeshDS()->nodesIterator();
437

438
    MeshCore::MeshPointArray verts;
439
    MeshCore::MeshFacetArray faces;
440
    verts.reserve(mesh->NbNodes());
441
    faces.reserve(mesh->NbFaces());
442

443
    int index = 0;
444
    std::map<const SMDS_MeshNode*, int> mapNodeIndex;
445
    for (; aNodeIter->more();) {
446
        const SMDS_MeshNode* aNode = aNodeIter->next();
447
        MeshCore::MeshPoint p;
448
        p.Set((float)aNode->X(), (float)aNode->Y(), (float)aNode->Z());
449
        verts.push_back(p);
450
        mapNodeIndex[aNode] = index++;
451
    }
452

453
    for (; aFaceIter->more();) {
454
        const SMDS_MeshFace* aFace = aFaceIter->next();
455
        if (aFace->NbNodes() == 3) {
456
            MeshCore::MeshFacet f;
457
            for (int i = 0; i < 3; i++) {
458
                const SMDS_MeshNode* node = aFace->GetNode(i);
459
                f._aulPoints[i] = mapNodeIndex[node];
460
            }
461
            faces.push_back(f);
462
        }
463
        else if (aFace->NbNodes() == 4) {
464
            MeshCore::MeshFacet f1, f2;
465
            const SMDS_MeshNode* node0 = aFace->GetNode(0);
466
            const SMDS_MeshNode* node1 = aFace->GetNode(1);
467
            const SMDS_MeshNode* node2 = aFace->GetNode(2);
468
            const SMDS_MeshNode* node3 = aFace->GetNode(3);
469

470
            f1._aulPoints[0] = mapNodeIndex[node0];
471
            f1._aulPoints[1] = mapNodeIndex[node1];
472
            f1._aulPoints[2] = mapNodeIndex[node2];
473

474
            f2._aulPoints[0] = mapNodeIndex[node0];
475
            f2._aulPoints[1] = mapNodeIndex[node2];
476
            f2._aulPoints[2] = mapNodeIndex[node3];
477

478
            faces.push_back(f1);
479
            faces.push_back(f2);
480
        }
481
        else if (aFace->NbNodes() == 6) {
482
            MeshCore::MeshFacet f1, f2, f3, f4;
483
            const SMDS_MeshNode* node0 = aFace->GetNode(0);
484
            const SMDS_MeshNode* node1 = aFace->GetNode(1);
485
            const SMDS_MeshNode* node2 = aFace->GetNode(2);
486
            const SMDS_MeshNode* node3 = aFace->GetNode(3);
487
            const SMDS_MeshNode* node4 = aFace->GetNode(4);
488
            const SMDS_MeshNode* node5 = aFace->GetNode(5);
489

490
            f1._aulPoints[0] = mapNodeIndex[node0];
491
            f1._aulPoints[1] = mapNodeIndex[node3];
492
            f1._aulPoints[2] = mapNodeIndex[node5];
493

494
            f2._aulPoints[0] = mapNodeIndex[node1];
495
            f2._aulPoints[1] = mapNodeIndex[node4];
496
            f2._aulPoints[2] = mapNodeIndex[node3];
497

498
            f3._aulPoints[0] = mapNodeIndex[node2];
499
            f3._aulPoints[1] = mapNodeIndex[node5];
500
            f3._aulPoints[2] = mapNodeIndex[node4];
501

502
            f4._aulPoints[0] = mapNodeIndex[node3];
503
            f4._aulPoints[1] = mapNodeIndex[node4];
504
            f4._aulPoints[2] = mapNodeIndex[node5];
505

506
            faces.push_back(f1);
507
            faces.push_back(f2);
508
            faces.push_back(f3);
509
            faces.push_back(f4);
510
        }
511
        else if (aFace->NbNodes() == 8) {
512
            MeshCore::MeshFacet f1, f2, f3, f4, f5, f6;
513
            const SMDS_MeshNode* node0 = aFace->GetNode(0);
514
            const SMDS_MeshNode* node1 = aFace->GetNode(1);
515
            const SMDS_MeshNode* node2 = aFace->GetNode(2);
516
            const SMDS_MeshNode* node3 = aFace->GetNode(3);
517
            const SMDS_MeshNode* node4 = aFace->GetNode(4);
518
            const SMDS_MeshNode* node5 = aFace->GetNode(5);
519
            const SMDS_MeshNode* node6 = aFace->GetNode(6);
520
            const SMDS_MeshNode* node7 = aFace->GetNode(7);
521

522
            f1._aulPoints[0] = mapNodeIndex[node0];
523
            f1._aulPoints[1] = mapNodeIndex[node4];
524
            f1._aulPoints[2] = mapNodeIndex[node7];
525

526
            f2._aulPoints[0] = mapNodeIndex[node1];
527
            f2._aulPoints[1] = mapNodeIndex[node5];
528
            f2._aulPoints[2] = mapNodeIndex[node4];
529

530
            f3._aulPoints[0] = mapNodeIndex[node2];
531
            f3._aulPoints[1] = mapNodeIndex[node6];
532
            f3._aulPoints[2] = mapNodeIndex[node5];
533

534
            f4._aulPoints[0] = mapNodeIndex[node3];
535
            f4._aulPoints[1] = mapNodeIndex[node7];
536
            f4._aulPoints[2] = mapNodeIndex[node6];
537

538
            // Two solutions are possible:
539
            // <4,6,7>, <4,5,6> or <4,5,7>, <5,6,7>
540
            Base::Vector3d v4(node4->X(), node4->Y(), node4->Z());
541
            Base::Vector3d v5(node5->X(), node5->Y(), node5->Z());
542
            Base::Vector3d v6(node6->X(), node6->Y(), node6->Z());
543
            Base::Vector3d v7(node7->X(), node7->Y(), node7->Z());
544
            double dist46 = Base::DistanceP2(v4, v6);
545
            double dist57 = Base::DistanceP2(v5, v7);
546
            if (dist46 > dist57) {
547
                f5._aulPoints[0] = mapNodeIndex[node4];
548
                f5._aulPoints[1] = mapNodeIndex[node6];
549
                f5._aulPoints[2] = mapNodeIndex[node7];
550

551
                f6._aulPoints[0] = mapNodeIndex[node4];
552
                f6._aulPoints[1] = mapNodeIndex[node5];
553
                f6._aulPoints[2] = mapNodeIndex[node6];
554
            }
555
            else {
556
                f5._aulPoints[0] = mapNodeIndex[node4];
557
                f5._aulPoints[1] = mapNodeIndex[node5];
558
                f5._aulPoints[2] = mapNodeIndex[node7];
559

560
                f6._aulPoints[0] = mapNodeIndex[node5];
561
                f6._aulPoints[1] = mapNodeIndex[node6];
562
                f6._aulPoints[2] = mapNodeIndex[node7];
563
            }
564

565
            faces.push_back(f1);
566
            faces.push_back(f2);
567
            faces.push_back(f3);
568
            faces.push_back(f4);
569
            faces.push_back(f5);
570
            faces.push_back(f6);
571
        }
572
        else {
573
            Base::Console().Warning("Face with %d nodes ignored\n", aFace->NbNodes());
574
        }
575
    }
576

577
    MeshCore::MeshKernel kernel;
578
    kernel.Adopt(verts, faces, true);
579

580
    Mesh::MeshObject* meshdata = new Mesh::MeshObject();
581
    meshdata->swap(kernel);
582
    return meshdata;
583
}
584