FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2011 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 <cassert>
26
# include <BRep_Tool.hxx>
27
# include <BRepAdaptor_Curve.hxx>
28
# include <BRepAdaptor_Surface.hxx>
29
# include <BRepBuilderAPI_MakeEdge.hxx>
30
# include <BRepBuilderAPI_MakeFace.hxx>
31
# include <BRepLProp_SLProps.hxx>
32
# include <BRepMesh_IncrementalMesh.hxx>
33
# include <CSLib.hxx>
34
# include <Geom_BSplineSurface.hxx>
35
# include <Geom_Line.hxx>
36
# include <Geom_Plane.hxx>
37
# include <Geom_Point.hxx>
38
# include <GeomAPI_IntSS.hxx>
39
# include <GeomAPI_ProjectPointOnSurf.hxx>
40
# include <GeomAdaptor_Curve.hxx>
41
# include <GeomLib.hxx>
42
# include <GeomLProp_SLProps.hxx>
43
# include <GeomPlate_BuildPlateSurface.hxx>
44
# include <GeomPlate_CurveConstraint.hxx>
45
# include <GeomPlate_MakeApprox.hxx>
46
# include <GeomPlate_PlateG0Criterion.hxx>
47
# include <GeomPlate_PointConstraint.hxx>
48
# include <gp_Lin.hxx>
49
# include <gp_Pln.hxx>
50
# include <gp_Quaternion.hxx>
51
# include <Poly_Connect.hxx>
52
# include <Poly_Triangulation.hxx>
53
# include <Precision.hxx>
54
# include <Standard_Mutex.hxx>
55
# include <Standard_TypeMismatch.hxx>
56
# include <Standard_Version.hxx>
57
# include <TColStd_ListIteratorOfListOfTransient.hxx>
58
# include <TColStd_ListOfTransient.hxx>
59
# include <TColgp_SequenceOfXY.hxx>
60
# include <TColgp_SequenceOfXYZ.hxx>
61
# include <TopoDS.hxx>
62
# if OCC_VERSION_HEX < 0x070600
63
# include <Adaptor3d_HCurveOnSurface.hxx>
64
# include <GeomAdaptor_HCurve.hxx>
65
# endif
66
#endif
67

68
#include <Base/Exception.h>
69
#include <Base/Vector3D.h>
70

71
#include "Tools.h"
72

73

74
void Part::closestPointsOnLines(const gp_Lin& lin1, const gp_Lin& lin2, gp_Pnt& p1, gp_Pnt& p2)
75
{
76
    // they might be the same point
77
    gp_Vec v1(lin1.Direction());
78
    gp_Vec v2(lin2.Direction());
79
    gp_Vec v3(lin2.Location(), lin1.Location());
80

81
    double a = v1*v1;
82
    double b = v1*v2;
83
    double c = v2*v2;
84
    double d = v1*v3;
85
    double e = v2*v3;
86
    double D = a*c - b*b;
87
    double s, t;
88

89
    // D = (v1 x v2) * (v1 x v2)
90
    if (D < Precision::Angular()){
91
        // the lines are considered parallel
92
        s = 0.0;
93
        t = (b>c ? d/b : e/c);
94
    }
95
    else {
96
        s = (b*e - c*d) / D;
97
        t = (a*e - b*d) / D;
98
    }
99

100
    p1 = lin1.Location().XYZ() + s * v1.XYZ();
101
    p2 = lin2.Location().XYZ() + t * v2.XYZ();
102
}
103

104
bool Part::intersect(const gp_Pln& pln1, const gp_Pln& pln2, gp_Lin& lin)
105
{
106
    bool found = false;
107
    Handle (Geom_Plane) gp1 = new Geom_Plane(pln1);
108
    Handle (Geom_Plane) gp2 = new Geom_Plane(pln2);
109

110
    GeomAPI_IntSS intSS(gp1, gp2, Precision::Confusion());
111
    if (intSS.IsDone()) {
112
        int numSol = intSS.NbLines();
113
        if (numSol > 0) {
114
            Handle(Geom_Curve) curve = intSS.Line(1);
115
            lin = Handle(Geom_Line)::DownCast(curve)->Lin();
116
            found = true;
117
        }
118
    }
119

120
    return found;
121
}
122

123
/*! The objects in \a theBoundaries must be of the type Adaptor3d_HCurveOnSurface or
124
GeomAdaptor_HCurve or Geom_Point indicating type of a constraint. Otherwise an exception
125
Standard_TypeMismatch is thrown.
126

127
If the \a theBoundaries list is empty then Standard_ConstructionError is thrown.
128

129
If the algorithm fails it returns a null surface.
130
\see http://opencascade.blogspot.com/2010/03/surface-modeling-part6.html
131
*/
132
Handle(Geom_Surface)
133
Part::Tools::makeSurface(const TColStd_ListOfTransient &theBoundaries,
134
                         const Standard_Real theTol,
135
                         const Standard_Integer theNbPnts,
136
                         const Standard_Integer theNbIter,
137
                         const Standard_Integer theMaxDeg)
138
{
139
    (void)theTol;
140
    //constants for algorithm
141
    const Standard_Integer aNbIter = theNbIter; //number of algorithm iterations
142
    const Standard_Integer aNbPnts = theNbPnts; //sample points per each constraint
143
    const Standard_Integer aDeg = 3; //requested surface degree ?
144
    const Standard_Integer aMaxDeg = theMaxDeg;
145
    const Standard_Integer aMaxSeg = 10000;
146
    const Standard_Real aTol3d = 1.e-04;
147
    const Standard_Real aTol2d = 1.e-05;
148
    const Standard_Real anAngTol = 1.e-02; //angular
149
    const Standard_Real aCurvTol = 1.e-01; //curvature
150

151
    Handle(Geom_Surface) aRes;
152
    GeomPlate_BuildPlateSurface aPlateBuilder (aDeg, aNbPnts, aNbIter, aTol2d, aTol3d, anAngTol, aCurvTol);
153

154
    TColStd_ListIteratorOfListOfTransient anIt (theBoundaries);
155
    if (anIt.More()) {
156
        int i = 1;
157
        for (; anIt.More(); anIt.Next(), i++) {
158
            const Handle(Standard_Transient)& aCur = anIt.Value();
159
            if (aCur.IsNull()) {
160
                assert (0);
161
                Standard_ConstructionError::Raise ("Tools::makeSurface()");
162
            }
163
#if OCC_VERSION_HEX >= 0x070600
164
            else if (aCur->IsKind (STANDARD_TYPE (Adaptor3d_CurveOnSurface))) {
165
                //G1 constraint
166
                Handle(Adaptor3d_CurveOnSurface) aHCOS (Handle(Adaptor3d_CurveOnSurface)::DownCast (aCur));
167
                Handle (GeomPlate_CurveConstraint) aConst = new GeomPlate_CurveConstraint (aHCOS, 1 /*GeomAbs_G1*/,aNbPnts, aTol3d, anAngTol, aCurvTol);
168
                aPlateBuilder.Add (aConst);
169
            }
170
            else if (aCur->IsKind (STANDARD_TYPE (GeomAdaptor_Curve))) {
171
                //G0 constraint
172
                Handle(GeomAdaptor_Curve) aHC (Handle(GeomAdaptor_Curve)::DownCast (aCur));
173
                Handle (GeomPlate_CurveConstraint) aConst = new GeomPlate_CurveConstraint (aHC, 0 /*GeomAbs_G0*/, aNbPnts, aTol3d);
174
                aPlateBuilder.Add (aConst);
175
            }
176
#else
177
            else if (aCur->IsKind (STANDARD_TYPE (Adaptor3d_HCurveOnSurface))) {
178
                //G1 constraint
179
                Handle(Adaptor3d_HCurveOnSurface) aHCOS (Handle(Adaptor3d_HCurveOnSurface)::DownCast (aCur));
180
                Handle (GeomPlate_CurveConstraint) aConst = new GeomPlate_CurveConstraint (aHCOS, 1 /*GeomAbs_G1*/,aNbPnts, aTol3d, anAngTol, aCurvTol);
181
                aPlateBuilder.Add (aConst);
182
            }
183
            else if (aCur->IsKind (STANDARD_TYPE (GeomAdaptor_HCurve))) {
184
                //G0 constraint
185
                Handle(GeomAdaptor_HCurve) aHC (Handle(GeomAdaptor_HCurve)::DownCast (aCur));
186
                Handle (GeomPlate_CurveConstraint) aConst = new GeomPlate_CurveConstraint (aHC, 0 /*GeomAbs_G0*/, aNbPnts, aTol3d);
187
                aPlateBuilder.Add (aConst);
188
            }
189
#endif
190
            else if (aCur->IsKind (STANDARD_TYPE (Geom_Point))) {
191
                //Point constraint
192
                Handle(Geom_Point) aGP (Handle(Geom_Point)::DownCast (aCur));
193
                Handle(GeomPlate_PointConstraint) aConst = new GeomPlate_PointConstraint(aGP->Pnt(),0);
194
                aPlateBuilder.Add(aConst);
195
            }
196
            else {
197
                Standard_TypeMismatch::Raise ("Tools::makeSurface()");
198
            }
199
        }
200
    }
201
    else {
202
        Standard_ConstructionError::Raise ("Tools::makeSurface()");
203
    }
204

205
    //construct
206
    aPlateBuilder.Perform();
207

208
    if (!aPlateBuilder.IsDone()) {
209
        return aRes;
210
    }
211

212
    const Handle(GeomPlate_Surface)& aPlate = aPlateBuilder.Surface();
213
    //approximation (see BRepFill_Filling - when no initial surface was given)
214
    Standard_Real aDMax = aPlateBuilder.G0Error();
215
    TColgp_SequenceOfXY aS2d;
216
    TColgp_SequenceOfXYZ aS3d;
217
    aPlateBuilder.Disc2dContour (4, aS2d);
218
    aPlateBuilder.Disc3dContour (4, 0, aS3d);
219
    Standard_Real aMax = Max (aTol3d, 10. * aDMax);
220
    GeomPlate_PlateG0Criterion aCriterion (aS2d, aS3d, aMax);
221
    {
222
        //data races in AdvApp2Var used by GeomApprox_Surface, use global mutex
223
        //Standard_Mutex::Sentry aSentry (theBSMutex);
224
        GeomPlate_MakeApprox aMakeApprox (aPlate, aCriterion, aTol3d, aMaxSeg, aMaxDeg);
225
        aRes = aMakeApprox.Surface();
226
    }
227

228
    return aRes;
229
}
230

231
bool Part::Tools::getTriangulation(const TopoDS_Face& face, std::vector<gp_Pnt>& points, std::vector<Poly_Triangle>& facets)
232
{
233
    TopLoc_Location loc;
234
    Handle(Poly_Triangulation) hTria = BRep_Tool::Triangulation(face, loc);
235
    if (hTria.IsNull())
236
        return false;
237

238
    // getting the transformation of the face
239
    gp_Trsf transf;
240
    bool identity = true;
241
    if (!loc.IsIdentity()) {
242
        identity = false;
243
        transf = loc.Transformation();
244
    }
245

246
    // check orientation
247
    TopAbs_Orientation orient = face.Orientation();
248

249
    Standard_Integer nbNodes = hTria->NbNodes();
250
    Standard_Integer nbTriangles = hTria->NbTriangles();
251
#if OCC_VERSION_HEX < 0x070600
252
    const TColgp_Array1OfPnt& nodes = hTria->Nodes();
253
    const Poly_Array1OfTriangle& triangles = hTria->Triangles();
254
#endif
255

256
    points.reserve(nbNodes);
257
    facets.reserve(nbTriangles);
258

259
    // cycling through the poly mesh
260
    //
261
    for (int i = 1; i <= nbNodes; i++) {
262
#if OCC_VERSION_HEX < 0x070600
263
        gp_Pnt p = nodes(i);
264
#else
265
        gp_Pnt p = hTria->Node(i);
266
#endif
267

268
        // transform the vertices to the location of the face
269
        if (!identity) {
270
            p.Transform(transf);
271
        }
272

273
        points.push_back(p);
274
    }
275

276
    for (int i = 1; i <= nbTriangles; i++) {
277
        // Get the triangle
278
        Standard_Integer n1,n2,n3;
279
#if OCC_VERSION_HEX < 0x070600
280
        triangles(i).Get(n1, n2, n3);
281
#else
282
        hTria->Triangle(i).Get(n1, n2, n3);
283
#endif
284
        --n1; --n2; --n3;
285

286
        // change orientation of the triangles
287
        if (orient != TopAbs_FORWARD) {
288
            std::swap(n1, n2);
289
        }
290

291
        facets.emplace_back(n1, n2, n3);
292
    }
293

294
    return true;
295
}
296

297
bool Part::Tools::getPolygonOnTriangulation(const TopoDS_Edge& edge, const TopoDS_Face& face, std::vector<gp_Pnt>& points)
298
{
299
    TopLoc_Location loc;
300
    Handle(Poly_Triangulation) hTria = BRep_Tool::Triangulation(face, loc);
301
    if (hTria.IsNull())
302
        return false;
303

304
    // this holds the indices of the edge's triangulation to the actual points
305
    Handle(Poly_PolygonOnTriangulation) hPoly = BRep_Tool::PolygonOnTriangulation(edge, hTria, loc);
306
    if (hPoly.IsNull())
307
        return false;
308

309
    // getting the transformation of the edge
310
    gp_Trsf transf;
311
    bool identity = true;
312
    if (!loc.IsIdentity()) {
313
        identity = false;
314
        transf = loc.Transformation();
315
    }
316

317
    // getting size and create the array
318
    Standard_Integer nbNodes = hPoly->NbNodes();
319
    points.reserve(nbNodes);
320
    const TColStd_Array1OfInteger& indices = hPoly->Nodes();
321
#if OCC_VERSION_HEX < 0x070600
322
    const TColgp_Array1OfPnt& Nodes = hTria->Nodes();
323
#endif
324

325
    // go through the index array
326
    for (Standard_Integer i = indices.Lower(); i <= indices.Upper(); i++) {
327
#if OCC_VERSION_HEX < 0x070600
328
        gp_Pnt p = Nodes(indices(i));
329
#else
330
        gp_Pnt p = hTria->Node(indices(i));
331
#endif
332
        if (!identity) {
333
            p.Transform(transf);
334
        }
335

336
        points.push_back(p);
337
    }
338

339
    return true;
340
}
341

342
bool Part::Tools::getPolygon3D(const TopoDS_Edge& edge, std::vector<gp_Pnt>& points)
343
{
344
    TopLoc_Location loc;
345
    Handle(Poly_Polygon3D) hPoly = BRep_Tool::Polygon3D(edge, loc);
346
    if (hPoly.IsNull())
347
        return false;
348

349
    // getting the transformation of the edge
350
    gp_Trsf transf;
351
    bool identity = true;
352
    if (!loc.IsIdentity()) {
353
        identity = false;
354
        transf = loc.Transformation();
355
    }
356

357
    // getting size and create the array
358
    Standard_Integer nbNodes = hPoly->NbNodes();
359
    points.reserve(nbNodes);
360
    const TColgp_Array1OfPnt& nodes = hPoly->Nodes();
361

362
    for (int i = 1; i <= nbNodes; i++) {
363
        gp_Pnt p = nodes(i);
364

365
        // transform the vertices to the location of the face
366
        if (!identity) {
367
            p.Transform(transf);
368
        }
369

370
        points.push_back(p);
371
    }
372

373
    return true;
374
}
375

376
void Part::Tools::getPointNormals(const std::vector<gp_Pnt>& points, const std::vector<Poly_Triangle>& facets, std::vector<gp_Vec>& vertexnormals)
377
{
378
    vertexnormals.resize(points.size());
379

380
    for (const auto& it : facets) {
381
        // Get the triangle
382
        Standard_Integer n1,n2,n3;
383
        it.Get(n1,n2,n3);
384

385
        // Calculate triangle normal
386
        gp_Vec v1(points[n1].XYZ());
387
        gp_Vec v2(points[n2].XYZ());
388
        gp_Vec v3(points[n3].XYZ());
389
        gp_Vec n = (v2 - v1) ^ (v3 - v1);
390

391
        // add the triangle normal to the vertex normal for all points of this triangle
392
        vertexnormals[n1] += n;
393
        vertexnormals[n2] += n;
394
        vertexnormals[n3] += n;
395
    }
396

397
    for (auto& it : vertexnormals)
398
        it.Normalize();
399
}
400

401
void Part::Tools::getPointNormals(const std::vector<gp_Pnt>& points, const TopoDS_Face& face, std::vector<gp_Vec>& vertexnormals)
402
{
403
    if (points.size() != vertexnormals.size())
404
        return;
405

406
    Handle(Geom_Surface) hSurface = BRep_Tool::Surface(face);
407
    if (hSurface.IsNull())
408
        return;
409

410
    // normalize all vertex normals
411
    for (std::size_t i = 0; i < points.size(); i++) {
412
        try {
413
            GeomAPI_ProjectPointOnSurf ProPntSrf(points[i], hSurface);
414
            Standard_Real u, v;
415
            ProPntSrf.Parameters(1, u, v);
416

417
            GeomLProp_SLProps propOfFace(hSurface, u, v, 2, gp::Resolution());
418

419
            gp_Dir normal = propOfFace.Normal();
420
            gp_Vec temp = normal;
421
            if (temp * vertexnormals[i] < 0.0)
422
                temp = -temp;
423
            vertexnormals[i] = temp;
424

425
        }
426
        catch (...) {
427
        }
428

429
        vertexnormals[i].Normalize();
430
    }
431
}
432

433
void Part::Tools::getPointNormals(const TopoDS_Face& theFace, Handle(Poly_Triangulation) aPolyTri, TColgp_Array1OfDir& theNormals)
434
{
435
#if OCC_VERSION_HEX < 0x070600
436
    const TColgp_Array1OfPnt& aNodes = aPolyTri->Nodes();
437

438
    if(aPolyTri->HasNormals())
439
    {
440
        // normals pre-computed in triangulation structure
441
        const TShort_Array1OfShortReal& aNormals = aPolyTri->Normals();
442
        const Standard_ShortReal*       aNormArr = &(aNormals.Value(aNormals.Lower()));
443

444
        for(Standard_Integer aNodeIter = aNodes.Lower(); aNodeIter <= aNodes.Upper(); ++aNodeIter)
445
        {
446
            const Standard_Integer anId = 3 * (aNodeIter - aNodes.Lower());
447
            const gp_Dir aNorm(aNormArr[anId + 0],
448
                               aNormArr[anId + 1],
449
                               aNormArr[anId + 2]);
450
            theNormals(aNodeIter) = aNorm;
451
        }
452

453
        if(theFace.Orientation() == TopAbs_REVERSED)
454
        {
455
            for(Standard_Integer aNodeIter = aNodes.Lower(); aNodeIter <= aNodes.Upper(); ++aNodeIter)
456
            {
457
                theNormals.ChangeValue(aNodeIter).Reverse();
458
            }
459
        }
460
    }
461
    else {
462
        // take in face the surface location
463
        Poly_Connect thePolyConnect(aPolyTri);
464
        const TopoDS_Face      aZeroFace = TopoDS::Face(theFace.Located(TopLoc_Location()));
465
        Handle(Geom_Surface)   aSurf     = BRep_Tool::Surface(aZeroFace);
466
        const Standard_Real    aTol      = Precision::Confusion();
467
        Handle(TShort_HArray1OfShortReal) aNormals = new TShort_HArray1OfShortReal(1, aPolyTri->NbNodes() * 3);
468
        const Poly_Array1OfTriangle& aTriangles = aPolyTri->Triangles();
469
        const TColgp_Array1OfPnt2d*  aNodesUV   = aPolyTri->HasUVNodes() && !aSurf.IsNull()
470
                ? &aPolyTri->UVNodes()
471
                : nullptr;
472
        Standard_Integer aTri[3];
473

474
        for(Standard_Integer aNodeIter = aNodes.Lower(); aNodeIter <= aNodes.Upper(); ++aNodeIter)
475
        {
476
            // try to retrieve normal from real surface first, when UV coordinates are available
477
            if (!aNodesUV || GeomLib::NormEstim(aSurf, aNodesUV->Value(aNodeIter), aTol, theNormals(aNodeIter)) > 1)
478
            {
479
                // compute flat normals
480
                gp_XYZ eqPlan(0.0, 0.0, 0.0);
481

482
                for(thePolyConnect.Initialize(aNodeIter); thePolyConnect.More(); thePolyConnect.Next())
483
                {
484
                    aTriangles(thePolyConnect.Value()).Get(aTri[0], aTri[1], aTri[2]);
485
                    const gp_XYZ v1(aNodes(aTri[1]).Coord() - aNodes(aTri[0]).Coord());
486
                    const gp_XYZ v2(aNodes(aTri[2]).Coord() - aNodes(aTri[1]).Coord());
487
                    const gp_XYZ vv = v1 ^ v2;
488
                    const Standard_Real aMod = vv.Modulus();
489

490
                    if(aMod >= aTol)
491
                    {
492
                        eqPlan += vv / aMod;
493
                    }
494
                }
495

496
                const Standard_Real aModMax = eqPlan.Modulus();
497
                theNormals(aNodeIter) = (aModMax > aTol) ? gp_Dir(eqPlan) : gp::DZ();
498
            }
499

500
            const Standard_Integer anId = (aNodeIter - aNodes.Lower()) * 3;
501
            aNormals->SetValue(anId + 1, (Standard_ShortReal)theNormals(aNodeIter).X());
502
            aNormals->SetValue(anId + 2, (Standard_ShortReal)theNormals(aNodeIter).Y());
503
            aNormals->SetValue(anId + 3, (Standard_ShortReal)theNormals(aNodeIter).Z());
504
        }
505

506
        aPolyTri->SetNormals(aNormals);
507

508
        if(theFace.Orientation() == TopAbs_REVERSED)
509
        {
510
            for(Standard_Integer aNodeIter = aNodes.Lower(); aNodeIter <= aNodes.Upper(); ++aNodeIter)
511
            {
512
                theNormals.ChangeValue(aNodeIter).Reverse();
513
            }
514
        }
515
    }
516
#else
517
    Standard_Integer numNodes = aPolyTri->NbNodes();
518

519
    if(aPolyTri->HasNormals())
520
    {
521
        for(Standard_Integer aNodeIter = 1; aNodeIter <= numNodes; ++aNodeIter)
522
        {
523
            theNormals(aNodeIter) = aPolyTri->Normal(aNodeIter);
524
        }
525

526
        if(theFace.Orientation() == TopAbs_REVERSED)
527
        {
528
            for(Standard_Integer aNodeIter = 1; aNodeIter <= numNodes; ++aNodeIter)
529
            {
530
                theNormals.ChangeValue(aNodeIter).Reverse();
531
            }
532
        }
533
    }
534
    else {
535
        // take in face the surface location
536
        Poly_Connect thePolyConnect(aPolyTri);
537
        const TopoDS_Face      aZeroFace = TopoDS::Face(theFace.Located(TopLoc_Location()));
538
        Handle(Geom_Surface)   aSurf     = BRep_Tool::Surface(aZeroFace);
539
        const Standard_Real    aTol      = Precision::Confusion();
540
        Standard_Boolean hasNodesUV      = aPolyTri->HasUVNodes() && !aSurf.IsNull();
541
        Standard_Integer aTri[3];
542

543
        aPolyTri->AddNormals();
544
        for(Standard_Integer aNodeIter = 1; aNodeIter <= numNodes; ++aNodeIter)
545
        {
546
            // try to retrieve normal from real surface first, when UV coordinates are available
547
            if (!hasNodesUV || GeomLib::NormEstim(aSurf, aPolyTri->UVNode(aNodeIter), aTol, theNormals(aNodeIter)) > 1)
548
            {
549
                // compute flat normals
550
                gp_XYZ eqPlan(0.0, 0.0, 0.0);
551

552
                for(thePolyConnect.Initialize(aNodeIter); thePolyConnect.More(); thePolyConnect.Next())
553
                {
554
                    aPolyTri->Triangle(thePolyConnect.Value()).Get(aTri[0], aTri[1], aTri[2]);
555
                    const gp_XYZ v1(aPolyTri->Node(aTri[1]).Coord() - aPolyTri->Node(aTri[0]).Coord());
556
                    const gp_XYZ v2(aPolyTri->Node(aTri[2]).Coord() - aPolyTri->Node(aTri[1]).Coord());
557
                    const gp_XYZ vv = v1 ^ v2;
558
                    const Standard_Real aMod = vv.Modulus();
559

560
                    if(aMod >= aTol)
561
                    {
562
                        eqPlan += vv / aMod;
563
                    }
564
                }
565

566
                const Standard_Real aModMax = eqPlan.Modulus();
567
                theNormals(aNodeIter) = (aModMax > aTol) ? gp_Dir(eqPlan) : gp::DZ();
568
            }
569

570
            aPolyTri->SetNormal(aNodeIter, theNormals(aNodeIter));
571
        }
572

573
        if(theFace.Orientation() == TopAbs_REVERSED)
574
        {
575
            for(Standard_Integer aNodeIter = 1; aNodeIter <= numNodes; ++aNodeIter)
576
            {
577
                theNormals.ChangeValue(aNodeIter).Reverse();
578
            }
579
        }
580
    }
581
#endif
582
}
583

584
void Part::Tools::getPointNormals(const TopoDS_Face& face, Handle(Poly_Triangulation) aPoly, std::vector<gp_Vec>& normals)
585
{
586
    TColgp_Array1OfDir dirs (1, aPoly->NbNodes());
587
    getPointNormals(face, aPoly, dirs);
588
    normals.reserve(aPoly->NbNodes());
589

590
    for (int i = dirs.Lower(); i <= dirs.Upper(); ++i) {
591
        normals.emplace_back(dirs(i).XYZ());
592
    }
593
}
594

595
void Part::Tools::applyTransformationOnNormals(const TopLoc_Location& loc, std::vector<gp_Vec>& normals)
596
{
597
    if (!loc.IsIdentity()) {
598
        gp_Trsf myTransf = loc.Transformation();
599

600
        for (auto& it : normals) {
601
            it.Transform(myTransf);
602
        }
603
    }
604
}
605

606
Handle (Poly_Triangulation) Part::Tools::triangulationOfFace(const TopoDS_Face& face)
607
{
608
    TopLoc_Location loc;
609
    Handle (Poly_Triangulation) mesh = BRep_Tool::Triangulation(face, loc);
610
    if (!mesh.IsNull())
611
        return mesh;
612

613
    // If no triangulation exists then the shape is probably infinite
614
    double u1{}, u2{}, v1{}, v2{};
615
    try {
616
        BRepAdaptor_Surface adapt(face);
617
        u1 = adapt.FirstUParameter();
618
        u2 = adapt.LastUParameter();
619
        v1 = adapt.FirstVParameter();
620
        v2 = adapt.LastVParameter();
621
    }
622
    catch (const Standard_Failure&) {
623
        return nullptr;
624
    }
625

626
    auto selectRange = [](double& p1, double& p2) {
627
        if (Precision::IsInfinite(p1) && Precision::IsInfinite(p2)) {
628
            p1 = -50.0;
629
            p2 =  50.0;
630
        }
631
        else if (Precision::IsInfinite(p1)) {
632
            p1 = p2 - 100.0;
633
        }
634
        else if (Precision::IsInfinite(p2)) {
635
            p2 = p1 + 100.0;
636
        }
637
    };
638

639
    // recreate a face with a clear boundary in case it's infinite
640
    selectRange(u1, u2);
641
    selectRange(v1, v2);
642

643
    Handle(Geom_Surface) surface = BRep_Tool::Surface(face);
644
    if ( surface.IsNull() ) {
645
        FC_THROWM(Base::CADKernelError, "Cannot create surface from face");
646
    }
647
    BRepBuilderAPI_MakeFace mkBuilder(surface, u1, u2, v1, v2, Precision::Confusion() );
648
    TopoDS_Shape shape = mkBuilder.Shape();
649
    shape.Location(loc);
650

651
    BRepMesh_IncrementalMesh(shape, 0.005, false, 0.1, true);
652
    return BRep_Tool::Triangulation(TopoDS::Face(shape), loc);
653
}
654

655
Handle(Poly_Polygon3D) Part::Tools::polygonOfEdge(const TopoDS_Edge& edge, TopLoc_Location& loc)
656
{
657
    BRepAdaptor_Curve adapt(edge);
658
    double u = adapt.FirstParameter();
659
    double v = adapt.LastParameter();
660
    Handle(Poly_Polygon3D) aPoly = BRep_Tool::Polygon3D(edge, loc);
661
    if (!aPoly.IsNull() && !Precision::IsInfinite(u) && !Precision::IsInfinite(v))
662
        return aPoly;
663

664
    // recreate an edge with a clear range
665
    u = std::max(-50.0, u);
666
    v = std::min( 50.0, v);
667

668
    double uv;
669
    Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, uv, uv);
670

671
    BRepBuilderAPI_MakeEdge mkBuilder(curve, u, v);
672
    TopoDS_Shape shape = mkBuilder.Shape();
673
    // why do we have to set the inverted location here?
674
    TopLoc_Location inv = loc.Inverted();
675
    shape.Location(inv);
676

677
    BRepMesh_IncrementalMesh(shape, 0.005, false, 0.1, true);
678
    TopLoc_Location tmp;
679
    return BRep_Tool::Polygon3D(TopoDS::Edge(shape), tmp);
680
}
681

682
// helper function to use in getNormal, here we pass the local properties
683
// of the surface given by the #LProp_SLProps objects
684
template <typename T>
685
void getNormalBySLProp(T& prop, double u, double v, Standard_Real lastU, Standard_Real lastV,
686
                     const Standard_Real tol, gp_Dir& dir, Standard_Boolean& done)
687
{
688
    if (prop.D1U().Magnitude() > tol &&
689
        prop.D1V().Magnitude() > tol &&
690
        prop.IsNormalDefined()) {
691
        dir = prop.Normal();
692
        done = Standard_True;
693
    }
694
    // use an alternative method in case of a null normal
695
    else {
696
        CSLib_NormalStatus stat;
697
        CSLib::Normal(prop.D1U(), prop.D1V(), prop.D2U(), prop.D2V(), prop.DUV(),
698
            tol, done, stat, dir);
699
        // at the right boundary, the normal is flipped with respect to the
700
        // normal on surrounding points.
701
        if (stat == CSLib_D1NuIsNull) {
702
            if (Abs(lastV - v) < tol)
703
                dir.Reverse();
704
        }
705
        else if (stat == CSLib_D1NvIsNull || stat == CSLib_D1NuIsParallelD1Nv) {
706
            if (Abs(lastU - u) < tol)
707
                dir.Reverse();
708
        }
709
    }
710
}
711

712
void Part::Tools::getNormal(const Handle(Geom_Surface)& surf, double u, double v,
713
                            const Standard_Real tol, gp_Dir& dir, Standard_Boolean& done)
714
{
715
    GeomLProp_SLProps prop(surf, u, v, 1, tol);
716
    Standard_Real u1,u2,v1,v2;
717
    surf->Bounds(u1,u2,v1,v2);
718

719
    getNormalBySLProp<GeomLProp_SLProps>(prop, u, v, u2, v2, tol, dir, done);
720
}
721

722
void Part::Tools::getNormal(const TopoDS_Face& face, double u, double v,
723
                            const Standard_Real tol, gp_Dir& dir, Standard_Boolean& done)
724
{
725
    BRepAdaptor_Surface adapt(face);
726
    BRepLProp_SLProps prop(adapt, u, v, 1, tol);
727
    Standard_Real u2 = adapt.LastUParameter();
728
    Standard_Real v2 = adapt.LastVParameter();
729

730
    getNormalBySLProp<BRepLProp_SLProps>(prop, u, v, u2, v2, tol, dir, done);
731

732
    if (face.Orientation() == TopAbs_REVERSED)
733
        dir.Reverse();
734
}
735

736
TopLoc_Location Part::Tools::fromPlacement(const Base::Placement& plm)
737
{
738
    Base::Rotation r = plm.getRotation();
739
    double q1, q2, q3, q4;
740
    r.getValue(q1, q2, q3, q4);
741
    Base::Vector3d t = plm.getPosition();
742

743
    gp_Trsf trf;
744
    trf.SetTransformation(gp_Quaternion(q1, q2, q3, q4), gp_Vec(t.x, t.y, t.z));
745
    return {trf};
746
}
747