From 2fa69852709a645ac2ca2c61bbcef270f0506aca Mon Sep 17 00:00:00 2001 From: vsk Date: Fri, 18 Oct 2024 14:16:26 +0200 Subject: [PATCH 01/24] Modifications triggered by RevEng --- .../include/GoTools/creators/ApproxSurf.h | 30 +- .../include/GoTools/geometry/BoundedCurve.h | 3 +- .../include/GoTools/geometry/BoundedUtils.h | 39 +- .../include/GoTools/geometry/Circle.h | 3 +- .../include/GoTools/geometry/CurveOnSurface.h | 5 +- .../GoTools/geometry/ElementarySurface.h | 8 +- .../include/GoTools/geometry/Ellipse.h | 3 +- .../GoTools/geometry/GoIntersections.h | 6 + .../include/GoTools/geometry/Hyperbola.h | 5 +- gotools-core/include/GoTools/geometry/Line.h | 5 +- .../include/GoTools/geometry/Parabola.h | 5 +- .../include/GoTools/geometry/ParamCurve.h | 9 +- gotools-core/include/GoTools/geometry/Plane.h | 5 +- .../include/GoTools/geometry/RectDomain.h | 14 + .../include/GoTools/geometry/SplineCurve.h | 3 +- .../GoTools/geometry/SplineDebugUtils.h | 10 + gotools-core/include/GoTools/geometry/Torus.h | 51 +- gotools-core/include/GoTools/utils/Array.h | 7 + .../include/GoTools/utils/DirectionCone.h | 4 + gotools-core/src/creators/ApproxCurve.C | 39 +- gotools-core/src/creators/ApproxSurf.C | 165 ++++- gotools-core/src/creators/CreatorsUtils.C | 3 +- gotools-core/src/geometry/BoundedCurve.C | 4 +- gotools-core/src/geometry/BoundedUtils.C | 649 +++++++++++++++++- gotools-core/src/geometry/Circle.C | 50 +- gotools-core/src/geometry/Cone.C | 35 +- gotools-core/src/geometry/CurveOnSurface.C | 111 +-- gotools-core/src/geometry/Cylinder.C | 64 +- gotools-core/src/geometry/Ellipse.C | 4 +- gotools-core/src/geometry/FileUtils.C | 2 +- gotools-core/src/geometry/GSCappendCurve.C | 4 +- gotools-core/src/geometry/GeometryTools.C | 91 ++- gotools-core/src/geometry/Hyperbola.C | 3 +- gotools-core/src/geometry/Line.C | 19 +- gotools-core/src/geometry/Parabola.C | 3 +- gotools-core/src/geometry/Plane.C | 190 +++++ gotools-core/src/geometry/RectDomain.C | 61 +- gotools-core/src/geometry/SplineDebugUtils.C | 70 ++ gotools-core/src/geometry/SurfaceTools.C | 3 +- gotools-core/src/geometry/Torus.C | 20 +- gotools-core/src/geometry/closestPtCurves.C | 123 +++- gotools-core/src/geometry/intersectcurves.C | 54 ++ gotools-core/src/utils/BoundingBox.C | 13 +- gotools-core/src/utils/DirectionCone.C | 46 +- gotools-data | 2 +- sisl | 2 +- .../GoTools/trivariate/CurveOnVolume.h | 3 +- trivariate/src/CurveOnVolume.C | 4 +- 48 files changed, 1845 insertions(+), 207 deletions(-) diff --git a/gotools-core/include/GoTools/creators/ApproxSurf.h b/gotools-core/include/GoTools/creators/ApproxSurf.h index ad673ef06..ca6188f3f 100644 --- a/gotools-core/include/GoTools/creators/ApproxSurf.h +++ b/gotools-core/include/GoTools/creators/ApproxSurf.h @@ -60,7 +60,10 @@ #include "GoTools/geometry/SplineCurve.h" #include - +enum + { + ACCURACY_MAXDIST, ACCURACY_AVDIST_FRAC + }; class SmoothSurf; @@ -203,6 +206,11 @@ class ApproxSurf /// Approximate C1 continuity with a given impartance 0 <= fac < 1 void setC1Approx(double fac1, double fac2); + void setMBA(bool mba) + { + mba_ = mba; + } + /// Fetch the approximating surface /// When everything else is set, this function can be used to run the @@ -227,6 +235,11 @@ class ApproxSurf /// match against the current surface. int reParam(); + std::vector getParvals() + { + return parvals_; + } + /// Check whether or not the spline space will be refined during /// the approximation iterations bool getDoRefine() @@ -241,6 +254,13 @@ class ApproxSurf refine_ = refine; } + /// Set accuracy criterion + void setAccuracyCrit(int criterion, double frac) + { + acc_criter_ = (criterion == 1) ? ACCURACY_AVDIST_FRAC : ACCURACY_MAXDIST; + acc_frac_ = frac; + } + protected: /// Default constructor @@ -264,6 +284,7 @@ class ApproxSurf bool close_belt_; bool repar_; bool refine_; + bool mba_; int dim_; std::vector points_; @@ -275,6 +296,8 @@ class ApproxSurf int constdir_; bool orig_; double c1fac1_, c1fac2_; + int acc_criter_; + double acc_frac_; /// Generate an initial curve representing the spline space int makeInitSurf(std::vector > &crvs, @@ -289,6 +312,9 @@ class ApproxSurf /// Generate a smoothing surface int makeSmoothSurf(); + /// Approximate using mba method + void mbaApprox(); + /// Check the accuracy of the current surface int checkAccuracy(std::vector& acc_outside_u, std::vector& nmb_outside_u, @@ -307,6 +333,8 @@ class ApproxSurf /// Define free and fixed coefficients void coefKnownFromPoints(); void setCoefKnown(); + + bool approxOK(); }; } // namespace Go diff --git a/gotools-core/include/GoTools/geometry/BoundedCurve.h b/gotools-core/include/GoTools/geometry/BoundedCurve.h index 388cdc29d..6799c0d34 100644 --- a/gotools-core/include/GoTools/geometry/BoundedCurve.h +++ b/gotools-core/include/GoTools/geometry/BoundedCurve.h @@ -134,7 +134,8 @@ class GO_API BoundedCurve : public ParamCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, diff --git a/gotools-core/include/GoTools/geometry/BoundedUtils.h b/gotools-core/include/GoTools/geometry/BoundedUtils.h index 5ab8011e4..868c57566 100644 --- a/gotools-core/include/GoTools/geometry/BoundedUtils.h +++ b/gotools-core/include/GoTools/geometry/BoundedUtils.h @@ -286,6 +286,11 @@ namespace BoundedUtils { std::vector< shared_ptr< CurveOnSurface > >& bnd_cvs, double eps); + std::vector > + intersectWithElementarySurface(shared_ptr& surf, + shared_ptr& elem, + double geom_tol); + /// Find the intersection curve(s) between a parametric surface and a given plane. /// The plane is defined by its normal and a point on the plane. /// \param surf the surface to intersect with the plane @@ -322,7 +327,21 @@ namespace BoundedUtils { intersectWithCylinder(shared_ptr& surf, Point pnt, Point vec, double radius, double geom_tol); - /// Find the intersction curve(s) between two spline surfaces + std::vector > + intersectWithSphere(shared_ptr& surf, + Point pnt, double radius, double geom_tol); + + std::vector > + intersectWithCone(shared_ptr& surf, + Point pnt, Point axispt, Point sfpt, double geom_tol); + + std::vector > + intersectWithTorus(shared_ptr& surf, + Point pnt, Point normal, double rad1, + double rad2, double geom_tol); + + /// Find the intersction curve(s) between two parametric surfaces. The surfaces + /// are represented as spline surface if this was not the case initially. /// \param sf1 the first surface to intersect /// \param sf2 the second surface to intersect /// \retval int_segments1 a vector of CurveOnSurface s, representing the intersection @@ -337,7 +356,23 @@ namespace BoundedUtils { double epsge); - /// Translate a given BoundedSurface + /// Find the intersction curve(s) between two parametric surfaces where at least + /// one is an elementary surface. + /// \param sf1 the first surface to intersect + /// \param sf2 the second surface to intersect + /// \retval int_segments1 a vector of CurveOnSurface s, representing the intersection + /// curves as lying on 'sf1'. + /// \retval int_segments2 a vector of CurveOnSurface s, representing the intersection + /// curves as lying on 'sf2'. + /// \param epsge geometrical tolerance to be used for intersection computations + void getIntersectionCurveElem(shared_ptr& sf1, + shared_ptr& sf2, + std::vector >& int_segments1, + std::vector >& int_segments2, + double epsge); + + + /// Translate a given BoundedSurface /// \param trans_vec the vector specifying the translation to apply to the surface /// \param bd_sf the surface to translate /// \param deg_eps an epsilon value used when determining degenerate boundary loops diff --git a/gotools-core/include/GoTools/geometry/Circle.h b/gotools-core/include/GoTools/geometry/Circle.h index d61ac1f2b..d2a6706da 100644 --- a/gotools-core/include/GoTools/geometry/Circle.h +++ b/gotools-core/include/GoTools/geometry/Circle.h @@ -139,7 +139,8 @@ class Circle : public ElementaryCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, diff --git a/gotools-core/include/GoTools/geometry/CurveOnSurface.h b/gotools-core/include/GoTools/geometry/CurveOnSurface.h index 7f0998300..60a0aa637 100644 --- a/gotools-core/include/GoTools/geometry/CurveOnSurface.h +++ b/gotools-core/include/GoTools/geometry/CurveOnSurface.h @@ -248,7 +248,8 @@ class GO_API CurveOnSurface : public ParamCurve // inherited from ParamCurve. NB: does not check whether the resulting ParamCurve // stays inside parameter domain (or that the space curve stays on surface). virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); /// Set the underlying surface to the one pointed to by the argument /// \param surface the pointer to the surface we will set as underlying for this @@ -441,6 +442,8 @@ class GO_API CurveOnSurface : public ParamCurve /// Check if the curve is a constant parameter curve with regard to the /// parameter tol + /// pardir = 1: constant in 1. parameter direction + /// pardir = 2: constant in 2. parameter direction bool isConstantCurve(double tol, int& pardir, double& parval) const; /// Fetch recorded constant curve information diff --git a/gotools-core/include/GoTools/geometry/ElementarySurface.h b/gotools-core/include/GoTools/geometry/ElementarySurface.h index 8d4cd63c4..f8626440c 100644 --- a/gotools-core/include/GoTools/geometry/ElementarySurface.h +++ b/gotools-core/include/GoTools/geometry/ElementarySurface.h @@ -134,10 +134,11 @@ class ElementarySurface : public ParamSurface virtual void setParameterBounds(double from_upar, double from_vpar, double to_upar, double to_vpar) = 0; + /// Fetch parameter bounds. NB! Not oriented virtual RectDomain getParameterBounds() const = 0; - virtual void turnOrientation(); + virtual void turnOrientation(); virtual void reverseParameterDirection(bool direction_is_u); virtual void swapParameterDirection(); @@ -166,6 +167,11 @@ class ElementarySurface : public ParamSurface return -1.0; } + virtual double radius2(double u, double v) const + { + return -1.0; + } + virtual Point location() const { Point dummy; diff --git a/gotools-core/include/GoTools/geometry/Ellipse.h b/gotools-core/include/GoTools/geometry/Ellipse.h index 386c7cbcd..c8e0da287 100644 --- a/gotools-core/include/GoTools/geometry/Ellipse.h +++ b/gotools-core/include/GoTools/geometry/Ellipse.h @@ -126,7 +126,8 @@ class Ellipse : public ElementaryCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, diff --git a/gotools-core/include/GoTools/geometry/GoIntersections.h b/gotools-core/include/GoTools/geometry/GoIntersections.h index 4dcb26043..e6003e3f7 100644 --- a/gotools-core/include/GoTools/geometry/GoIntersections.h +++ b/gotools-core/include/GoTools/geometry/GoIntersections.h @@ -115,6 +115,12 @@ void intersectParamCurves(ParamCurve* cv1, ParamCurve* cv2, double epsge, double epsge, std::vector& intpar, std::vector >& int_cvs); +/// Intersect a curve with a cylinder + void intersectCurveCylinder(ParamCurve* cv, const Point& pos, + const Point&axis, double radius, + double epsge, std::vector& intpar, + std::vector >& int_cvs); + ///Intersect two spline curves. Collect intersection parameters. /// \param cv1 pointer to the first spline curve /// \param cv2 pointer to the second spline curve diff --git a/gotools-core/include/GoTools/geometry/Hyperbola.h b/gotools-core/include/GoTools/geometry/Hyperbola.h index 6f63e373d..d24ed3917 100644 --- a/gotools-core/include/GoTools/geometry/Hyperbola.h +++ b/gotools-core/include/GoTools/geometry/Hyperbola.h @@ -126,7 +126,8 @@ class Hyperbola : public ElementaryCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, @@ -151,7 +152,7 @@ class Hyperbola : public ElementaryCurve /// Query if parametrization is bounded. Both upper and lower /// parameter bounds must be finite for this to be true. /// \return \a true if bounded, \a false otherwise - bool isBounded() const; + virtual bool isBounded() const; /// If the curve is 2 dimensional, x and y coordinates will be swapped. /// Used when curve is a parameter curve. diff --git a/gotools-core/include/GoTools/geometry/Line.h b/gotools-core/include/GoTools/geometry/Line.h index 120389234..3c9b86c0a 100644 --- a/gotools-core/include/GoTools/geometry/Line.h +++ b/gotools-core/include/GoTools/geometry/Line.h @@ -132,7 +132,8 @@ class Line : public ElementaryCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, @@ -165,7 +166,7 @@ class Line : public ElementaryCurve /// Query if parametrization is bounded. Both upper and lower /// parameter bounds must be finite for this to be true. /// \return \a true if bounded, \a false otherwise - bool isBounded() const; + virtual bool isBounded() const; /// Confirm that the curve is linear virtual bool isLinear(Point& dir, double tol); diff --git a/gotools-core/include/GoTools/geometry/Parabola.h b/gotools-core/include/GoTools/geometry/Parabola.h index b310d885f..a9b7f3eca 100644 --- a/gotools-core/include/GoTools/geometry/Parabola.h +++ b/gotools-core/include/GoTools/geometry/Parabola.h @@ -126,7 +126,8 @@ class Parabola : public ElementaryCurve virtual void appendCurve(ParamCurve* cv, bool reparam=true); virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); virtual void closestPoint(const Point& pt, double tmin, @@ -148,7 +149,7 @@ class Parabola : public ElementaryCurve /// Query if parametrization is bounded. Both upper and lower /// parameter bounds must be finite for this to be true. /// \return \a true if bounded, \a false otherwise - bool isBounded() const; + virtual bool isBounded() const; // Translate the curve along a given vector virtual void translateCurve(const Point& dir); diff --git a/gotools-core/include/GoTools/geometry/ParamCurve.h b/gotools-core/include/GoTools/geometry/ParamCurve.h index 41608fcdc..db5fd109e 100644 --- a/gotools-core/include/GoTools/geometry/ParamCurve.h +++ b/gotools-core/include/GoTools/geometry/ParamCurve.h @@ -237,8 +237,10 @@ class GO_API ParamCurve : public GeomObject /// \param dist a measure of the local distorsion around the transition in order /// to achieve the specified continuity. /// \param reparam specify whether or not there should be reparametrization + /// \param tol relevant only for CurveOnSurface virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true) = 0; + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4) = 0; /// Estimate the length of the curve, by sampling it at a certain number of points /// and calculating the linear approximation to the curve through these points. @@ -331,6 +333,11 @@ class GO_API ParamCurve : public GeomObject // bounded curves and some elementary curves } + virtual bool isBounded() const + { + return true; // Is overridden when an unbounded curve is possible + } + /// Check if the curve is linear virtual bool isLinear(Point& dir, double tol); diff --git a/gotools-core/include/GoTools/geometry/Plane.h b/gotools-core/include/GoTools/geometry/Plane.h index 097fc4593..a8282f555 100644 --- a/gotools-core/include/GoTools/geometry/Plane.h +++ b/gotools-core/include/GoTools/geometry/Plane.h @@ -166,6 +166,9 @@ class Plane : public ElementarySurface /// Check for paralell and anti paralell partial derivatives in surface corners virtual void getDegenerateCorners(std::vector& deg_corners, double tol) const; + virtual shared_ptr + getElementaryParamCurve(ElementaryCurve* space_crv, double tol, + const Point* start_par_pt = NULL, const Point* end_par_pt = NULL) const; // --- Functions specific to Plane --- /// Point in plane @@ -228,7 +231,7 @@ class Plane : public ElementarySurface /// \return \a true if bounded, \a false otherwise virtual bool isBounded() const; - /// Check if the plane is closed. Virtual function - always false. + /// Check if the plane is closed. Virtual function - always false. virtual bool isClosed(bool& closed_dir_u, bool& closed_dir_v) const; /// Return the result from intersecting the unbounded plane with a diff --git a/gotools-core/include/GoTools/geometry/RectDomain.h b/gotools-core/include/GoTools/geometry/RectDomain.h index 588ba8a81..2d6b65a9e 100644 --- a/gotools-core/include/GoTools/geometry/RectDomain.h +++ b/gotools-core/include/GoTools/geometry/RectDomain.h @@ -82,6 +82,16 @@ class GO_API RectDomain : public Domain virtual bool isOnBoundary(const Array& point, double tolerance) const; + /// bd = -1 : Not a boundary point + /// bd = 0: umin + /// bd = 1: umax + /// bd = 2: vmin + /// bd = 3: vmax + /// bd2 = -1 : Not a corner point + /// otherwise: as bd + bool isOnBoundary(const Array& point, + double tolerance, int& bd, int& bd2) const; + /// Check if a given parameter pair lies on a corner in the domain within /// the given tolerance bool isOnCorner(const Array& point, @@ -144,6 +154,10 @@ class GO_API RectDomain : public Domain /// included bool overlap(const RectDomain& rd, double tol); + /// Check if two domains overlap, boundary overlap within tolerance + /// dependent on parameter direction included + bool overlap(const RectDomain& rd, double tol1, double tol2); + /// Get the 'lower left' corner of this RectDomain. /// \return a 2D array containing the 'lower left' corner of this RectDomain Array lowerLeft() const { return ll_; } diff --git a/gotools-core/include/GoTools/geometry/SplineCurve.h b/gotools-core/include/GoTools/geometry/SplineCurve.h index aeca2e476..769cbabb7 100644 --- a/gotools-core/include/GoTools/geometry/SplineCurve.h +++ b/gotools-core/include/GoTools/geometry/SplineCurve.h @@ -294,7 +294,8 @@ class GO_API SplineCurve : public ParamCurve virtual void appendCurve(ParamCurve* other_curve, int continuity, double& dist, - bool repar=true); + bool repar=true, + double tol = 1.0e-4); /// Inherited from ParamCurve. /// Short hand function to call \ref appendCurve with C^1 diff --git a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h index a2f7c0f4f..21e9945ba 100644 --- a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h +++ b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h @@ -44,6 +44,7 @@ #include "GoTools/geometry/SplineSurface.h" #include "GoTools/geometry/SplineCurve.h" #include "GoTools/geometry/Line.h" +#include "GoTools/geometry/Circle.h" #include #include "GoTools/utils/config.h" #include "GoTools/geometry/BoundedSurface.h" @@ -65,9 +66,18 @@ namespace SplineDebugUtils void GO_API writeSpaceParamCurve(const SplineCurve& pcurve, std::ostream& os, double z = 0.0); + void GO_API writeSpaceParamCurve(shared_ptr pcurve, + std::ostream& os, double z = 0.0); + void GO_API writeSpaceParamCurve(const Line& pline, std::ostream& os, double z = 0.0); + void GO_API writeSpaceParamCurve(const Circle& pcirc, + std::ostream& os, double z = 0.0); + + void GO_API writeSpace1DCurve(const SplineCurve& pcurve, + std::ostream& os, double z = 0.0); + /// Write the parameter curve (if existing) and the space curve (if /// existing) to the output stream. Both curves are written as 3D curves /// extending 2D curves with the given z-value. diff --git a/gotools-core/include/GoTools/geometry/Torus.h b/gotools-core/include/GoTools/geometry/Torus.h index 6a81bbdc7..0eba59b8c 100644 --- a/gotools-core/include/GoTools/geometry/Torus.h +++ b/gotools-core/include/GoTools/geometry/Torus.h @@ -186,7 +186,37 @@ class Torus : public ElementarySurface /// Check for paralell and anti paralell partial derivatives in surface corners virtual void getDegenerateCorners(std::vector& deg_corners, double tol) const; - // --- Functions specific to Torus --- + virtual double radius(double u, double v) const + { + return major_radius_; + } + + virtual double radius2(double u, double v) const + { + return minor_radius_; + } + + virtual Point location() const + { + return location_; + } + + virtual Point direction() const + { + return z_axis_; + } + + virtual Point direction2() const + { + return x_axis_; + } + + virtual void enlarge(double len1, double len2, double len3, double len4); + + virtual void translate(const Point& vec); + + + // --- Functions specific to Torus --- /// Fetch major radius double getMajorRadius() const @@ -278,25 +308,6 @@ class Torus : public ElementarySurface /// returned. shared_ptr getMinorCircle(double upar) const; - virtual Point location() const - { - return location_; - } - - virtual Point direction() const - { - return z_axis_; - } - - virtual Point direction2() const - { - return x_axis_; - } - - virtual void enlarge(double len1, double len2, double len3, double len4); - - virtual void translate(const Point& vec); - protected: double major_radius_; diff --git a/gotools-core/include/GoTools/utils/Array.h b/gotools-core/include/GoTools/utils/Array.h index 38051d536..668b92c6f 100644 --- a/gotools-core/include/GoTools/utils/Array.h +++ b/gotools-core/include/GoTools/utils/Array.h @@ -281,6 +281,13 @@ namespace Go (*this) /= length(); } + void normalize_checked() + { + double len = length(); + if (len > 1.0e-12) + (*this) /= len; + } + /// The sum of two vectors. Array operator + (const Array &v) const { diff --git a/gotools-core/include/GoTools/utils/DirectionCone.h b/gotools-core/include/GoTools/utils/DirectionCone.h index 9f3c4f1e7..bb743bc72 100644 --- a/gotools-core/include/GoTools/utils/DirectionCone.h +++ b/gotools-core/include/GoTools/utils/DirectionCone.h @@ -125,6 +125,10 @@ class GO_API DirectionCone { /// it covers the direction represented by 'pt'. void addUnionWith(const Point& pt); + /// Test run. Which angle would an added vector result in. Angles larger +/// than M_PI is reported as M_PI. +double unionAngle(const Point& pt) const; + /// If necessary, increase the angle of the DirectionCone so that /// it covers all directions covered by 'cone'. void addUnionWith(const DirectionCone& cone); diff --git a/gotools-core/src/creators/ApproxCurve.C b/gotools-core/src/creators/ApproxCurve.C index ca183338c..da19b00a2 100644 --- a/gotools-core/src/creators/ApproxCurve.C +++ b/gotools-core/src/creators/ApproxCurve.C @@ -423,7 +423,7 @@ void ApproxCurve::checkAccuracy(std::vector& newknots, int uniform) curr_crv_->writeStandardHeader(of); curr_crv_->write(of); #endif - bool reparam = true; + bool reparam = (dim_ == 1) ? false : true; // double par_tol = 0.000000000001; maxdist_ = -10000.0; @@ -445,20 +445,29 @@ void ApproxCurve::checkAccuracy(std::vector& newknots, int uniform) double *pt, *param; int ki; for (ki=0, pt=&points_[0], param=&parvals_[0]; kistartparam(); - tb = curr_crv_->endparam(); - pos.setValue(pt); - try { - curr_crv_->closestPoint(pos, ta, tb, - par, clpos, dist, param); - } catch (...) { - MESSAGE("Failed finding closest point."); - par = *param; // We're using our input value. - // This should at least be an upper boundary of closest dist. + ki++, pt+=dim_, param++) + { + if (reparam || dim_ > 1) + { + // Compute closest point + ta = curr_crv_->startparam(); + tb = curr_crv_->endparam(); + pos.setValue(pt); + try { + curr_crv_->closestPoint(pos, ta, tb, + par, clpos, dist, param); + } catch (...) { + MESSAGE("Failed finding closest point."); + par = *param; // We're using our input value. + // This should at least be an upper boundary of closest dist. + dist = (pos - curr_crv_->ParamCurve::point(par)).length(); + } + } + else + { + par = *param; dist = (pos - curr_crv_->ParamCurve::point(par)).length(); - } + } left = curr_crv_->basis().knotInterval(par); if (reparam) parvals_[ki] = par; @@ -548,7 +557,7 @@ int ApproxCurve::doApprox(int max_iter) { coef[ki] = points_[ki]; coef[(in-1)*dim_+ki] = - points_[(nmbpoints-1)*dim_+ki]; + points_[(nmbpoints-1)*dim_+ki]; } // Approximate diff --git a/gotools-core/src/creators/ApproxSurf.C b/gotools-core/src/creators/ApproxSurf.C index 1a0113f88..cf38194a5 100644 --- a/gotools-core/src/creators/ApproxSurf.C +++ b/gotools-core/src/creators/ApproxSurf.C @@ -42,6 +42,7 @@ #include "GoTools/creators/SmoothSurf.h" #include "GoTools/geometry/CurveLoop.h" #include "GoTools/creators/CoonsPatchGen.h" +#include "GoTools/utils/Array.h" #include #include #include @@ -88,8 +89,10 @@ ApproxSurf::ApproxSurf() orig_ = false; repar_ = true; refine_ = true; + mba_ = false; c1fac1_ = 0.0; c1fac2_ = 0.0; + acc_criter_ = ACCURACY_MAXDIST; } //*************************************************************************** @@ -128,8 +131,10 @@ ApproxSurf::ApproxSurf(std::vector >& crvs, orig_ = false; repar_ = repar; refine_ = true; + mba_ = false; c1fac1_ = 0.0; c1fac2_ = 0.0; + acc_criter_ = ACCURACY_MAXDIST; makeInitSurf(crvs, domain); @@ -173,8 +178,10 @@ ApproxSurf::ApproxSurf(shared_ptr& srf, orig_ = approx_orig; repar_ = repar; refine_ = true; + mba_ = false; c1fac1_ = 0.0; c1fac2_ = 0.0; + acc_criter_ = ACCURACY_MAXDIST; points_ = points; parvals_ = parvals; @@ -626,7 +633,10 @@ int ApproxSurf::doApprox(int max_iter, int keep_init) prev_srf_ = shared_ptr(curr_srf_->clone()); // Store last version prevdist_ = maxdist_; prevav_ = avdist_; - stat = makeSmoothSurf(); + if (mba_) + mbaApprox(); + else + stat = makeSmoothSurf(); if (stat < 0) return stat; @@ -646,10 +656,12 @@ int ApproxSurf::doApprox(int max_iter, int keep_init) cout << avdist_ << " # out " << outsideeps_ << endl; #endif - if (maxdist_ < aepsge_) + bool OK = approxOK(); + if (OK) break; - if (maxdist_ > prevdist_) + //if (maxdist_ > prevdist_) + if (avdist_ > prevav_) { curr_srf_ = prev_srf_; maxdist_ = prevdist_; @@ -687,6 +699,140 @@ int ApproxSurf::doApprox(int max_iter, int keep_init) return 0; } +//*************************************************************************** +void add_contribution(int dim, vector >& target, + int ix, double nom[], double denom) +{ + for (int ki=0; kiorder_u(); + int order_v = curr_srf_->order_v(); + int ncoef_u = curr_srf_->numCoefs_u(); + int ncoef_v = curr_srf_->numCoefs_v(); + + for (int iter=0; iter<2; ++iter) + { + // Vector to accumulate numerator and denominator to compute final coefficient value + // for each b-spline + Array nom_denom0(0.0); + vector > nom_denom(ncoef_u*ncoef_v, nom_denom0); + + // Temporary vector to store weights associated with a given data point + vector tmpvec(dim_); + + int ka, kb, kc, ix1, ix2; + const BsplineBasis bspline_u = curr_srf_->basis_u(); + const BsplineBasis bspline_v = curr_srf_->basis_v(); + vector tmp_wgts(order_u*order_v, 0.0); + for (int ki=0; ki basis_u = bspline_u.computeBasisValues(parvals_[2*ki]); + vector basis_v = bspline_v.computeBasisValues(parvals_[2*ki+1]); + const int uleft = bspline_u.lastKnotInterval(); + const int vleft = bspline_v.lastKnotInterval(); + + // Compute surface position + Point pos(dim_); + pos.setValue(0.0); + std::vector::iterator it2 = curr_srf_->coefs_begin(); + std::vector::iterator it1; + for (kb=0, it2+=(ncoef_u*(vleft-order_v+1)+uleft-order_u+1)*dim_; kbParamSurface::point(parvals_[2*ki], parvals_[2*ki+1]); + + Point pnt(points_.begin()+ki*dim_, points_.begin()+(ki+1)*dim_); + Point distvec = pnt - pos; + + // Computing weights for the current data point + it2 = curr_srf_->coefs_begin(); + std::vector::iterator itknown2 = coef_known_.begin(); + std::vector::iterator itknown1; + double total_squared_inv = 0; + for (kb=0, it2+=(ncoef_u*(vleft-order_v+1)+uleft-order_u+1)*dim_, + itknown2+=(ncoef_u*(vleft-order_v+1)+uleft-order_u+1); kbcoefs_begin(); + int ix1, ix2; + for (kb=0, it2+=(ncoef_u*(vleft-order_v+1)+uleft-order_u+1)*dim_, + ix2=ncoef_u*(vleft-order_v+1)+uleft-order_u+1; kb::iterator it1 = curr_srf_->coefs_begin(); + for (kb=0, ix1=0; kb curr = nom_denom[ix1]; + for (kc=0; kcreplaceCoefficient(ix1, coef+coef2); + } + int stop_break = 1; + } +} + + //*************************************************************************** shared_ptr @@ -936,3 +1082,16 @@ void ApproxSurf::setC1Approx(double fac1, double fac2) c1fac1_ = fac1; c1fac2_ = fac2; } + + +bool ApproxSurf::approxOK() +{ + if (acc_criter_ == ACCURACY_AVDIST_FRAC) + { + size_t nmb_pts = points_.size()/dim_; + bool OK = (avdist_ <= aepsge_ && (double)outsideeps_/(double)nmb_pts > acc_frac_); + return OK; + } + else + return (maxdist_ <= aepsge_); +} diff --git a/gotools-core/src/creators/CreatorsUtils.C b/gotools-core/src/creators/CreatorsUtils.C index 8c8820ebe..407a1004d 100644 --- a/gotools-core/src/creators/CreatorsUtils.C +++ b/gotools-core/src/creators/CreatorsUtils.C @@ -197,7 +197,8 @@ createCrossTangent(const Go::CurveOnSurface& cv, smooth_cv.attach(basis_space_cv, &coef_known[0]); } else { //smooth_cv.attach(space_cv, seem, &coef_known[0]); - smooth_cv.attach(space_cv, &coef_known[0]); + shared_ptr space_cv2(space_cv->clone()); + smooth_cv.attach(space_cv2, &coef_known[0]); } double smoothweight = 0.000000001; diff --git a/gotools-core/src/geometry/BoundedCurve.C b/gotools-core/src/geometry/BoundedCurve.C index 3c473478c..afe34c55c 100644 --- a/gotools-core/src/geometry/BoundedCurve.C +++ b/gotools-core/src/geometry/BoundedCurve.C @@ -309,8 +309,8 @@ void BoundedCurve::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== -void BoundedCurve::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) +void BoundedCurve::appendCurve(ParamCurve* cv, int continuity, double& dist, + bool reparam, double tol) //=========================================================================== { MESSAGE("Not implemented!"); diff --git a/gotools-core/src/geometry/BoundedUtils.C b/gotools-core/src/geometry/BoundedUtils.C index 4baede2fb..33c6b893c 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -37,7 +37,7 @@ * written agreement between you and SINTEF ICT. */ -//#define DEBUG1 +#define DEBUG1 //#define SBR_DBG @@ -61,6 +61,12 @@ #include "GoTools/geometry/SurfaceTools.h" #include "GoTools/geometry/ClosestPoint.h" #include "GoTools/geometry/SurfaceOfLinearExtrusion.h" +#include "GoTools/geometry/Plane.h" +#include "GoTools/geometry/Cylinder.h" +#include "GoTools/geometry/Cone.h" +#include "GoTools/geometry/Sphere.h" +#include "GoTools/geometry/Torus.h" +#include "GoTools/geometry/ClosestPoint.h" #include "GoTools/creators/CoonsPatchGen.h" @@ -719,12 +725,28 @@ BoundedUtils::getSurfaceIntersections(const shared_ptr& surf1, under_sf2->writeStandardHeader(debug0); under_sf2->write(debug0); #endif - try { - getIntersectionCurve(under_sf1, under_sf2, int_segments1, int_segments2, - epsge); - } catch (...) { - THROW("Failed intersecting spline surface with plane."); - } + shared_ptr elem1 = + dynamic_pointer_cast(under_sf1); + shared_ptr elem2 = + dynamic_pointer_cast(under_sf2); + if (elem1.get() || elem2.get()) + { + try { + getIntersectionCurveElem(under_sf1, under_sf2, int_segments1, int_segments2, + epsge); + } catch (...) { + THROW("Failed intersecting surfaces where at least one is an elementary surface."); + } + } + else + { + try { + getIntersectionCurve(under_sf1, under_sf2, int_segments1, int_segments2, + epsge); + } catch (...) { + THROW("Failed intersecting spline surfaces."); + } + } // Ensure correct surface pointer for (size_t kj=0; kj > >& } +//=========================================================================== +vector > +BoundedUtils::intersectWithElementarySurface(shared_ptr& surf, + shared_ptr& elem, + double geom_tol) +//=========================================================================== +{ + vector > int_seg; + if (elem->instanceType() == Class_Plane) + { + shared_ptr plane = dynamic_pointer_cast(elem); + Point pnt = plane->getPoint(); + Point normal = plane->getNormal(); + int_seg = intersectWithPlane(surf, pnt, normal, geom_tol); + } + else if (elem->instanceType() == Class_Cylinder) + { + shared_ptr cyl = dynamic_pointer_cast(elem); + Point loc = cyl->getLocation(); + Point dir = cyl->getAxis(); + double rad = cyl->getRadius(); + int_seg = intersectWithCylinder(surf, loc, dir, rad, geom_tol); + } + else if (elem->instanceType() == Class_Cone) + { + shared_ptr cone = dynamic_pointer_cast(elem); + Point loc = cone->getLocation(); + Point dir = cone->getAxis(); + Point dir2 = cone->direction2(); + double rad = cone->getRadius(); + double alpha = cone->getConeAngle(); + if (rad < 0.01) + { + RectDomain dom = cone->getParameterBounds(); + double upar = 0.5*(dom.umin()+dom.umax()); + double vpar = 0.5*(dom.vmin()+dom.vmax()); + rad = cone->radius(upar, vpar); + } + double dd = rad/tan(alpha); + Point top = loc - dd*dir; + Point axispt = loc + rad*dir; + Point sfpt = loc + rad*dir2; + + int_seg = intersectWithCone(surf, top, loc, sfpt, geom_tol); + } + else if (elem->instanceType() == Class_Sphere) + { + shared_ptr sph = dynamic_pointer_cast(elem); + Point centre = sph->getLocation(); + double rad = sph->getRadius(); + int_seg = intersectWithSphere(surf, centre, rad, geom_tol); + } + else if (elem->instanceType() == Class_Torus) + { + shared_ptr tor = dynamic_pointer_cast(elem); + Point centre = tor->getLocation(); + Point dir = tor->direction(); + double rad1 = tor->getMajorRadius(); + double rad2 = tor->getMinorRadius(); + int_seg = intersectWithTorus(surf, centre, dir, rad1, rad2, geom_tol); + } + // else + // { + // vector > dummy; + // return dummy; + // } + + if (int_seg.size() > 0 && elem->isBounded()) + { + // Trim the intersection curves with the boundaries of the + // elementary surface + const Point sf_epspar = SurfaceTools::getParEpsilon(*elem, geom_tol); + const double epspar = std::min(sf_epspar[0], sf_epspar[1]); + double angtol = 0.01; + CurveLoop cvloop = elem->outerBoundaryLoop(); + vector > int_seg2; + vector > bdcvs = cvloop.getCurves(); + for (size_t ki=0; ki bd_par; + for (size_t kj=0; kjstartparam(); + double end = int_seg[ki]->endparam(); + double par1, par2, dist; + Point ptc1, ptc2; + ClosestPoint::closestPtCurves(int_seg[ki].get(), bdcvs[kj].get(), + par1, par2, dist, ptc1, ptc2); + if (dist < geom_tol && par1-start > epspar && end-par1 > epspar) + bd_par.push_back(par1); + } + if (bd_par.size() > 0) + std::sort(bd_par.begin(), bd_par.end()); + if (bd_par.size() == 0 || bd_par[0] > int_seg[ki]->startparam()+epspar) + bd_par.insert(bd_par.begin(), int_seg[ki]->startparam()); + if (bd_par.size() == 0 || + bd_par[bd_par.size()-1] < int_seg[ki]->endparam()-epspar) + bd_par.push_back(int_seg[ki]->endparam()); + + for (size_t kj=bd_par.size()-1; kj>0; --kj) + if (bd_par[kj] - bd_par[kj-1] < epspar) + { + bd_par[kj-1] = 0.5*(bd_par[kj-1]+bd_par[kj]); + bd_par.erase(bd_par.begin()+kj); + } + + vector inside(bd_par.size()-1, true); + for (int ka=(int)bd_par.size()-2; ka>=0; --ka) + { + Point midp; + double par = 0.5*(bd_par[ka] + bd_par[ka+1]); + int_seg[ki]->point(midp, par); + double upar, vpar, dist; + Point close; + elem->closestPoint(midp, upar, vpar, close, dist, epspar); + if (dist > geom_tol) + { + Point norm; + elem->normal(norm, upar, vpar); + Point vec = close - midp; + double ang = vec.angle(norm); + ang = std::min(ang, M_PI-ang); + if (ang > angtol) + inside[ka] = false; + } + } + + size_t kr; + for (size_t kj=0; kjint_seg[ki]->startparam()+epspar || + bd_par[kr]endparam()-epspar) + { + shared_ptr sub(int_seg[ki]->subCurve(bd_par[kj], + bd_par[kr])); + int_seg2.push_back(sub); + } + else + int_seg2.push_back(int_seg[ki]); + } + } + int stop_break = 1; + } + return int_seg2; + } + else + return int_seg; +} + //=========================================================================== vector > BoundedUtils::intersectWithPlane(shared_ptr& surf, @@ -2363,6 +2537,467 @@ BoundedUtils::intersectWithCylinder(shared_ptr& surf, return curves; } +//=========================================================================== +vector > +BoundedUtils::intersectWithSphere(shared_ptr& surf, + Point pnt, double radius, + double geom_tol) +//=========================================================================== +{ + vector > curves; + + // Convert the surface to a SISLSurf in order to use SISL functions + // on it. The "false" argument dictates that the SISLSurf will only + // copy pointers to arrays, not the arrays themselves. + shared_ptr tmp_sf; + SplineSurface *splinesf = surf->getSplineSurface(); + if (!splinesf) + { + tmp_sf = shared_ptr(surf->asSplineSurface()); + splinesf = tmp_sf.get(); + } + ALWAYS_ERROR_IF(splinesf == 0, + "Requiringsurface to be a SplineSurface."); + SISLSurf* sislsf = GoSurf2SISL(*splinesf, false); + int dim = 3; + double epsco = 1e-15; // Not used +// double epsge = 1e-6; + int numintpt = 0; + double* pointpar = 0; + int numintcr = 0; + SISLIntcurve** intcurves = 0; + int stat = 0; + // Find the topology of the intersection + s1852(sislsf, pnt.begin(), radius, dim, epsco, geom_tol, + &numintpt, &pointpar, &numintcr, &intcurves, &stat); + // @@sbr Not sure this is the right solution. Maybe stat!=0 because of warning. + ALWAYS_ERROR_IF(stat!=0, + "s1852 returned code: " << stat); + // pointpar is not used any further + free(pointpar); + double maxstep = 0.0; + int makecurv = 2; // Make both geometric and parametric curves + int graphic = 0; // Do not draw the curve +// epsge = tol_.neighbour; + for (int i = 0; i < numintcr; ++i) { + // March out the intersection curves + s1315(sislsf,pnt.begin(), radius, dim, epsco, geom_tol, + maxstep, intcurves[i], makecurv, graphic, &stat); + SISLCurve* sc = intcurves[i]->pgeom; + if (sc == 0) { + MESSAGE("s1315 returned code: " << stat << ", returning."); + continue; + // freeIntcrvlist(intcurves, numintcr); + // freeSurf(sislsf); + // return curves; + } + double* t = sc->et; + double* c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + // Convert the geometric curve to Go format + shared_ptr gcv(new SplineCurve(sc->in, sc->ik, t, c, 3)); + sc = intcurves[i]->ppar1; + t = sc->et; + c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + shared_ptr pcv(new SplineCurve(sc->in, sc->ik, t, c, 2)); + // We prefer parameter curves. + // curves.push_back(shared_ptr + // (new CurveOnSurface(surf, pcv, gcv, true))); + curves.push_back(shared_ptr + (new CurveOnSurface(surf, pcv, gcv, false))); + + // // We would like the curve to have direction in consistency with trimming. + // vector derivs(2); + // gcv->point(derivs, gcv->startparam(), 1); + // Point par_pt = pcv->point(pcv->startparam()); + // Point surf_normal; + // surf->normal(surf_normal, par_pt[0], par_pt[1]); + // Point dir_ind = surf_normal%derivs[1]; // Cross product. + // // If angle(dir_ind, normal) < 90, turn curve (should be in (90, 180). + // double inner_product = dir_ind*vec; + // if (inner_product > 0) + // curves[curves.size()-1]->reverseParameterDirection(); + } + freeIntcrvlist(intcurves, numintcr); + freeSurf(sislsf); + + return curves; +} + +//=========================================================================== +vector > +BoundedUtils::intersectWithCone(shared_ptr& surf, + Point pnt, Point axispt, Point sfpt, + double geom_tol) +//=========================================================================== +{ + vector > curves; + + // Convert the surface to a SISLSurf in order to use SISL functions + // on it. The "false" argument dictates that the SISLSurf will only + // copy pointers to arrays, not the arrays themselves. + shared_ptr tmp_sf; + SplineSurface *splinesf = surf->getSplineSurface(); + if (!splinesf) + { + tmp_sf = shared_ptr(surf->asSplineSurface()); + splinesf = tmp_sf.get(); + } + ALWAYS_ERROR_IF(splinesf == 0, + "Requiringsurface to be a SplineSurface."); + SISLSurf* sislsf = GoSurf2SISL(*splinesf, false); + int dim = 3; + double epsco = 1e-15; // Not used +// double epsge = 1e-6; + int numintpt; + double* pointpar = 0; + int numintcr; + SISLIntcurve** intcurves = 0; + int stat; + // Find the topology of the intersection + s1854(sislsf, pnt.begin(), axispt.begin(), sfpt.begin(), dim, epsco, + geom_tol, &numintpt, &pointpar, &numintcr, &intcurves, &stat); + // @@sbr Not sure this is the right solution. Maybe stat!=0 because of warning. + ALWAYS_ERROR_IF(stat!=0, + "s1853 returned code: " << stat); + // pointpar is not used any further + free(pointpar); + double maxstep = 0.0; + int makecurv = 2; // Make both geometric and parametric curves + int graphic = 0; // Do not draw the curve +// epsge = tol_.neighbour; + for (int i = 0; i < numintcr; ++i) { + // March out the intersection curves + s1317(sislsf, pnt.begin(), axispt.begin(), sfpt.begin(), dim, epsco, + geom_tol, maxstep, intcurves[i], makecurv, graphic, &stat); + SISLCurve* sc = intcurves[i]->pgeom; + if (sc == 0) { + MESSAGE("s1317 returned code: " << stat << ", returning."); + continue; + // freeIntcrvlist(intcurves, numintcr); + // freeSurf(sislsf); + // return curves; + } + double* t = sc->et; + double* c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + // Convert the geometric curve to Go format + shared_ptr gcv(new SplineCurve(sc->in, sc->ik, t, c, 3)); + sc = intcurves[i]->ppar1; + t = sc->et; + c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + shared_ptr pcv(new SplineCurve(sc->in, sc->ik, t, c, 2)); + // We prefer parameter curves. + // curves.push_back(shared_ptr + // (new CurveOnSurface(surf, pcv, gcv, true))); + curves.push_back(shared_ptr + (new CurveOnSurface(surf, pcv, gcv, false))); + + // We would like the curve to have direction in consistency with trimming. + vector derivs(2); + gcv->point(derivs, gcv->startparam(), 1); + Point par_pt = pcv->point(pcv->startparam()); + Point surf_normal; + surf->normal(surf_normal, par_pt[0], par_pt[1]); + Point dir_ind = surf_normal%derivs[1]; // Cross product. + // If angle(dir_ind, normal) < 90, turn curve (should be in (90, 180). + Point vec = axispt - pnt; + double inner_product = dir_ind*vec; + if (inner_product > 0) + curves[curves.size()-1]->reverseParameterDirection(); + } + freeIntcrvlist(intcurves, numintcr); + freeSurf(sislsf); + + return curves; +} + +//=========================================================================== +vector > +BoundedUtils::intersectWithTorus(shared_ptr& surf, + Point pnt, Point normal, double rad1, + double rad2, double geom_tol) +//=========================================================================== +{ + vector > curves; + + // Convert the surface to a SISLSurf in order to use SISL functions + // on it. The "false" argument dictates that the SISLSurf will only + // copy pointers to arrays, not the arrays themselves. + shared_ptr tmp_sf; + SplineSurface *splinesf = surf->getSplineSurface(); + if (!splinesf) + { + tmp_sf = shared_ptr(surf->asSplineSurface()); + splinesf = tmp_sf.get(); + } + ALWAYS_ERROR_IF(splinesf == 0, + "Requiringsurface to be a SplineSurface."); + SISLSurf* sislsf = GoSurf2SISL(*splinesf, false); + int dim = 3; + double epsco = 1e-15; // Not used +// double epsge = 1e-6; + int numintpt; + double* pointpar = 0; + int numintcr; + SISLIntcurve** intcurves = 0; + int stat; + // Find the topology of the intersection + s1369(sislsf, pnt.begin(), normal.begin(), rad1, rad2, dim, epsco, + geom_tol, &numintpt, &pointpar, &numintcr, &intcurves, &stat); + // @@sbr Not sure this is the right solution. Maybe stat!=0 because of warning. + ALWAYS_ERROR_IF(stat!=0, + "s1853 returned code: " << stat); + // pointpar is not used any further + free(pointpar); + double maxstep = 0.0; + int makecurv = 2; // Make both geometric and parametric curves + int graphic = 0; // Do not draw the curve +// epsge = tol_.neighbour; + for (int i = 0; i < numintcr; ++i) { + // March out the intersection curves + s1318(sislsf,pnt.begin(), normal.begin(), rad1, rad2, dim, epsco, + geom_tol, maxstep, intcurves[i], makecurv, graphic, &stat); + SISLCurve* sc = intcurves[i]->pgeom; + if (sc == 0) { + MESSAGE("s1318 returned code: " << stat << ", returning."); + continue; + // freeIntcrvlist(intcurves, numintcr); + // freeSurf(sislsf); + // return curves; + } + double* t = sc->et; + double* c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + // Convert the geometric curve to Go format + shared_ptr gcv(new SplineCurve(sc->in, sc->ik, t, c, 3)); + sc = intcurves[i]->ppar1; + t = sc->et; + c = (sc->ikind==2 || sc->ikind==4)? sc->rcoef : sc->ecoef; + shared_ptr pcv(new SplineCurve(sc->in, sc->ik, t, c, 2)); + // We prefer parameter curves. + // curves.push_back(shared_ptr + // (new CurveOnSurface(surf, pcv, gcv, true))); + curves.push_back(shared_ptr + (new CurveOnSurface(surf, pcv, gcv, false))); + + // We would like the curve to have direction in consistency with trimming. + vector derivs(2); + gcv->point(derivs, gcv->startparam(), 1); + Point par_pt = pcv->point(pcv->startparam()); + Point surf_normal; + surf->normal(surf_normal, par_pt[0], par_pt[1]); + Point dir_ind = surf_normal%derivs[1]; // Cross product. + // If angle(dir_ind, normal) < 90, turn curve (should be in (90, 180). + double inner_product = dir_ind*normal; + if (inner_product > 0) + curves[curves.size()-1]->reverseParameterDirection(); + } + freeIntcrvlist(intcurves, numintcr); + freeSurf(sislsf); + + return curves; +} + +//=========================================================================== +void +BoundedUtils::getIntersectionCurveElem(shared_ptr& sf1, + shared_ptr& sf2, + vector >& int_segments1, + vector >& int_segments2, + double epsgeo) +//=========================================================================== +{ + // @@sbr epsgeo should be close to 1e-06; + ALWAYS_ERROR_IF((int_segments2.size() != 0) || (int_segments2.size() != 0), + "Segment vectors must be empty!"); + + shared_ptr elem1 = + dynamic_pointer_cast(sf1); + shared_ptr elem2 = + dynamic_pointer_cast(sf2); + + // Check if the second surface is closed if elementary + bool close_u1=false, close_v1=false; + bool close_u=false, close_v=false; + if (elem1.get()) + elem1->isClosed(close_u1, close_v1); + if (elem2.get()) + elem2->isClosed(close_u, close_v); + + bool not_swap = false; + if (elem1.get() && sf1->instanceType() != Class_Plane && close_u1 == false && + close_v1 == false) + not_swap = true; + // if (sf2->instanceType() == Class_SplineSurface || + // (sf2->instanceType() == Class_Plane && not_swap == false && + // (sf1->instanceType() != Class_SplineSurface && sf1->instanceType() != Class_Plane))) + if ((sf2->instanceType() == Class_SplineSurface || + sf1->instanceType() == Class_Plane) && + sf1->instanceType() != sf2->instanceType()) + { + getIntersectionCurveElem(sf2, sf1, int_segments2, int_segments1, epsgeo); + return; + } + + shared_ptr surf1 = dynamic_pointer_cast(sf1); + if (!surf1.get()) + { + surf1 = shared_ptr(sf1->asSplineSurface()); + } + + if (!(surf1 && elem2.get())) + THROW("Unexpected surface types in intersection"); + + shared_ptr surf1_2 = surf1; + vector > int_seg1 = + intersectWithElementarySurface(surf1_2, elem2, epsgeo); + + const Point sf_epspar = SurfaceTools::getParEpsilon(*elem2, epsgeo); + const double epspar = std::min(sf_epspar[0], sf_epspar[1]); + const double angtol = 0.01; + if (int_seg1.size() > 0) + { + // Split intersection curve is it traverses across a seam + RectDomain dom = elem2->containingDomain(); + if (close_u) + { + // Pick seam curve + double tpar = dom.umin(); + vector > cvs = elem2->constParamCurves(tpar, false); + if (cvs.size() == 1) + { + size_t nsegs = int_seg1.size(); + for (size_t ki=0; kiisClosed(); + double start = int_seg1[ki]->startparam(); + double end = int_seg1[ki]->endparam(); + double par1, par2, dist; + Point ptc1, ptc2; + ClosestPoint::closestPtCurves(int_seg1[ki].get(), cvs[0].get(), + par1, par2, dist, ptc1, ptc2); + if (dist < epsgeo && par1-start > epspar && end-par1 > epspar) + { + // Split intersection curve + vector > sub_cvs = int_seg1[ki]->split(par1); + if (sub_cvs.size() == 2 && (!cv_close)) + { + // Check if the curves can be merged in the other end + vector der1(2), der2(2); + sub_cvs[1]->point(der1, sub_cvs[1]->endparam(), 1); + sub_cvs[0]->point(der2, sub_cvs[0]->startparam(), 1); + double dsub = der1[0].dist(der2[0]); + double asub = der1[1].angle(der2[1]); + if (dsub <= epsgeo && asub <= angtol) + { + shared_ptr tmp_cv(dynamic_pointer_cast(sub_cvs[1])->clone()); + double dist2; + tmp_cv->appendCurve(sub_cvs[0].get(), 0, dist2, false); + if (dist2 <= epsgeo) + { + int_seg1[ki] = tmp_cv; + sub_cvs.clear(); + } + } + } + + if (sub_cvs.size() > 1) + { + int_seg1[ki] = dynamic_pointer_cast(sub_cvs[0]); + for (size_t kj=1; kj(sub_cvs[kj])); + } + } + } + } + } + if (close_v) + { + double tpar = dom.vmin(); + vector > cvs = elem2->constParamCurves(tpar, true); + if (cvs.size() == 1) + { + size_t nsegs = int_seg1.size(); + for (size_t ki=0; kiisClosed(); + double start = int_seg1[ki]->startparam(); + double end = int_seg1[ki]->endparam(); + double par1, par2, dist; + Point ptc1, ptc2; + ClosestPoint::closestPtCurves(int_seg1[ki].get(), cvs[0].get(), + par1, par2, dist, ptc1, ptc2); + if (dist < epsgeo && par1-start > epspar && end-par1 > epspar) + { + // Split intersection curve + vector > sub_cvs = int_seg1[ki]->split(par1); + if (sub_cvs.size() == 2 && (!cv_close)) + { + // Check if the curves can be merged in the other end + vector der1(2), der2(2); + sub_cvs[1]->point(der1, sub_cvs[1]->endparam(), 1); + sub_cvs[0]->point(der2, sub_cvs[0]->startparam(), 1); + double dsub = der1[0].dist(der2[0]); + double asub = der1[1].angle(der2[1]); + if (dsub <= epsgeo && asub <= angtol) + { + shared_ptr tmp_cv(dynamic_pointer_cast(sub_cvs[1])->clone()); + double dist2; + tmp_cv->appendCurve(sub_cvs[0].get(), 0, dist2, false); + if (dist2 <= epsgeo) + { + int_seg1[ki] = tmp_cv; + sub_cvs.clear(); + } + } + } + + if (sub_cvs.size() > 1) + { + + int_seg1[ki] = dynamic_pointer_cast(sub_cvs[0]); + for (size_t kj=1; kj(sub_cvs[kj])); + } + } + } + } + } + + vector > int_seg2(int_seg1.size()); + for (size_t ki=0; kihasSpaceCurve()) + THROW("Non-existing intersection curve in geometry space"); + + // Create intersection curve for surface two + int_seg2[ki] = shared_ptr(new CurveOnSurface(sf2, + int_seg1[ki]->spaceCurve(), + false)); + int_seg2[ki]->ensureParCrvExistence(epsgeo); + + // Ensure consistence for surface one + if (surf1.get() != sf1.get()) + { + if (sf1->instanceType() == Class_Plane) + int_seg1[ki] = shared_ptr(new CurveOnSurface(sf1, + int_seg1[ki]->parameterCurve(), + int_seg1[ki]->spaceCurve(), + false)); + else + { + int_seg1[ki] = shared_ptr(new CurveOnSurface(sf1, + int_seg1[ki]->spaceCurve(), + false)); + int_seg1[ki]->ensureParCrvExistence(epsgeo); + } + } + } + int_segments1 = int_seg1; + int_segments2 = int_seg2; + } +} + //=========================================================================== void BoundedUtils::getIntersectionCurve(shared_ptr& sf1, diff --git a/gotools-core/src/geometry/Circle.C b/gotools-core/src/geometry/Circle.C index 615d8be5c..42b89c32a 100644 --- a/gotools-core/src/geometry/Circle.C +++ b/gotools-core/src/geometry/Circle.C @@ -539,6 +539,8 @@ SplineCurve* Circle::createNonRationalSpline(double eps) const SplineCurve *crv = new SplineCurve(sislcrv->in, sislcrv->ik, sislcrv->et, sislcrv->ecoef, centre_.dimension()); crv->setParameterInterval(startparam_, endparam_); + if (isReversed_) + crv->reverseParameterDirection(); freeCurve(sislcrv); return crv; @@ -566,8 +568,8 @@ void Circle::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== -void Circle::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) +void Circle::appendCurve(ParamCurve* cv, int continuity, double& dist, + bool reparam, double tol) //=========================================================================== { // Check input @@ -686,7 +688,7 @@ void Circle::closestPoint(const Point& pt, // If input is on the "centre line", we arbitrarily assign the // point with t = tmin. Point vec = pt - centre_; - Point tmp = vec.cross(normal_); + Point tmp = (dimension() == 3) ? vec.cross(normal_) : vec; if (tmp.length() == 0.0) { clo_t = (seed != NULL) ? *seed : tmin; clo_pt = centre_ + radius_*(cos(clo_t)*vec1_ + sin(clo_t)*vec2_); @@ -800,33 +802,33 @@ void Circle::setParamBounds(double startpar, double endpar) //=========================================================================== { if (fabs(startpar) < ptol_) - startpar = 0.0; + startpar = 0.0; else if (fabs(2.0*M_PI-startpar) < ptol_) startpar = 2.0*M_PI; if (fabs(endpar) < ptol_) - endpar = 0.0; + endpar = 0.0; else if (fabs(2.0*M_PI-endpar) < ptol_) endpar = 2.0*M_PI; - if (startpar > -2.0 * M_PI - ptol_ && startpar < -2.0 * M_PI) - startpar = -2.0 * M_PI; - if (endpar < 2.0 * M_PI + ptol_ && endpar >2.0 * M_PI) - endpar = 2.0 * M_PI; - if (startpar >= endpar) - THROW("First parameter must be strictly less than second."); - if (startpar < -2.0 * M_PI || endpar > 2.0 * M_PI) - THROW("Parameters must be in [-2pi, 2pi]."); - if (endpar - startpar > 2.0 * M_PI) - THROW("(endpar - startpar) must not exceed 2pi."); - - double start = - parbound1_ + (startpar-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_); - double end = - parbound1_ + (endpar-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_); - parbound1_ = startpar; - parbound2_ = endpar; - startparam_ = start; - endparam_ = end; + if (startpar > -2.0 * M_PI - ptol_ && startpar < -2.0 * M_PI) + startpar = -2.0 * M_PI; + if (endpar < 2.0 * M_PI + ptol_ && endpar >2.0 * M_PI) + endpar = 2.0 * M_PI; + if (startpar >= endpar) + THROW("First parameter must be strictly less than second."); + if (startpar < -2.0 * M_PI || endpar > 2.0 * M_PI) + THROW("Parameters must be in [-2pi, 2pi]."); + if (endpar - startpar > 2.0 * M_PI) + THROW("(endpar - startpar) must not exceed 2pi."); + + double start = + parbound1_ + (startpar-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_); + double end = + parbound1_ + (endpar-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_); + parbound1_ = startpar; + parbound2_ = endpar; + startparam_ = start; + endparam_ = end; } //=========================================================================== diff --git a/gotools-core/src/geometry/Cone.C b/gotools-core/src/geometry/Cone.C index 48c4e6af2..8714baa10 100644 --- a/gotools-core/src/geometry/Cone.C +++ b/gotools-core/src/geometry/Cone.C @@ -696,10 +696,16 @@ void Cone::closestPoint(const Point& pt, // Bottom - a circle rad = radius_ + vmin * tan(cone_angle_); loc = location_ + vmin * z_axis_; - circle = Circle(rad, loc, z_axis_, x_axis_); + circle = Circle(fabs(rad), loc, z_axis_, x_axis_); circle.setParamBounds(umin, umax); circle.closestPoint(pt, umin, umax, clo_u, clo_pt, clo_dist); clo_v = vmin; + if (rad < 0) + { + clo_u += M_PI; + if (clo_u > umax) + clo_u -= (2.0*M_PI); + } if (clo_dist < epsilon) { clo_u = domain_.umin() + (clo_u - parbound_.umin())*(domain_.umax()-domain_.umin())/(parbound_.umax()-parbound_.umin()); @@ -713,10 +719,16 @@ void Cone::closestPoint(const Point& pt, // Top - a circle rad = radius_ + vmax * tan(cone_angle_); loc = location_ + vmax * z_axis_; - circle = Circle(rad, loc, z_axis_, x_axis_); + circle = Circle(fabs(rad), loc, z_axis_, x_axis_); circle.setParamBounds(umin, umax); circle.closestPoint(pt, umin, umax, tmp_clo_u, tmp_clo_pt, tmp_clo_dist); tmp_clo_v = vmax; + if (rad < 0) + { + tmp_clo_u += M_PI; + if (tmp_clo_u > umax) + tmp_clo_u -= (2.0*M_PI); + } if (tmp_clo_dist < clo_dist) { clo_u = tmp_clo_u; clo_v = tmp_clo_v; @@ -890,7 +902,7 @@ void Cone::getDegenerateParam(double& par, int& dir) const double eps = 1.0e-10; par = 0.0; double ang_tan = tan(cone_angle_); - if (ang_tan < eps) + if (fabs(ang_tan) < eps) dir = 0; else { @@ -1340,6 +1352,8 @@ Cone::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, const Point* start_par_pt, const Point* end_par_pt) const //=========================================================================== { + double eps = 1.0e-9; + // Default is not simple elementary parameter curve exists shared_ptr dummy; @@ -1405,7 +1419,7 @@ Cone::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, if (isClosed(dummy_u, dummy_v)) { // Extra check at the seam - double ptol = 1.0e-4; + double ptol = (tol < 0.001) ? 1.0e-4 : 1.0e-2; if (u1 < ptol) u1 = 2.0*M_PI; else if (u3 > 2.0*M_PI - ptol) @@ -1447,8 +1461,19 @@ Cone::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, } } + double domain_max = (idx == 1) ? parbound_.umax() : parbound_.vmax(); + bool isreversed = false; + if (closed && fabs(domain_max-par1[1-idx]) < eps) + { + double domain_min = (idx == 1) ? parbound_.umin() : parbound_.vmin(); + par1[1-idx] = domain_min; + par2[1-idx] = domain_max; + isreversed = true; + } + shared_ptr param_cv(new Line(par1, par2, - space_crv->startparam(), space_crv->endparam())); + space_crv->startparam(), space_crv->endparam(), + isreversed)); // TEST Point p1 = param_cv->ParamCurve::point(param_cv->startparam()); diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index a3b7637fe..4a3c0f9a9 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -58,6 +58,7 @@ #include #include + using namespace Go; using std::vector; using std::max; @@ -149,7 +150,12 @@ CurveOnSurface::CurveOnSurface(shared_ptr surf, Point pnt6 = curve->point(t4); Point pnt7 = surf->point(pt3[0], pt3[1]); Point pnt8 = surf->point(pt4[0], pt4[1]); - if (pnt5.dist(pnt7) + pnt6.dist(pnt8) > + if (constdir_ && (pt4-pt3)*(pt2-pt1) < 0.0) + { + same_orientation_ = false; + std::swap(pt1,pt2); + } + else if (pnt5.dist(pnt7) + pnt6.dist(pnt8) > pnt5.dist(pnt8) + pnt6.dist(pnt7)) { same_orientation_ = false; @@ -808,15 +814,16 @@ void CurveOnSurface::closestPoint(const Point& pt, std::vector > par_and_dist; // More values may be push_backed here in order to make a // better starting guess - int nmb_sample_pts = 3; + int psample = 0, ssample = 0; if ((pcurve_.get() != 0) && (pcurve_->instanceType() == Class_SplineCurve)) { shared_ptr tempsc(dynamic_pointer_cast(pcurve_)); - nmb_sample_pts = max(nmb_sample_pts, tempsc->numCoefs()); + psample = tempsc->numCoefs(); } if ((spacecurve_.get() != 0) && (spacecurve_->instanceType() == Class_SplineCurve)) { - nmb_sample_pts = max(nmb_sample_pts, - (dynamic_pointer_cast(spacecurve_))->numCoefs()); + ssample = dynamic_pointer_cast(spacecurve_)->numCoefs(); } + int nmb_sample_pts = (prefer_parameter_ && psample > 0) ? psample : ssample; + nmb_sample_pts = std::min(std::max(3, nmb_sample_pts), 50); par_and_dist.reserve(nmb_sample_pts); // VSK. 0517. Avoid guess points in the curve ends @@ -864,7 +871,8 @@ void CurveOnSurface::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== void CurveOnSurface::appendCurve(ParamCurve* other_curve, - int continuity, double& dist, bool reparam) + int continuity, double& dist, bool reparam, + double tol) //=========================================================================== { CurveOnSurface* other_cv = dynamic_cast(other_curve); @@ -887,7 +895,7 @@ void CurveOnSurface::appendCurve(ParamCurve* other_curve, "Mismatch between curve represetations of the two objects."); - double tol = 1.0e-4; + //double tol = 1.0e-4; if (prefer_parameter_ && (pcurve_.get() != NULL)) { // Save input curves @@ -1489,6 +1497,13 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, elem_cv = shared_ptr(elem_cv2->clone()); elem_cv->setParamBounds(bd_cv->startparam(), bd_cv->endparam()); } + if (!elem_cv) + { + double tol = epsgeo; + Point dir; + bool linear = spacecurve_->isLinear(dir, epsgeo); + int stop_break = 1; + } } if (elem_sf.get() && elem_cv.get()) @@ -1508,47 +1523,47 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, { Point spacecv_startpt = spacecurve_->point(spacecurve_->startparam()); Point spacecv_endpt = spacecurve_->point(spacecurve_->endparam()); - Point lifted_startpt = surface_->point(startpt[0], startpt[1]); - Point lifted_endpt = surface_->point(endpt[0], endpt[1]); - double dist_start = spacecv_startpt.dist(lifted_startpt); - double dist_end = spacecv_endpt.dist(lifted_endpt); - if (dist_start > epspar || dist_end > epspar) - { - MESSAGE("Mismatch between space_cv and proj end pts, dist_start: " - << dist_start << ", dist_end: " << dist_end); + // Point lifted_startpt = surface_->point(startpt[0], startpt[1]); + // Point lifted_endpt = surface_->point(endpt[0], endpt[1]); + // double dist_start = spacecv_startpt.dist(lifted_startpt); + // double dist_end = spacecv_endpt.dist(lifted_endpt); + // if (dist_start > epspar || dist_end > epspar) + // { + // MESSAGE("Mismatch between space_cv and proj end pts, dist_start: " + // << dist_start << ", dist_end: " << dist_end); // We write to file the iso-curves in these parameters. - std::ofstream debug_sf("tmp/under_sf.g2"); + std::ofstream debug_sf("under_sf.g2"); surface_->writeStandardHeader(debug_sf); surface_->write(debug_sf); - std::ofstream debug("tmp/iso_cvs.g2"); + std::ofstream debug("iso_cvs.g2"); spacecurve_->writeStandardHeader(debug); spacecurve_->write(debug); - vector pts; - pts.insert(pts.end(), spacecv_startpt.begin(), spacecv_startpt.end()); - pts.insert(pts.end(), lifted_startpt.begin(), lifted_startpt.end()); - pts.insert(pts.end(), spacecv_endpt.begin(), spacecv_endpt.end()); - pts.insert(pts.end(), lifted_endpt.begin(), lifted_endpt.end()); - PointCloud3D pt_cl(pts.begin(), 4); - pt_cl.writeStandardHeader(debug); - pt_cl.write(debug); - vector iso_par(4); - iso_par[0] = startpt[0]; - iso_par[1] = startpt[1]; - iso_par[2] = endpt[0]; - iso_par[3] = endpt[1]; - for (size_t ki = 0; ki < iso_par.size(); ++ki) - { - bool pardir_is_u = (ki%2 == 1); // Direction of moving parameter. - vector > const_cvs = surface_->constParamCurves(iso_par[ki], pardir_is_u); - for (size_t kj = 0; kj < const_cvs.size(); ++kj) - { - const_cvs[kj]->writeStandardHeader(debug); - const_cvs[kj]->write(debug); - } - } + // vector pts; + // pts.insert(pts.end(), spacecv_startpt.begin(), spacecv_startpt.end()); + // pts.insert(pts.end(), lifted_startpt.begin(), lifted_startpt.end()); + // pts.insert(pts.end(), spacecv_endpt.begin(), spacecv_endpt.end()); + // pts.insert(pts.end(), lifted_endpt.begin(), lifted_endpt.end()); + // PointCloud3D pt_cl(pts.begin(), 4); + // pt_cl.writeStandardHeader(debug); + // pt_cl.write(debug); + // vector iso_par(4); + // iso_par[0] = startpt[0]; + // iso_par[1] = startpt[1]; + // iso_par[2] = endpt[0]; + // iso_par[3] = endpt[1]; + // for (size_t ki = 0; ki < iso_par.size(); ++ki) + // { + // bool pardir_is_u = (ki%2 == 1); // Direction of moving parameter. + // vector > const_cvs = surface_->constParamCurves(iso_par[ki], pardir_is_u); + // for (size_t kj = 0; kj < const_cvs.size(); ++kj) + // { + // const_cvs[kj]->writeStandardHeader(debug); + // const_cvs[kj]->write(debug); + // } + // } double debug_val = 0.0; - } + // } } #endif @@ -2365,7 +2380,19 @@ int CurveOnSurface::whichBoundary(double tol, bool& same_orientation) const if (constdir_ == 1 || constdir_ == 2) { same_orientation = same_orientation_; - return at_bd_; + if (at_bd_ < 0) + { + BoundingBox box2D = pcurve_->boundingBox(); + Point low = box2D.low(); + Point high = box2D.high(); + RectDomain domain = surface_->containingDomain(); + Array p1(low[0], low[1]); + Array p2(high[0], high[1]); + int bd = domain.whichBoundary(p1, p2, tol); + return bd; + } + else + return at_bd_; } else { diff --git a/gotools-core/src/geometry/Cylinder.C b/gotools-core/src/geometry/Cylinder.C index 7f733e32a..8dd1e8828 100644 --- a/gotools-core/src/geometry/Cylinder.C +++ b/gotools-core/src/geometry/Cylinder.C @@ -649,12 +649,62 @@ void Cylinder::closestBoundaryPoint(const Point& pt, double *seed) const //=========================================================================== { - // This is a bit like cheating... + if (!isBounded()) + { + MESSAGE("closestBoundaryPoints does not make sense for unbounded surfaces"); + clo_dist = -1.0; + return; + } + + RectDomain domain = containingDomain(); + if (!rd) + rd = &domain; + shared_ptr bdcrv; + clo_dist = 1.0e10; // Initialize with a large number + Point cpt; + double cdist, cpar; + + // Checking closest point on the bottom boundary + bdcrv = constParamCurve(rd->vmin(), true, rd->umin(), rd->umax()); + bdcrv->closestPoint(pt, rd->umin(), rd->umax(), + clo_u, clo_pt, clo_dist, seed); + clo_v = rd->vmin(); + + // Checking the right boundary + bdcrv = constParamCurve(rd->umax(), false, rd->vmin(), rd->vmax()); + bdcrv->closestPoint(pt, rd->vmin(), rd->vmax(), cpar, cpt, cdist, + (seed == 0) ? seed : seed+1); + if (cdist < clo_dist) + { + clo_pt = cpt; + clo_u = rd->umax(); + clo_v = cpar; + clo_dist = cdist; + } - SplineSurface* sf = geometrySurface(); - sf->closestBoundaryPoint(pt, clo_u, clo_v, clo_pt, clo_dist, epsilon, - rd, seed); - delete sf; + // Checking the upper boundary + bdcrv = constParamCurve(rd->vmax(), true, rd->umin(), rd->umax()); + bdcrv->closestPoint(pt, rd->umin(), rd->umax(), cpar, cpt, cdist, + seed); + if (cdist < clo_dist) + { + clo_pt = cpt; + clo_u = cpar; + clo_v = rd->vmax(); + clo_dist = cdist; + } + + // Checking the left boundary + bdcrv = constParamCurve(rd->umin(), false, rd->vmin(), rd->vmax()); + bdcrv->closestPoint(pt, rd->vmin(), rd->vmax(), cpar, cpt, cdist, + (seed == 0) ? seed : seed+1); + if (cdist < clo_dist) + { + clo_pt = cpt; + clo_u = rd->umin(); + clo_v = cpar; + clo_dist = cdist; + } } @@ -688,6 +738,8 @@ Cylinder::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, const Point* start_par_pt, const Point* end_par_pt) const //=========================================================================== { + double eps = 1.0e-9; + // Default is not simple elementary parameter curve exists shared_ptr dummy; @@ -730,7 +782,7 @@ Cylinder::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, if (d1 > tol || d2 > tol) return dummy; - Point par1(2), par2(2); + Point par1(2), par2(2); par1[idx] = par2[idx] = 0.5*(parval1[idx] + parval2[idx]); par1[1-idx] = parval1[1-idx]; par2[1-idx] = parval2[1-idx]; diff --git a/gotools-core/src/geometry/Ellipse.C b/gotools-core/src/geometry/Ellipse.C index ba49f11bb..51564fb82 100644 --- a/gotools-core/src/geometry/Ellipse.C +++ b/gotools-core/src/geometry/Ellipse.C @@ -544,8 +544,8 @@ void Ellipse::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== -void Ellipse::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) +void Ellipse::appendCurve(ParamCurve* cv, int continuity, + double& dist, bool reparam, double tol) //=========================================================================== { // Check input diff --git a/gotools-core/src/geometry/FileUtils.C b/gotools-core/src/geometry/FileUtils.C index 3d90a5162..05bc02f48 100644 --- a/gotools-core/src/geometry/FileUtils.C +++ b/gotools-core/src/geometry/FileUtils.C @@ -133,7 +133,7 @@ void FileUtils::readTxtPointFile(std::ifstream& is, int del, for (int ki=1; ki> xx; - if (xx != ',') + if (xx != ',' && xx != ';') is.putback(xx); is >> tmp; data.push_back(tmp); diff --git a/gotools-core/src/geometry/GSCappendCurve.C b/gotools-core/src/geometry/GSCappendCurve.C index 8d984849d..cbafbcc9f 100644 --- a/gotools-core/src/geometry/GSCappendCurve.C +++ b/gotools-core/src/geometry/GSCappendCurve.C @@ -55,8 +55,8 @@ namespace Go{ //=========================================================================== -void SplineCurve::appendCurve(ParamCurve* other_curve, - int continuity, double& dist, bool repar) +void SplineCurve::appendCurve(ParamCurve* other_curve, int continuity, + double& dist, bool repar, double tol) //=========================================================================== { SplineCurve* other_cv = dynamic_cast(other_curve); diff --git a/gotools-core/src/geometry/GeometryTools.C b/gotools-core/src/geometry/GeometryTools.C index 1c807fb0c..9027050ed 100644 --- a/gotools-core/src/geometry/GeometryTools.C +++ b/gotools-core/src/geometry/GeometryTools.C @@ -1076,17 +1076,18 @@ void GeometryTools::translateSplineCurve(const Point& trans_vec, SplineCurve& cv //------------------------------------------------------------------- { int ki; - ASSERT(trans_vec.dimension() == 3); // We're working in 3D space. - int dim = 3 + cv.rational(); + int dim0 = trans_vec.dimension(); + ASSERT(dim0 == 2 || dim0 == 3); // We're working in 2D or 3D space. + int dim = dim0 + cv.rational(); std::vector::iterator iter = cv.rational() ? cv.rcoefs_begin() : cv.coefs_begin(); std::vector::iterator end_iter = cv.rational() ? cv.rcoefs_end() : cv.coefs_end(); std::vector::iterator coef_iter = cv.coefs_begin(); while (iter != end_iter) { // @@ A faster approach would be to use rotation matrix directly. - double weight = cv.rational() ? iter[3] : 1.0; - for (ki = 0; ki < 3; ++ki) + double weight = cv.rational() ? iter[dim0] : 1.0; + for (ki = 0; ki < dim0; ++ki) iter[ki] = (iter[ki]/weight + trans_vec[ki])*weight; if (cv.rational()) { - for (ki = 0; ki < 3; ++ki) + for (ki = 0; ki < dim0; ++ki) coef_iter[ki] = iter[ki]/weight; coef_iter += dim - 1; } @@ -1138,15 +1139,16 @@ void GeometryTools::rotateSplineCurve(Point rot_axis, double alpha, SplineCurve& //------------------------------------------------------------------- { rot_axis.normalize(); - ASSERT(rot_axis.dimension() == 3); // We're working in 3D space. - int dim = 3 + cv.rational(); + int dim0 = rot_axis.dimension(); + ASSERT(dim0 == 2 || dim0 == 3); // We're working in 2D or 3D space. + int dim = dim0 + cv.rational(); std::vector::iterator iter = cv.rational() ? cv.rcoefs_begin() : cv.coefs_begin(); std::vector::iterator end_iter = cv.rational() ? cv.rcoefs_end() : cv.coefs_end(); std::vector::iterator coef_iter = cv.coefs_begin(); while (iter != end_iter) { // @@ A faster approach would be to use rotation matrix directly. rotatePoint(rot_axis, alpha, &*iter); if (cv.rational()) { - for (int ki = 0; ki < 3; ++ki) + for (int ki = 0; ki < dim0; ++ki) coef_iter[ki] = iter[ki]/iter[3]; coef_iter += dim - 1; } @@ -1177,17 +1179,21 @@ void GeometryTools::rotateLineCloud(Point rot_axis, double alpha, LineCloud& lc) void GeometryTools::rotatePoint(Point rot_axis, double alpha, double* space_pt) //------------------------------------------------------------------- { + int dim = rot_axis.dimension(); + ASSERT(dim == 2 || dim == 3); + // We're working in 2D or 3D space. The axis content is used only in 3D + if (dim == 3) rot_axis.normalize(); - ASSERT(rot_axis.dimension() == 3); // We're working in 3D space. - std::vector rot_mat = GeometryTools::getRotationMatrix(rot_axis, alpha); - int ki, kj; - Point rotated_pt(0.0, 0.0, 0.0); - for (ki = 0; ki < 3; ++ki) - for (kj = 0; kj < 3; ++kj) - rotated_pt[ki] += rot_mat[ki*3+kj]*space_pt[kj]; - for (ki = 0; ki < 3; ++ki) - space_pt[ki] = rotated_pt[ki]; + std::vector rot_mat = GeometryTools::getRotationMatrix(rot_axis, alpha); + int ki, kj; + Point rotated_pt(dim); + for (ki = 0; ki < dim; ++ki) + for (kj = 0; kj < dim; ++kj) + rotated_pt[ki] += rot_mat[ki*dim+kj]*space_pt[kj]; + + for (ki = 0; ki < dim; ++ki) + space_pt[ki] = rotated_pt[ki]; } @@ -1203,25 +1209,38 @@ void GeometryTools::rotatePoint(Point rot_axis, double alpha, Point& space_pt) std::vector GeometryTools::getRotationMatrix(const Point& unit_rot_axis, double alpha) //------------------------------------------------------------------- { - ASSERT(unit_rot_axis.dimension() == 3); // We're working in 3D space. - std::vector return_matrix(9); - // Using the notation in 'Computational Gemoetry for Design and Manufacture'. - double u1 = unit_rot_axis[0]; - double u2 = unit_rot_axis[1]; - double u3 = unit_rot_axis[2]; - double cos_a = cos(alpha); - double sin_a = sin(alpha); - - return_matrix[0] = u1*u1 + cos_a*(1 - u1*u1); - return_matrix[1] = u1*u2*(1 - cos_a) - u3*sin_a; - return_matrix[2] = u3*u1*(1 - cos_a) + u2*sin_a; - return_matrix[3] = u1*u2*(1 - cos_a) + u3*sin_a; - return_matrix[4] = u2*u2 + cos_a*(1 - u2*u2); - return_matrix[5] = u2*u3*(1 - cos_a) - u1*sin_a; - return_matrix[6] = u3*u1*(1 - cos_a) - u2*sin_a; - return_matrix[7] = u2*u3*(1 - cos_a) + u1*sin_a; - return_matrix[8] = u3*u3 + cos_a*(1 - u3*u3); - + int dim = unit_rot_axis.dimension(); + ASSERT(dim == 2 || dim == 3); + // We're working in 2D or 3D space. The axis content is used only in 3D + std::vector return_matrix; + double cos_a = cos(alpha); + double sin_a = sin(alpha); + if (dim == 2) + { + return_matrix.resize(4); + return_matrix[0] = cos_a; + return_matrix[1] = -sin_a; + return_matrix[2] = sin_a; + return_matrix[3] = cos_a; + } + else + { + return_matrix.resize(9); + // Using the notation in 'Computational Gemoetry for Design and Manufacture'. + double u1 = unit_rot_axis[0]; + double u2 = unit_rot_axis[1]; + double u3 = unit_rot_axis[2]; + + return_matrix[0] = u1*u1 + cos_a*(1 - u1*u1); + return_matrix[1] = u1*u2*(1 - cos_a) - u3*sin_a; + return_matrix[2] = u3*u1*(1 - cos_a) + u2*sin_a; + return_matrix[3] = u1*u2*(1 - cos_a) + u3*sin_a; + return_matrix[4] = u2*u2 + cos_a*(1 - u2*u2); + return_matrix[5] = u2*u3*(1 - cos_a) - u1*sin_a; + return_matrix[6] = u3*u1*(1 - cos_a) - u2*sin_a; + return_matrix[7] = u2*u3*(1 - cos_a) + u1*sin_a; + return_matrix[8] = u3*u3 + cos_a*(1 - u3*u3); + } return return_matrix; } diff --git a/gotools-core/src/geometry/Hyperbola.C b/gotools-core/src/geometry/Hyperbola.C index 9e35b7bae..ad13e8694 100644 --- a/gotools-core/src/geometry/Hyperbola.C +++ b/gotools-core/src/geometry/Hyperbola.C @@ -435,7 +435,8 @@ void Hyperbola::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== void Hyperbola::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) + int continuity, double& dist, bool reparam, + double tol) //=========================================================================== { MESSAGE("appendCurve() not implemented!"); diff --git a/gotools-core/src/geometry/Line.C b/gotools-core/src/geometry/Line.C index c9294ac12..c31e98ec3 100644 --- a/gotools-core/src/geometry/Line.C +++ b/gotools-core/src/geometry/Line.C @@ -417,19 +417,20 @@ Line* Line::subCurve(double from_par, double to_par, { double start = endparam_ - (to_par - startparam_); double end = startparam_ + (endparam_ - from_par); + if (start > end) - { - std::swap(start, end); - } + { + std::swap(start, end); + } double bound1 = bounded ? parbound1_ + - (start-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_) : - start; + (start-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_) : + start; double bound2 = bounded ? parbound1_ + - (end-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_) : - end; + (end-startparam_)*(parbound2_-parbound1_)/(endparam_-startparam_) : + end; line->setParamBounds(bound1, bound2); if (from_par > to_par) - std::swap(from_par, to_par); + std::swap(from_par, to_par); line->setParameterInterval(from_par, to_par); } else @@ -467,7 +468,7 @@ void Line::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== void Line::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) + int continuity, double& dist, bool reparam, double tol) //=========================================================================== { // Check input diff --git a/gotools-core/src/geometry/Parabola.C b/gotools-core/src/geometry/Parabola.C index 4cbb6b408..1793144d1 100644 --- a/gotools-core/src/geometry/Parabola.C +++ b/gotools-core/src/geometry/Parabola.C @@ -428,7 +428,8 @@ void Parabola::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== void Parabola::appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam) + int continuity, double& dist, bool reparam, + double tol) //=========================================================================== { MESSAGE("appendCurve() not implemented!"); diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index d1670d0c8..25cf0fc57 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -39,6 +39,7 @@ #include "GoTools/geometry/Plane.h" #include "GoTools/geometry/Line.h" +#include "GoTools/geometry/Circle.h" #include "GoTools/geometry/SplineSurface.h" #include #include @@ -50,6 +51,7 @@ using std::numeric_limits; using std::streamsize; using std::swap; +#define DEBUG namespace Go { @@ -741,6 +743,193 @@ bool Plane::isDegenerate(bool& b, bool& r, return false; } +//=========================================================================== +shared_ptr +Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, + const Point* start_par_pt, const Point* end_par_pt) const +//=========================================================================== +{ + // Default is not simple elementary parameter curve exists + shared_ptr param_cv; + + // Bookkeeping related to swapped parameters + int ind1 = 0; + int ind2 = 1; + if (isSwapped()) + swap(ind1, ind2); + + double t1, t2; + int idx; + bool closed = false; + if (space_crv->instanceType() == Class_Line) + { + if (!((Line*)(space_crv))->isBounded()) + return param_cv; // Project endpoints onto the surface + t1 = space_crv->startparam(); + t2 = space_crv->endparam(); + idx = ind1; // 0 + } + else if (space_crv->instanceType() == Class_Circle) + { + t1 = space_crv->startparam(); + closed = ((Circle*)(space_crv))->isClosed(); + t2 = (closed) ? 0.5*(t1 + space_crv->endparam()) : + space_crv->endparam(); + idx = ind2; // 1 + } + else + return param_cv; + + double fac = (parbound_.umax()-parbound_.umin())/(domain_.umax()-domain_.umin()); + double parval1[2], parval2[2]; + double d1, d2; + Point close1, close2; + Point pos1 = space_crv->ParamCurve::point(t1); + Point pos2 = space_crv->ParamCurve::point(t2); + closestPoint(pos1, parval1[0], parval1[1], close1, d1, tol); + closestPoint(pos2, parval2[0], parval2[1], close2, d2, tol); + if (d1 > tol || d2 > tol) + return param_cv; + + Point par1(parval1[0], parval1[1]); + Point par2(parval2[0], parval2[1]); + if (space_crv->instanceType() == Class_Line) + { + Point pos = (t2*par1 - t1*par2)/(t2 - t1); + Point dir = (par2 - par1); + dir.normalize(); + + param_cv = shared_ptr(new Line(pos, dir)); + param_cv->setParamBounds(t1, t2); + } + else + { + Point mid = space_crv->ParamCurve::point(0.5*(t1+t2)); + double parval3[2]; + double d3; + Point close3; + closestPoint(mid, parval3[0], parval3[1], close3, d3, tol); + if (d3 > tol) + return param_cv; // Should not happen + + Point pmid(parval3[0], parval3[1]); + double alpha = t2 - t1; + double len = 0.5*par1.dist(par2); + double radius = len/sin(0.5*alpha); + Point vec = par1 + par2 - 2*pmid; + vec.normalize(); + Point centre = pmid + radius*vec; + Point param_cv_axis(0.0, 0.0); // A dimension two circle is not supported for every + // functionality + + Point centre_3d = ((Circle*)(space_crv))->getCentre(); + Point xvec_3d = ((Circle*)(space_crv))->getXAxis(); + Point vec1_3d = pos1 - centre_3d; + Point vec2_3d = pos2 - centre_3d; + double beta_3d = vec1_3d.angle(vec2_3d); + + Point vec1 = par1 - centre; + Point vec2 = par2 - centre; + double beta = vec1.angle(vec2); + + double gamma = 2*M_PI - t1; + double gamma2 = gamma; + if (gamma2 > 0.5*M_PI) + gamma2 -= 0.5*M_PI; + Point xvec(cos(gamma2)*(centre[0]-par1[0])-sin(gamma2)*(centre[1]-par1[1]), + sin(gamma2)*(centre[0]-par1[0])+cos(gamma2)*(centre[1]-par1[1])); + xvec.normalize(); + + Point xvec2(cos(t1)*(par1[0]-centre[0])+sin(t1)*(par1[1]-centre[1]), + -sin(t1)*(par1[0]-centre[0])+cos(t1)*(par1[1]-centre[1])); + xvec2.normalize(); + shared_ptr circ1(new Circle(radius, centre, param_cv_axis, + xvec, false)); + Point p1_1 = circ1->ParamCurve::point(t1); + Point p1_2 = circ1->ParamCurve::point(t2); + shared_ptr circ2(new Circle(radius, centre, param_cv_axis, + xvec2, false)); + Point p2_1 = circ2->ParamCurve::point(t1); + Point p2_2 = circ2->ParamCurve::point(t2); + Point p2_3 = circ2->ParamCurve::point(0.5*(t1+t2)); + + Point p3_1, p3_2; + point(p3_1, p2_1[0], p2_1[1]); + point(p3_2, p2_2[0], p2_2[1]); + double dd1 = pos1.dist(p3_1) + pos2.dist(p3_2); + if (dd1 > std::max(d1 + d2, tol)) + { + Point xvec3(cos(t2)*(par1[0]-centre[0])-sin(t2)*(par1[0]-centre[1]), + sin(t2)*(par1[1]-centre[0])+cos(t2)*(par1[1]-centre[1])); + xvec3.normalize(); + shared_ptr circ3(new Circle(radius, centre, param_cv_axis, + xvec3, true)); + Point p3_1 = circ3->ParamCurve::point(t1); + Point p3_2 = circ3->ParamCurve::point(t2); + Point p3_3 = circ3->ParamCurve::point(0.5*(t1+t2)); + + Point xvec4(xvec2[1],-xvec2[0]); + shared_ptr circ4(new Circle(radius, centre, param_cv_axis, + xvec4, false)); + Point p4_1 = circ4->ParamCurve::point(t1); + Point p4_2 = circ4->ParamCurve::point(t2); + Point p4_3 = circ4->ParamCurve::point(0.5*(t1+t2)); + + double t3 = space_crv->startparam() + space_crv->endparam() - t1; + Point xvec5(cos(t3)*(par1[0]-centre[0])-sin(t3)*(par1[0]-centre[1]), + sin(t3)*(par1[1]-centre[0])+cos(t3)*(par1[1]-centre[1])); + xvec5.normalize(); + shared_ptr circ5(new Circle(radius, centre, param_cv_axis, + xvec5, true)); + Point p5_1 = circ5->ParamCurve::point(t1); + Point p5_2 = circ5->ParamCurve::point(t2); + Point p5_3 = circ5->ParamCurve::point(0.5*(t1+t2)); + + + int stop_break_circ = 1; + } + + if (gamma2 < gamma) + { + std::swap(xvec[0], xvec[1]); + xvec[1] *= -1; + } + + double sgn1 = xvec_3d*vec2_3d; + double sgn2 = xvec*vec2; + bool reversed = false; + if (sgn1*sgn2 < 0.0) + { + // Opposite orientation of curve in geometry and parameter space. Redo xvec + // computation + gamma = 2*M_PI - t2; + gamma2 = gamma; + if (gamma2 > 0.5*M_PI) + gamma2 -= 0.5*M_PI; + xvec = Point(cos(gamma2)*(centre[0]-par1[0])-sin(gamma2)*(centre[1]-par1[1]), + sin(gamma2)*(centre[0]-par1[0])+cos(gamma2)*(centre[1]-par1[1])); + xvec.normalize(); + if (gamma2 < gamma) + { + std::swap(xvec[0], xvec[1]); + xvec[1] *= -1; + } + reversed = true; + } + param_cv = (dd1 <= tol) ? circ2 : + shared_ptr(new Circle(radius, centre, param_cv_axis, + xvec, reversed)); + param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); + } +#ifdef DEBUG + // TEST + Point p1 = param_cv->ParamCurve::point(param_cv->startparam()); + Point p2 = param_cv->ParamCurve::point(param_cv->endparam()); + int stop_break = 1; +#endif + + return param_cv; +} //=========================================================================== void Plane::getDegenerateCorners(vector& deg_corners, double tol) const @@ -980,6 +1169,7 @@ bool Plane::isBounded() const return false; return true; } + //=========================================================================== Plane* Plane::intersect(const RotatedBox& bd_box) const //=========================================================================== diff --git a/gotools-core/src/geometry/RectDomain.C b/gotools-core/src/geometry/RectDomain.C index eea7e0486..ad5466c7b 100644 --- a/gotools-core/src/geometry/RectDomain.C +++ b/gotools-core/src/geometry/RectDomain.C @@ -137,6 +137,47 @@ bool RectDomain::isOnBoundary(const Array& point, return false; } +//=========================================================================== +bool RectDomain::isOnBoundary(const Array& point, + double tolerance, int& bd, int& bd2) const +//=========================================================================== +{ + // Being on the boundary implies being in the domain + bd = bd2 = -1; + if (!isInDomain(point, tolerance)) + return false; + + if (point[0] < ll_[0] + tolerance) + { + bd = 0; + if (point[1] < ll_[1] + tolerance) + bd2 = 2; + else if (point[1] > ur_[1] - tolerance) + bd2 = 3; + return true; + } + if (point[0] > ur_[0] - tolerance) + { + bd = 1; + if (point[1] < ll_[1] + tolerance) + bd2 = 2; + else if (point[1] > ur_[1] - tolerance) + bd2 = 3; + return true; + } + if (point[1] < ll_[1] + tolerance) + { + bd = 2; + return true; + } + if (point[1] > ur_[1] - tolerance) + { + bd = 3; + return true; + } + return false; +} + //=========================================================================== bool RectDomain::isOnCorner(const Array& point, @@ -159,6 +200,7 @@ bool RectDomain::isOnCorner(const Array& point, } + //=========================================================================== int RectDomain::whichBoundary(const Array& point1, const Array& point2, @@ -173,10 +215,10 @@ int RectDomain::whichBoundary(const Array& point1, return 1; // umax if (fabs(point1[1] - ll_[1]) < tolerance && fabs(point2[1] - ll_[1]) < tolerance) - return 2; // umax + return 2; // vmin if (fabs(point1[1] - ur_[1]) < tolerance && fabs(point2[1] - ur_[1]) < tolerance) - return 3; // umax + return 3; // vmax return -1; // Not a common boundary } @@ -273,4 +315,19 @@ void RectDomain::intersectWith(const RectDomain& rd) // return true; // return false; } + +//=========================================================================== +bool RectDomain::overlap(const RectDomain& rd, double tol1, double tol2) +//=========================================================================== +{ + if (ll_[0] > rd.ur_[0]+tol1) + return false; + if (ur_[0] < rd.ll_[0]-tol1) + return false; + if (ll_[1] > rd.ur_[1]+tol2) + return false; + if (ur_[1] < rd.ll_[1]-tol2) + return false; + return true; +} } // namespace Go diff --git a/gotools-core/src/geometry/SplineDebugUtils.C b/gotools-core/src/geometry/SplineDebugUtils.C index eedf6ffdd..d03cc6d3b 100644 --- a/gotools-core/src/geometry/SplineDebugUtils.C +++ b/gotools-core/src/geometry/SplineDebugUtils.C @@ -77,6 +77,30 @@ void SplineDebugUtils::writeSpaceParamCurve(const SplineCurve& pcurve, std::ostr space_pcurve.write(os); } +//=========================================================================== +void SplineDebugUtils::writeSpace1DCurve(const SplineCurve& curve, std::ostream& os, double z) +//=========================================================================== +{ + ALWAYS_ERROR_IF(curve.dimension() != 1, + "Expecting input of 1D-curve."); + + std::vector coefs; + std::vector coefs_in(curve.coefs_begin(), curve.coefs_end()); + for (int i = 0; i < curve.numCoefs(); ++i) { + double gp = curve.basis().grevilleParameter(i); + coefs.push_back(gp); + coefs.push_back(coefs_in[i]); + coefs.push_back(z); + } + + SplineCurve space_curve = + SplineCurve(curve.numCoefs(), curve.order(), + curve.basis().begin(), coefs.begin(), 3); + space_curve.writeStandardHeader(os); + space_curve.write(os); +} + + //=========================================================================== void SplineDebugUtils::writeSpaceParamCurve(const Line& pline, std::ostream& os, double z) @@ -100,6 +124,52 @@ void SplineDebugUtils::writeSpaceParamCurve(const Line& pline, std::ostream& os, } +//=========================================================================== +void SplineDebugUtils::writeSpaceParamCurve(const Circle& pcirc, std::ostream& os, double z) +//=========================================================================== +{ + ALWAYS_ERROR_IF(pcirc.dimension() != 2, + "Expecting input of 2D-curve."); + + Point centre = pcirc.getCentre(); + Point vec1 = pcirc.getXAxis(); + Point normal(0.0, 0.0, 1.0); + Point centre_3d(centre[0], centre[1], z); + Point vec1_3d(vec1[0], vec1[1], 0.0); + double radius = pcirc.getRadius(); + Circle circ_3d(radius, centre_3d, normal, vec1_3d); + circ_3d.setParamBounds(pcirc.startparam(), pcirc.endparam()); + + circ_3d.writeStandardHeader(os); + circ_3d.write(os); +} + + +//=========================================================================== +void SplineDebugUtils::writeSpaceParamCurve(shared_ptr pcurve, + std::ostream& os, double z) +//=========================================================================== +{ + if (pcurve->instanceType() == Class_SplineCurve) + { + shared_ptr spline_cv = + dynamic_pointer_cast(pcurve); + writeSpaceParamCurve(*spline_cv, os, z); + } + else if (pcurve->instanceType() == Class_Line) + { + shared_ptr line = + dynamic_pointer_cast(pcurve); + writeSpaceParamCurve(*line, os, z); + } + else if (pcurve->instanceType() == Class_Circle) + { + shared_ptr circle = + dynamic_pointer_cast(pcurve); + writeSpaceParamCurve(*circle, os, z); + } +} + //=========================================================================== void SplineDebugUtils::writeTrimmedInfo(BoundedSurface& bd_sf, std::ostream& os, double z) diff --git a/gotools-core/src/geometry/SurfaceTools.C b/gotools-core/src/geometry/SurfaceTools.C index 7c6bcfb48..6203ab283 100644 --- a/gotools-core/src/geometry/SurfaceTools.C +++ b/gotools-core/src/geometry/SurfaceTools.C @@ -186,7 +186,8 @@ SurfaceTools::absolutelyAllBoundarySfLoops(shared_ptr surf, // Use a negative degeneracy tolarance to tell that also degenerate // boundaries must be included in the loop std::vector cvloopvec; - cvloopvec.push_back(SurfaceTools::outerBoundarySfLoop(surf, 0.0)); + //cvloopvec.push_back(SurfaceTools::outerBoundarySfLoop(surf, 0.0)); + cvloopvec.push_back(SurfaceTools::outerBoundarySfLoop(surf, -0.1)); return cvloopvec; } } diff --git a/gotools-core/src/geometry/Torus.C b/gotools-core/src/geometry/Torus.C index 84aef2769..42e688d60 100644 --- a/gotools-core/src/geometry/Torus.C +++ b/gotools-core/src/geometry/Torus.C @@ -796,7 +796,25 @@ void Torus::setParameterBounds(double from_upar, double from_vpar, getOrientedParameters(from_upar, from_vpar); getOrientedParameters(to_upar, to_vpar); - // NOTE: If parameters are swapped, from_upar and from_vpar are swapped. + if (fabs(from_upar) < ptol_) + from_upar = 0.0; + else if (fabs(2.0*M_PI-from_upar) < ptol_) + from_upar = 2.0*M_PI; + if (fabs(from_vpar) < ptol_) + from_vpar = 0.0; + else if (fabs(2.0*M_PI-from_vpar) < ptol_) + from_vpar = 2.0*M_PI; + + if (fabs(to_upar) < ptol_) + to_upar = 0.0; + else if (fabs(2.0*M_PI-to_upar) < ptol_) + to_upar = 2.0*M_PI; + if (fabs(to_vpar) < ptol_) + to_vpar = 0.0; + else if (fabs(2.0*M_PI-to_vpar) < ptol_) + to_vpar = 2.0*M_PI; + + // NOTE: If parameters are swapped, from_upar and from_vpar are swapped. // Ditto for to_upar/to_vpar. if (from_upar > -2.0 * M_PI - ptol_ && from_upar < -2.0 * M_PI) from_upar = -2.0 * M_PI; diff --git a/gotools-core/src/geometry/closestPtCurves.C b/gotools-core/src/geometry/closestPtCurves.C index e7289ff2e..641436b4c 100644 --- a/gotools-core/src/geometry/closestPtCurves.C +++ b/gotools-core/src/geometry/closestPtCurves.C @@ -42,6 +42,7 @@ using std::vector; #include "GoTools/geometry/GeometryTools.h" #include "GoTools/geometry/ClosestPoint.h" +#include "GoTools/geometry/CurveOnSurface.h" #include "GoTools/utils/Values.h" // MAXDOUBLE //*************************************************************************** @@ -87,6 +88,9 @@ namespace void computeSeedCvCv(const SplineCurve* pc1, const SplineCurve* pc2, double& seed1, double& seed2); +void computeSeedCvCv2(const ParamCurve* cv1, const ParamCurve* cv2, + double& seed1, double& seed2); + void insideParamDomain(double& delta, double acoef, double astart, double aend); @@ -127,8 +131,9 @@ void ClosestPoint::closestPtCurves(const ParamCurve* cv1, const ParamCurve* cv2, computeSeedCvCv(pc1, pc2, seed1, seed2); } else { - seed1 = 0.5*(tmin1+tmax1); - seed2 = 0.5*(tmin2+tmax2); + computeSeedCvCv2(cv1, cv2, seed1, seed2); + // seed1 = 0.5*(tmin1+tmax1); + // seed2 = 0.5*(tmin2+tmax2); } // Iterate for closest point @@ -370,6 +375,120 @@ void computeSeedCvCv(const SplineCurve* cv1, const SplineCurve* cv2, } +//*************************************************************************** +void computeSeedCvCv2(const ParamCurve* cv1, const ParamCurve* cv2, + double& seed1, double& seed2) +//*************************************************************************** +{ + const int dim = cv1->dimension(); + DEBUG_ERROR_IF(dim!=cv2->dimension(), "Dimension mismatch."); + + // Make guess point to the iteration. + // Find position of closest vertices + std::vector::const_iterator co1; + std::vector::const_iterator co2; + std::vector::const_iterator co3; + std::vector::const_iterator co12; + std::vector::const_iterator co22; + + const SplineCurve *pc1 = dynamic_cast(cv1); + const CurveOnSurface *sfc1 = dynamic_cast(cv1); + if (sfc1) + pc1 = dynamic_cast(sfc1->spaceCurve().get()); + const SplineCurve *pc2 = dynamic_cast(cv2); + const CurveOnSurface *sfc2 = dynamic_cast(cv2); + if (sfc2) + pc2 = dynamic_cast(sfc2->spaceCurve().get()); + + int num_sample = 5; + vector pts1; + vector pts2; + vector par1(num_sample); + vector par2(num_sample); + if (pc1 && pc1->numCoefs() > 3) + { + co1 = pc1->coefs_begin(); + co12 = pc1->coefs_end(); + } + else + { + double t1 = cv1->startparam(); + double t2 = cv1->endparam(); + double tdel = (t2 - t1)/(double)(num_sample+1); + double tpar; + int ka; + for (ka=0, tpar=t1+0.5*tdel; kapoint(tpar); + pts1.insert(pts1.end(), pt.begin(), pt.end()); + par1[ka] = tpar; + } + co1 = pts1.begin(); + co12 = pts1.end(); + } + + if (pc2 && pc2->numCoefs() > 3) + { + co1 = pc2->coefs_begin(); + co22 = pc2->coefs_end(); + } + else + { + double t1 = cv2->startparam(); + double t2 = cv2->endparam(); + double tdel = (t2 - t1)/(double)(num_sample+1); + double tpar; + int ka; + for (ka=0, tpar=t1+0.5*tdel; kapoint(tpar); + pts2.insert(pts2.end(), pt.begin(), pt.end()); + par2[ka] = tpar; + } + co2 = pts2.begin(); + co22 = pts2.end(); + } + + double td, tmin=1.0e8; + int minidx1=0, minidx2=0; + int ki, k1, k2; + for (k1=0; co1numCoefs() > 3) + { + std::vector::const_iterator st; + int kk = pc1->order(); + for (k1=minidx1+1, st=pc1->basis().begin(), seed1=0.0; + k1numCoefs() > 3) + { + std::vector::const_iterator st; + int kk = pc2->order(); + for (k1=minidx2+1, st=pc2->basis().begin(), seed2=0.0; + k1& intpar, + vector >& int_cvs) + //************************************************************************ + // + // Intersect a curve with a plane. Collect intersection parameters. + // + //*********************************************************************** +{ + shared_ptr spline_cv(cv->geometryCurve()); + + // Make sisl curves and call sisl. + SISLCurve *pc = Curve2SISL(*spline_cv, false); + + int kntrack = 0; + int trackflag = 0; // Do not make tracks. + SISLTrack **track =0; + int knpt=0, kncrv=0; + double *par=0; + int *pretop = 0; + SISLIntcurve **intcrvs = 0; + int stat = 0; + sh1372(pc, (double*)pos.begin(), (double*)axis.begin(), + radius, pos.dimension(), 0.0, epsge, + trackflag, &kntrack, &track, + &knpt, &par, &pretop, &kncrv, &intcrvs, &stat); + + ALWAYS_ERROR_IF(stat<0,"Error in intersection, code: " << stat); + + + // Remember intersections points. + int ki; + for (ki=0; kiipoint; + int_cvs.push_back(std::make_pair(intcrvs[ki]->epar1[0], + intcrvs[ki]->epar1[np-1])); + } + + if (kncrv > 0) + freeIntcrvlist(intcrvs, kncrv); + + if (par != 0) free(par); + if (pretop != 0) free(pretop); + if (pc != 0) freeCurve(pc); +} + } // namespace Go diff --git a/gotools-core/src/utils/BoundingBox.C b/gotools-core/src/utils/BoundingBox.C index cab0bf2ec..e82294500 100644 --- a/gotools-core/src/utils/BoundingBox.C +++ b/gotools-core/src/utils/BoundingBox.C @@ -170,8 +170,17 @@ void BoundingBox::addUnionWith(const Point& pt) void BoundingBox::addUnionWith(const BoundingBox& box) //=========================================================================== { - addUnionWith(box.low_); - addUnionWith(box.high_); + if (valid_) + { + addUnionWith(box.low_); + addUnionWith(box.high_); + } + else + { + low_ = box.low_; + high_ = box.high_; + valid_ = true; + } } //=========================================================================== diff --git a/gotools-core/src/utils/DirectionCone.C b/gotools-core/src/utils/DirectionCone.C index 34d4784f2..be915471c 100644 --- a/gotools-core/src/utils/DirectionCone.C +++ b/gotools-core/src/utils/DirectionCone.C @@ -202,6 +202,50 @@ void DirectionCone::addUnionWith(const Point& pt) } +//=========================================================================== +double DirectionCone::unionAngle(const Point& pt) const +//=========================================================================== +{ + ALWAYS_ERROR_IF(greater_than_pi_ < 0, + "Must initialize DirectionCone"); + ALWAYS_ERROR_IF(centre_.dimension() != pt.dimension(), + "Dimension mismatch"); + + Point t = pt; + double tlength = t.length(); + if (tlength == 0.0) { + MESSAGE("Ignoring vector of zero length"); + return angle_; + } + t /= tlength; + double theta = centre_ * t; + if (theta < -1.0) // This test and the following is needed for + theta = -1.0; // some reason... @jbt + if (theta > 1.0) + theta = 1.0; + theta = acos(theta); + if (theta <= angle_) + return angle_; + if (theta + angle_ >= M_PI + zero_tol_) { + greater_than_pi_ = 1; + return M_PI; + } + + // New code + double t1 = sin(0.5 * (theta + angle_)); + double t2 = sin(0.5 * (theta - angle_)); + Point centre = centre_ * t1 + t * t2; + if (centre.length() == 0.0) + // The cone spans 180 degrees. + // Happens for instance if first 2 vectors point in opposite directions. + return M_PI; + + double angle = 0.5 * (theta + angle_); + + return angle; +} + + //=========================================================================== void DirectionCone::addUnionWith(const DirectionCone& cone) //=========================================================================== @@ -237,7 +281,7 @@ void DirectionCone::addUnionWith(const DirectionCone& cone) void DirectionCone::check_angle() const //=========================================================================== { - ALWAYS_ERROR_IF(angle_ < 0.0, "Check failed -- negative angle."); + //ALWAYS_ERROR_IF(angle_ < 0.0, "Check failed -- negative angle."); if (angle_ <= M_PI + zero_tol_) greater_than_pi_ = 0; diff --git a/gotools-data b/gotools-data index 50cf659dc..77ebb337b 160000 --- a/gotools-data +++ b/gotools-data @@ -1 +1 @@ -Subproject commit 50cf659dc358cc77ad8a02456a6d67f90c7dc7cf +Subproject commit 77ebb337b5ed30fae965686fccb80d8cbdd54d10 diff --git a/sisl b/sisl index 885ade18a..66dbe489e 160000 --- a/sisl +++ b/sisl @@ -1 +1 @@ -Subproject commit 885ade18ac0e2ba852869da60de796522292a142 +Subproject commit 66dbe489ec58378f6d55d7976f4d3a2e50981485 diff --git a/trivariate/include/GoTools/trivariate/CurveOnVolume.h b/trivariate/include/GoTools/trivariate/CurveOnVolume.h index 13c920ffc..c5c7b6ec1 100644 --- a/trivariate/include/GoTools/trivariate/CurveOnVolume.h +++ b/trivariate/include/GoTools/trivariate/CurveOnVolume.h @@ -188,7 +188,8 @@ class GO_API CurveOnVolume : public ParamCurve // inherited from ParamCurve. NB: does not check whether the resulting ParamCurve // stays inside parameter domain (or that the space curve stays on surface). virtual void appendCurve(ParamCurve* cv, - int continuity, double& dist, bool reparam=true); + int continuity, double& dist, bool reparam=true, + double tol = 1.0e-4); /// Set the underlying surface to the one pointed to by the argument /// \param surface the pointer to the surface we will set as underlying for this diff --git a/trivariate/src/CurveOnVolume.C b/trivariate/src/CurveOnVolume.C index ceeb8b713..60e2b6b61 100644 --- a/trivariate/src/CurveOnVolume.C +++ b/trivariate/src/CurveOnVolume.C @@ -59,6 +59,7 @@ CurveOnVolume::CurveOnVolume() { } + //=========================================================================== CurveOnVolume::CurveOnVolume(shared_ptr vol, shared_ptr curve, @@ -717,7 +718,8 @@ void CurveOnVolume::appendCurve(ParamCurve* cv, bool reparam) //=========================================================================== void CurveOnVolume::appendCurve(ParamCurve* other_curve, - int continuity, double& dist, bool reparam) + int continuity, double& dist, bool reparam, + double tol) //=========================================================================== { CurveOnVolume *vol_cv = dynamic_cast(other_curve); From cd5e3abf3bef093a0fcc361d237241ff5f1562b2 Mon Sep 17 00:00:00 2001 From: vsk Date: Mon, 21 Oct 2024 10:11:06 +0200 Subject: [PATCH 02/24] updated refs --- gotools-data | 2 +- sisl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/gotools-data b/gotools-data index 77ebb337b..50cf659dc 160000 --- a/gotools-data +++ b/gotools-data @@ -1 +1 @@ -Subproject commit 77ebb337b5ed30fae965686fccb80d8cbdd54d10 +Subproject commit 50cf659dc358cc77ad8a02456a6d67f90c7dc7cf diff --git a/sisl b/sisl index 66dbe489e..885ade18a 160000 --- a/sisl +++ b/sisl @@ -1 +1 @@ -Subproject commit 66dbe489ec58378f6d55d7976f4d3a2e50981485 +Subproject commit 885ade18ac0e2ba852869da60de796522292a142 From 3ad40a91f528365186920ffdd3d02f6ffe107aa7 Mon Sep 17 00:00:00 2001 From: vsk Date: Mon, 21 Oct 2024 14:39:54 +0200 Subject: [PATCH 03/24] intersectWithTorus --- gotools-core/include/GoTools/geometry/BoundedUtils.h | 4 ++-- gotools-core/src/geometry/BoundedUtils.C | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/BoundedUtils.h b/gotools-core/include/GoTools/geometry/BoundedUtils.h index 868c57566..91c92fc5e 100644 --- a/gotools-core/include/GoTools/geometry/BoundedUtils.h +++ b/gotools-core/include/GoTools/geometry/BoundedUtils.h @@ -337,8 +337,8 @@ namespace BoundedUtils { std::vector > intersectWithTorus(shared_ptr& surf, - Point pnt, Point normal, double rad1, - double rad2, double geom_tol); + Point pnt, Point normal, double radius1, + double radius2, double geom_tol); /// Find the intersction curve(s) between two parametric surfaces. The surfaces /// are represented as spline surface if this was not the case initially. diff --git a/gotools-core/src/geometry/BoundedUtils.C b/gotools-core/src/geometry/BoundedUtils.C index 33c6b893c..469f582ab 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -2713,8 +2713,8 @@ BoundedUtils::intersectWithCone(shared_ptr& surf, //=========================================================================== vector > BoundedUtils::intersectWithTorus(shared_ptr& surf, - Point pnt, Point normal, double rad1, - double rad2, double geom_tol) + Point pnt, Point normal, double radius1, + double radius2, double geom_tol) //=========================================================================== { vector > curves; @@ -2741,7 +2741,7 @@ BoundedUtils::intersectWithTorus(shared_ptr& surf, SISLIntcurve** intcurves = 0; int stat; // Find the topology of the intersection - s1369(sislsf, pnt.begin(), normal.begin(), rad1, rad2, dim, epsco, + s1369(sislsf, pnt.begin(), normal.begin(), radius1, radius2, dim, epsco, geom_tol, &numintpt, &pointpar, &numintcr, &intcurves, &stat); // @@sbr Not sure this is the right solution. Maybe stat!=0 because of warning. ALWAYS_ERROR_IF(stat!=0, @@ -2754,7 +2754,7 @@ BoundedUtils::intersectWithTorus(shared_ptr& surf, // epsge = tol_.neighbour; for (int i = 0; i < numintcr; ++i) { // March out the intersection curves - s1318(sislsf,pnt.begin(), normal.begin(), rad1, rad2, dim, epsco, + s1318(sislsf,pnt.begin(), normal.begin(), radius1, radius2, dim, epsco, geom_tol, maxstep, intcurves[i], makecurv, graphic, &stat); SISLCurve* sc = intcurves[i]->pgeom; if (sc == 0) { From 422920bb71a4b7172e61500ea40067b58d796ee3 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Fri, 25 Oct 2024 14:39:14 +0200 Subject: [PATCH 04/24] Disabled some new code to allow regression test to pass. --- gotools-core/src/geometry/Cone.C | 4 ++++ gotools-core/src/geometry/Plane.C | 8 ++++++++ 2 files changed, 12 insertions(+) diff --git a/gotools-core/src/geometry/Cone.C b/gotools-core/src/geometry/Cone.C index 8714baa10..a371f1e7e 100644 --- a/gotools-core/src/geometry/Cone.C +++ b/gotools-core/src/geometry/Cone.C @@ -700,12 +700,14 @@ void Cone::closestPoint(const Point& pt, circle.setParamBounds(umin, umax); circle.closestPoint(pt, umin, umax, clo_u, clo_pt, clo_dist); clo_v = vmin; +#if 0 if (rad < 0) { clo_u += M_PI; if (clo_u > umax) clo_u -= (2.0*M_PI); } +#endif if (clo_dist < epsilon) { clo_u = domain_.umin() + (clo_u - parbound_.umin())*(domain_.umax()-domain_.umin())/(parbound_.umax()-parbound_.umin()); @@ -723,12 +725,14 @@ void Cone::closestPoint(const Point& pt, circle.setParamBounds(umin, umax); circle.closestPoint(pt, umin, umax, tmp_clo_u, tmp_clo_pt, tmp_clo_dist); tmp_clo_v = vmax; +#if 0 if (rad < 0) { tmp_clo_u += M_PI; if (tmp_clo_u > umax) tmp_clo_u -= (2.0*M_PI); } +#endif if (tmp_clo_dist < clo_dist) { clo_u = tmp_clo_u; clo_v = tmp_clo_v; diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 25cf0fc57..14a6f011d 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -795,15 +795,23 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, Point par2(parval2[0], parval2[1]); if (space_crv->instanceType() == Class_Line) { + return param_cv; + Point pos = (t2*par1 - t1*par2)/(t2 - t1); Point dir = (par2 - par1); dir.normalize(); param_cv = shared_ptr(new Line(pos, dir)); param_cv->setParamBounds(t1, t2); + + if (space_crv->isReversed()) { + param_cv->reverseParameterDirection(); + } } else { + return param_cv; + Point mid = space_crv->ParamCurve::point(0.5*(t1+t2)); double parval3[2]; double d3; From eef72442804194d414560a667f159abd9c3e0261 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Fri, 25 Oct 2024 15:22:44 +0200 Subject: [PATCH 05/24] Fixed issue with Line definition (now using correct function values for param bounds). --- gotools-core/src/geometry/Plane.C | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 14a6f011d..ac1c83d40 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -795,11 +795,9 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, Point par2(parval2[0], parval2[1]); if (space_crv->instanceType() == Class_Line) { - return param_cv; - + // We want L(t1) = par1 && L(t2) >= par2. Point pos = (t2*par1 - t1*par2)/(t2 - t1); - Point dir = (par2 - par1); - dir.normalize(); + Point dir = (par2 - par1)/(t2 - t1); param_cv = shared_ptr(new Line(pos, dir)); param_cv->setParamBounds(t1, t2); From e35d839dc721fe3d0f4592c70648ff4cb40661f2 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Mon, 28 Oct 2024 10:20:03 +0100 Subject: [PATCH 06/24] Elementary curve (circle) mods. --- gotools-core/src/geometry/Plane.C | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index ac1c83d40..8b7acb76b 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -808,7 +808,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, } else { - return param_cv; + //return param_cv; Point mid = space_crv->ParamCurve::point(0.5*(t1+t2)); double parval3[2]; @@ -926,6 +926,11 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, shared_ptr(new Circle(radius, centre, param_cv_axis, xvec, reversed)); param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); + + // if (space_crv->isReversed()) { + // param_cv->reverseParameterDirection(); + // } + } #ifdef DEBUG // TEST From bb077e665873aa51b8be4eff557f009ede30dd32 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Mon, 28 Oct 2024 18:08:11 +0100 Subject: [PATCH 07/24] Fixed typo. --- gotools-core/src/geometry/Plane.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 8b7acb76b..81a5aa7c2 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -795,7 +795,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, Point par2(parval2[0], parval2[1]); if (space_crv->instanceType() == Class_Line) { - // We want L(t1) = par1 && L(t2) >= par2. + // We want L(t1) = par1 && L(t2) = par2. Point pos = (t2*par1 - t1*par2)/(t2 - t1); Point dir = (par2 - par1)/(t2 - t1); From e4c1c1321fbca2a383d722f6ee9bf3ad329ba419 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Tue, 29 Oct 2024 11:55:50 +0100 Subject: [PATCH 08/24] Fixed issue with deriving par eps from space eps for plane. Now using the directional vectors. --- gotools-core/include/GoTools/geometry/Plane.h | 2 +- gotools-core/src/geometry/SurfaceTools.C | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/gotools-core/include/GoTools/geometry/Plane.h b/gotools-core/include/GoTools/geometry/Plane.h index a8282f555..2366774f3 100644 --- a/gotools-core/include/GoTools/geometry/Plane.h +++ b/gotools-core/include/GoTools/geometry/Plane.h @@ -182,7 +182,7 @@ class Plane : public ElementarySurface { return normal_; } /// Vectors in plane - void getSpanningVectors(Point& axis1, Point& axis2) + void getSpanningVectors(Point& axis1, Point& axis2) const { axis1 = vec1_; axis2 = vec2_; diff --git a/gotools-core/src/geometry/SurfaceTools.C b/gotools-core/src/geometry/SurfaceTools.C index 6203ab283..e4fcc9756 100644 --- a/gotools-core/src/geometry/SurfaceTools.C +++ b/gotools-core/src/geometry/SurfaceTools.C @@ -43,6 +43,7 @@ #include "GoTools/geometry/SurfaceOfRevolution.h" #include "GoTools/geometry/Cylinder.h" #include "GoTools/geometry/Cone.h" +#include "GoTools/geometry/Plane.h" #include "GoTools/geometry/SurfaceOfLinearExtrusion.h" using std::vector; @@ -817,6 +818,20 @@ Point SurfaceTools::getParEpsilon(const ParamSurface& sf, double epsgeo) std::swap(sf_epspar[0], sf_epspar[1]); } } + else if (sf.instanceType() == Class_Plane) + { + const Plane& plane = dynamic_cast(sf); + Point axis1, axis2; + plane.getSpanningVectors(axis1, axis2); + double length1 = axis1.length(); + double length2 = axis2.length(); + sf_epspar[0] = epsgeo*length1; + sf_epspar[1] = epsgeo*length2; + if (plane.isSwapped()) + { + std::swap(sf_epspar[0], sf_epspar[1]); + } + } else { // Set number of sampling points based on surface size From 87884e954c2ca8f17869a636abe6fb0925ff08c2 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Wed, 30 Oct 2024 16:23:09 +0100 Subject: [PATCH 09/24] Various debug and testing code. --- .../include/GoTools/geometry/Circle.h | 20 ++ .../include/GoTools/geometry/CurveLoop.h | 48 ++- gotools-core/src/geometry/Circle.C | 1 + gotools-core/src/geometry/CurveLoop.C | 1 + gotools-core/src/geometry/CurveOnSurface.C | 4 + gotools-core/src/geometry/Plane.C | 273 +++++++++--------- 6 files changed, 210 insertions(+), 137 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/Circle.h b/gotools-core/include/GoTools/geometry/Circle.h index d2a6706da..b17999d66 100644 --- a/gotools-core/include/GoTools/geometry/Circle.h +++ b/gotools-core/include/GoTools/geometry/Circle.h @@ -108,6 +108,15 @@ class Circle : public ElementaryCurve virtual double startparam() const; virtual double endparam() const; + double param_bounds_1() const + { + return parbound1_; + } + double param_bounds_2() const + { + return parbound2_; + } + //virtual void reverseParameterDirection(bool switchparam = false); /// The full circle is always parametrized on [0, 2*M_PI), it does @@ -169,6 +178,17 @@ class Circle : public ElementaryCurve return vec1_; } + Point getYAxis() const + { + return vec2_; + } + + void setYAxis(const Point& vec2) + { + vec2_ = vec2; + vec2_.normalize(); + } + double getRadius() const { return radius_; diff --git a/gotools-core/include/GoTools/geometry/CurveLoop.h b/gotools-core/include/GoTools/geometry/CurveLoop.h index 69fc9fb0d..f4ac5cc6d 100644 --- a/gotools-core/include/GoTools/geometry/CurveLoop.h +++ b/gotools-core/include/GoTools/geometry/CurveLoop.h @@ -45,6 +45,7 @@ #include #include "GoTools/utils/config.h" #include "GoTools/geometry/ParamCurve.h" +#include "GoTools/geometry/ElementaryCurve.h" namespace Go @@ -83,15 +84,50 @@ inline double computeLoopGap(const std::vector< PtrToCurveType >& curves) double maxdist = -1.0; double dist; for (i = 1; i < n; ++i) { - curves[i-1]->point(endp, curves[i-1]->endparam()); - curves[i]->point(startp, curves[i]->startparam()); + double endpar = curves[i-1]->endparam(); + shared_ptr elem_cv_prev = + dynamic_pointer_cast(curves[i-1]); + // if (elem_cv_prev && elem_cv_prev->isReversed()) { + // endpar = curves[i-1]->startparam(); + // } + curves[i-1]->point(endp, endpar); + + double startpar = curves[i]->startparam(); + shared_ptr elem_cv_curr = + dynamic_pointer_cast(curves[i]); + // if (elem_cv_curr && elem_cv_curr->isReversed()) { + // startpar = curves[i]->endparam(); + // } + curves[i]->point(startp, startpar); dist = endp.dist(startp); - if (dist > maxdist) maxdist = dist; + if (dist > maxdist) { + if (dist > 0.1) + std::cout << "i: " << i << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; + maxdist = dist; + } } - curves[n-1]->point(endp, curves[n-1]->endparam()); - curves[0]->point(startp, curves[0]->startparam()); + double endpar = curves[n-1]->endparam(); + shared_ptr elem_cv_prev = + dynamic_pointer_cast(curves[n-1]); + // if (elem_cv_prev && elem_cv_prev->isReversed()) { + // endpar = curves[n-1]->startparam(); + // } + curves[n-1]->point(endp, endpar); + + double startpar = curves[0]->startparam(); + shared_ptr elem_cv_curr = + dynamic_pointer_cast(curves[0]); + // if (elem_cv_curr && elem_cv_curr->isReversed()) { + // startpar = curves[0]->endparam(); + // } + curves[0]->point(startp, startpar); dist = endp.dist(startp); - if (dist > maxdist) maxdist = dist; + if (dist > maxdist) { + if (dist > 0.1) + std::cout << "i: " << 0 << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; + maxdist = dist; + } + std::cout << "maxdist: " << maxdist << std::endl; return maxdist; } diff --git a/gotools-core/src/geometry/Circle.C b/gotools-core/src/geometry/Circle.C index 42b89c32a..8bdafa512 100644 --- a/gotools-core/src/geometry/Circle.C +++ b/gotools-core/src/geometry/Circle.C @@ -254,6 +254,7 @@ Circle* Circle::clone() const isReversed_); circle->setParamBounds(parbound1_, parbound2_); circle->setParameterInterval(startparam_, endparam_); + circle->setYAxis(circle->getYAxis()); return circle; } diff --git a/gotools-core/src/geometry/CurveLoop.C b/gotools-core/src/geometry/CurveLoop.C index 70b8ce7dc..05a992945 100644 --- a/gotools-core/src/geometry/CurveLoop.C +++ b/gotools-core/src/geometry/CurveLoop.C @@ -424,6 +424,7 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) computeParLoopSpaceGap(cvs_on_sf) : -1.0; // The distance is lifted to space. if ((max_gap < space_epsilon_) && (par_loop_gap > 10.0*space_epsilon_)) { + std::cout << "max_gap: " << max_gap << ", space_epsilon_: " << space_epsilon_ << ", par_loop_gap: " << par_loop_gap << std::endl; return false; } diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 4a3c0f9a9..cdcc76ed1 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -1511,6 +1511,10 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, // The function returns a curve only if the configuration is simple pcurve_ = elem_sf->getElementaryParamCurve(elem_cv.get(), epspar, start_par_pt, end_par_pt); + if (pcurve_) { + bool debug_same_trace = sameTrace(epspar); + std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; + } } } diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 81a5aa7c2..c60e0bccf 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -745,13 +745,13 @@ bool Plane::isDegenerate(bool& b, bool& r, //=========================================================================== shared_ptr -Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, +Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, const Point* start_par_pt, const Point* end_par_pt) const //=========================================================================== { // Default is not simple elementary parameter curve exists shared_ptr param_cv; - + //std::cout << "epspar: " << epspar << std::endl; // Bookkeeping related to swapped parameters int ind1 = 0; int ind2 = 1; @@ -773,8 +773,10 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, { t1 = space_crv->startparam(); closed = ((Circle*)(space_crv))->isClosed(); - t2 = (closed) ? 0.5*(t1 + space_crv->endparam()) : - space_crv->endparam(); + //std::cout << "closed: " << closed << std::endl; + // t2 = (closed) ? 0.5*(t1 + space_crv->endparam()) : + // space_crv->endparam(); + t2 = space_crv->endparam(); idx = ind2; // 1 } else @@ -786,15 +788,14 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, Point close1, close2; Point pos1 = space_crv->ParamCurve::point(t1); Point pos2 = space_crv->ParamCurve::point(t2); - closestPoint(pos1, parval1[0], parval1[1], close1, d1, tol); - closestPoint(pos2, parval2[0], parval2[1], close2, d2, tol); - if (d1 > tol || d2 > tol) + closestPoint(pos1, parval1[0], parval1[1], close1, d1, epspar); + closestPoint(pos2, parval2[0], parval2[1], close2, d2, epspar); + if (d1 > epspar || d2 > epspar) return param_cv; Point par1(parval1[0], parval1[1]); Point par2(parval2[0], parval2[1]); - if (space_crv->instanceType() == Class_Line) - { + if (space_crv->instanceType() == Class_Line) { // We want L(t1) = par1 && L(t2) = par2. Point pos = (t2*par1 - t1*par2)/(t2 - t1); Point dir = (par2 - par1)/(t2 - t1); @@ -808,128 +809,138 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double tol, } else { - //return param_cv; - - Point mid = space_crv->ParamCurve::point(0.5*(t1+t2)); - double parval3[2]; - double d3; - Point close3; - closestPoint(mid, parval3[0], parval3[1], close3, d3, tol); - if (d3 > tol) - return param_cv; // Should not happen - - Point pmid(parval3[0], parval3[1]); - double alpha = t2 - t1; - double len = 0.5*par1.dist(par2); - double radius = len/sin(0.5*alpha); - Point vec = par1 + par2 - 2*pmid; - vec.normalize(); - Point centre = pmid + radius*vec; - Point param_cv_axis(0.0, 0.0); // A dimension two circle is not supported for every - // functionality - - Point centre_3d = ((Circle*)(space_crv))->getCentre(); - Point xvec_3d = ((Circle*)(space_crv))->getXAxis(); - Point vec1_3d = pos1 - centre_3d; - Point vec2_3d = pos2 - centre_3d; - double beta_3d = vec1_3d.angle(vec2_3d); - - Point vec1 = par1 - centre; - Point vec2 = par2 - centre; - double beta = vec1.angle(vec2); - - double gamma = 2*M_PI - t1; - double gamma2 = gamma; - if (gamma2 > 0.5*M_PI) - gamma2 -= 0.5*M_PI; - Point xvec(cos(gamma2)*(centre[0]-par1[0])-sin(gamma2)*(centre[1]-par1[1]), - sin(gamma2)*(centre[0]-par1[0])+cos(gamma2)*(centre[1]-par1[1])); - xvec.normalize(); - - Point xvec2(cos(t1)*(par1[0]-centre[0])+sin(t1)*(par1[1]-centre[1]), - -sin(t1)*(par1[0]-centre[0])+cos(t1)*(par1[1]-centre[1])); - xvec2.normalize(); - shared_ptr circ1(new Circle(radius, centre, param_cv_axis, - xvec, false)); - Point p1_1 = circ1->ParamCurve::point(t1); - Point p1_2 = circ1->ParamCurve::point(t2); - shared_ptr circ2(new Circle(radius, centre, param_cv_axis, - xvec2, false)); - Point p2_1 = circ2->ParamCurve::point(t1); - Point p2_2 = circ2->ParamCurve::point(t2); - Point p2_3 = circ2->ParamCurve::point(0.5*(t1+t2)); - - Point p3_1, p3_2; - point(p3_1, p2_1[0], p2_1[1]); - point(p3_2, p2_2[0], p2_2[1]); - double dd1 = pos1.dist(p3_1) + pos2.dist(p3_2); - if (dd1 > std::max(d1 + d2, tol)) - { - Point xvec3(cos(t2)*(par1[0]-centre[0])-sin(t2)*(par1[0]-centre[1]), - sin(t2)*(par1[1]-centre[0])+cos(t2)*(par1[1]-centre[1])); - xvec3.normalize(); - shared_ptr circ3(new Circle(radius, centre, param_cv_axis, - xvec3, true)); - Point p3_1 = circ3->ParamCurve::point(t1); - Point p3_2 = circ3->ParamCurve::point(t2); - Point p3_3 = circ3->ParamCurve::point(0.5*(t1+t2)); - - Point xvec4(xvec2[1],-xvec2[0]); - shared_ptr circ4(new Circle(radius, centre, param_cv_axis, - xvec4, false)); - Point p4_1 = circ4->ParamCurve::point(t1); - Point p4_2 = circ4->ParamCurve::point(t2); - Point p4_3 = circ4->ParamCurve::point(0.5*(t1+t2)); - - double t3 = space_crv->startparam() + space_crv->endparam() - t1; - Point xvec5(cos(t3)*(par1[0]-centre[0])-sin(t3)*(par1[0]-centre[1]), - sin(t3)*(par1[1]-centre[0])+cos(t3)*(par1[1]-centre[1])); - xvec5.normalize(); - shared_ptr circ5(new Circle(radius, centre, param_cv_axis, - xvec5, true)); - Point p5_1 = circ5->ParamCurve::point(t1); - Point p5_2 = circ5->ParamCurve::point(t2); - Point p5_3 = circ5->ParamCurve::point(0.5*(t1+t2)); - - - int stop_break_circ = 1; - } - - if (gamma2 < gamma) - { - std::swap(xvec[0], xvec[1]); - xvec[1] *= -1; - } - - double sgn1 = xvec_3d*vec2_3d; - double sgn2 = xvec*vec2; - bool reversed = false; - if (sgn1*sgn2 < 0.0) - { - // Opposite orientation of curve in geometry and parameter space. Redo xvec - // computation - gamma = 2*M_PI - t2; - gamma2 = gamma; - if (gamma2 > 0.5*M_PI) - gamma2 -= 0.5*M_PI; - xvec = Point(cos(gamma2)*(centre[0]-par1[0])-sin(gamma2)*(centre[1]-par1[1]), - sin(gamma2)*(centre[0]-par1[0])+cos(gamma2)*(centre[1]-par1[1])); - xvec.normalize(); - if (gamma2 < gamma) - { - std::swap(xvec[0], xvec[1]); - xvec[1] *= -1; - } - reversed = true; - } - param_cv = (dd1 <= tol) ? circ2 : - shared_ptr(new Circle(radius, centre, param_cv_axis, - xvec, reversed)); - param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); - - // if (space_crv->isReversed()) { - // param_cv->reverseParameterDirection(); - // } + //return param_cv; + + // For the circle we must check: + // 1) Space circle normal vs plane normal (may be flipped). + // 2) Plane swapped_. + // 3) Space curve: Is reversed. + + double ang_tol = 1e-04; + Point space_cv_normal = ((Circle*)(space_crv))->getNormal(); + double normal_ang = normal_.angle(space_cv_normal); + if ((normal_ang > ang_tol) && (normal_ang < M_PI - ang_tol)) { + return param_cv; + } + bool normal_flipped = (fabs(M_PI - normal_ang) < ang_tol); + //std::cout << "reversed: " << space_crv->isReversed() << std::endl; + + Point centre_2d(2); + Point centre_3d = ((Circle*)(space_crv))->getCentre(); + double d4; + Point close4; + closestPoint(centre_3d, centre_2d[0], centre_2d[1], close4, d4, epspar); + if (d4 > epspar) + return param_cv; + + double rad_par = centre_2d.dist(par1); + + //std::cout << "rad_par: " << rad_par << ", radius: " << radius << std::endl; + Point x_axis = ((Circle*)(space_crv))->getXAxis(); + Point x_axis_end = centre_3d + x_axis; + double d5; + Point close5; + Point x_axis_par_end(2); + closestPoint(x_axis_end, x_axis_par_end[0], x_axis_par_end[1], close5, d5, epspar); + if (d5 > epspar) + return param_cv; + Point x_axis_par = x_axis_par_end - centre_2d; + + Point y_axis = ((Circle*)(space_crv))->getYAxis(); + Point y_axis_end = centre_3d + y_axis; + double d6; + Point close6; + Point y_axis_par_end(2); + closestPoint(y_axis_end, y_axis_par_end[0], y_axis_par_end[1], close6, d6, epspar); + if (d6 > epspar) + return param_cv; + Point y_axis_par_proj = y_axis_par_end - centre_2d; + y_axis_par_proj.normalize(); + + Point y_axis_par(-x_axis_par[1], x_axis_par[0]); + //= ((Circle*)(param_cv.get()))->getYAxis(); + double ang_y_axis_par = y_axis_par.angle(y_axis_par_proj); + bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); + + bool reversed = (space_crv->isReversed()); + //double sign = (reversed) ? -1.0 : 1.0; + + Point param_cv_axis(0.0, 0.0); + param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, x_axis_par, reversed)); + + param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); + //param_cv->setParameterInterval(space_crv->startparam(), space_crv->endparam()); + + if (y_axis_reversed)// && (sign > 0)) + { + std::cout << "Reversing y axis." << std::endl; + //((Circle*)(param_cv.get()))->setYAxis(sign*y_axis_par_proj); + ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); + std::cout << "Done reversing y axis." << std::endl; + } + + //std::cout << "Setting the Circle!" << std::endl; + // We calculate the dist from the lifted point to the space cv point. + Point par_start = param_cv->point(param_cv->startparam()); + double dist_par1 = par1.dist(par_start); + Point par_end = param_cv->point(param_cv->endparam()); + double dist_par2 = par2.dist(par_end); + + // Calculate intermediate parameter. + double tpar = (1.0*param_cv->startparam() + 2.0*param_cv->endparam())/3.0; + Point tpar_3d = space_crv->point(tpar); + Point tpar_2d = param_cv->point(tpar); + Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); + double tpar_dist = tpar_3d.dist(tpar_sf_3d); + + if (start_par_pt != NULL) { + double d0 = start_par_pt->dist(par_start); + if (d0 > epspar) + std::cout << "d0: " << d0 << std::endl; + } + if (end_par_pt != NULL) { + double d1 = end_par_pt->dist(par_end); + if (d1 > epspar) + std::cout << "d1: " << d1 << std::endl; + } + + std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; + std::cout << "tpar_dist: " << tpar_dist << std::endl; + if (std::max(dist_par1, dist_par2) > epspar) + { + + double parb1 = ((Circle*)(space_crv))->param_bounds_1(); + double parb2 = ((Circle*)(space_crv))->param_bounds_2(); + std::cout << "\n\nparb1: " << parb1 << ", parb2: " << parb2 << std::endl; + std::cout << "Something wrong with the calculated circle! Mis-match at end point." << std::endl; + std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; + std::cout << "par1: " << par1 << ", par2: " << par2 << std::endl; + std::cout << "par_start: " << par_start << ", par_end: " << par_end << std::endl; + std::cout << "normal_: " << normal_ << std::endl; + std::cout << "normal_flipped: " << normal_flipped << std::endl; + std::cout << "param_cv->startparam(): " << param_cv->startparam() << std::endl; + std::cout << "param_cv->endparam(): " << param_cv->endparam() << std::endl; + std::cout << "swapped: " << isSwapped() << std::endl; + if (start_par_pt != NULL) { + std::cout << "start_par_pt: " << *start_par_pt << std::endl; + } + if (end_par_pt != NULL) { + std::cout << "end_par_pt: " << *end_par_pt << std::endl; + } + // std::cout << "xvec: " << xvec << std::endl; + // std::cout << "xvec2: " << xvec2 << std::endl; + std::cout << "x_axis_par: " << x_axis_par << std::endl; + std::cout << "space_cv reversed: " << reversed << std::endl; + std::cout << "space_cv normal: " << ((Circle*)(space_crv))->getNormal() << std::endl; + std::cout << "normal_ang: " << normal_ang << std::endl; + std::cout << "y_axis_par: " << y_axis_par << std::endl; + std::cout << "y_axis_par_proj: " << y_axis_par_proj << std::endl; + if (ang_y_axis_par > 0.5*M_PI) { + std::cout << "Flipped: ang_y_axis_par: " << ang_y_axis_par << std::endl; + } else { + std::cout << "Not flipped: ang_y_axis_par: " << ang_y_axis_par << std::endl; + } + } } #ifdef DEBUG From b74ccb5a4ea1680940e0796637d46335011c475a Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Wed, 30 Oct 2024 17:15:16 +0100 Subject: [PATCH 10/24] Removed cout. --- gotools-core/include/GoTools/geometry/CurveLoop.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gotools-core/include/GoTools/geometry/CurveLoop.h b/gotools-core/include/GoTools/geometry/CurveLoop.h index f4ac5cc6d..784561d16 100644 --- a/gotools-core/include/GoTools/geometry/CurveLoop.h +++ b/gotools-core/include/GoTools/geometry/CurveLoop.h @@ -127,7 +127,7 @@ inline double computeLoopGap(const std::vector< PtrToCurveType >& curves) std::cout << "i: " << 0 << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; maxdist = dist; } - std::cout << "maxdist: " << maxdist << std::endl; + //std::cout << "maxdist: " << maxdist << std::endl; return maxdist; } From bef5bbb7bb551ddbf7bbe1a659650c2ffc6f4aa2 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Wed, 30 Oct 2024 20:18:47 +0100 Subject: [PATCH 11/24] Fixed tolerance issue for debug write outer boundary loop. --- gotools-core/include/GoTools/geometry/SplineDebugUtils.h | 1 + gotools-core/src/geometry/BoundedUtils.C | 3 ++- gotools-core/src/geometry/SplineDebugUtils.C | 3 ++- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h index 21e9945ba..a7fa2733c 100644 --- a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h +++ b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h @@ -85,6 +85,7 @@ namespace SplineDebugUtils std::ostream& os, double z = 0.0); void GO_API writeOuterBoundaryLoop(ParamSurface& sf, + double epsgeo, std::ostream& os); void GO_API writeBoundary(BoundedSurface& sf, std::ostream& os); diff --git a/gotools-core/src/geometry/BoundedUtils.C b/gotools-core/src/geometry/BoundedUtils.C index 469f582ab..bc3fc45f0 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -5372,7 +5372,8 @@ bool BoundedUtils::createMissingParCvs(CurveLoop& bd_loop, bool loop_is_ccw) #ifndef NDEBUG { std::ofstream debug3("tmp/undersf_outer_loop.g2"); - Go::SplineDebugUtils::writeOuterBoundaryLoop(*under_sf, debug3); + Go::SplineDebugUtils::writeOuterBoundaryLoop(*under_sf, bd_loop.getSpaceEpsilon(), + debug3); double debug_val = 0.0; } #endif // NDEBUG diff --git a/gotools-core/src/geometry/SplineDebugUtils.C b/gotools-core/src/geometry/SplineDebugUtils.C index d03cc6d3b..a171b222b 100644 --- a/gotools-core/src/geometry/SplineDebugUtils.C +++ b/gotools-core/src/geometry/SplineDebugUtils.C @@ -299,11 +299,12 @@ void SplineDebugUtils::writeBoundary(BoundedSurface& bd_sf, //=========================================================================== void SplineDebugUtils::writeOuterBoundaryLoop(ParamSurface& sf, + double epsgeo, std::ostream& os) //=========================================================================== { // We also write the boundary loops of the underlying surface. - CurveLoop outer_loop = sf.outerBoundaryLoop(); + CurveLoop outer_loop = sf.outerBoundaryLoop(epsgeo); for (size_t ki = 0; ki < outer_loop.size(); ++ki) { outer_loop[ki]->writeStandardHeader(os); From 78752acb2ad20c3ccce89715662fbfe7e6d3dc2b Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Thu, 31 Oct 2024 14:59:22 +0100 Subject: [PATCH 12/24] Various debug output and experimental code. --- gotools-core/src/geometry/CurveLoop.C | 7 +++++++ gotools-core/src/geometry/CurveOnSurface.C | 8 +++++++- gotools-core/src/geometry/Plane.C | 18 ++++++++++-------- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/gotools-core/src/geometry/CurveLoop.C b/gotools-core/src/geometry/CurveLoop.C index 05a992945..75bf4a68b 100644 --- a/gotools-core/src/geometry/CurveLoop.C +++ b/gotools-core/src/geometry/CurveLoop.C @@ -549,6 +549,12 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) } double maxgap_par = computeLoopGap(par_cvs); double maxgap_space = computeLoopGap(space_cvs); + if (maxgap_space > space_epsilon_) { + std::cout << "WARN: maxgap_space: " << maxgap_space << std::endl; + } + if (maxgap_par > space_epsilon_) { + std::cout << "WARN: maxgap_par: " << maxgap_par << std::endl; + } if ((maxgap_space >= 0.0) && (maxgap_space < space_epsilon_)) { for (int k = 0; k < (int)curves.size(); ++k) @@ -585,6 +591,7 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) max_gap2 = computeLoopGap(curves); if (max_gap2 < max_gap) { + std::cout << "WARN: Replacing the loop curves." << std::endl; curves_ = curves; max_gap = max_gap2; } diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index cdcc76ed1..eeacf3b4d 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -55,6 +55,7 @@ #include "GoTools/creators/HermiteAppS.h" #include "GoTools/creators/CurveCreators.h" #include "GoTools/creators/CoonsPatchGen.h" +#include "GoTools/utils/Logger.h" #include #include @@ -1513,7 +1514,12 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, start_par_pt, end_par_pt); if (pcurve_) { bool debug_same_trace = sameTrace(epspar); - std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; + if (!debug_same_trace) { + LOG_DEBUG("debug_same_trace: " + debug_same_trace); + std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl;; + } else { + std::cout << "DEBUG: OK: debug_same_trace: " << debug_same_trace << std::endl; + } } } } diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index c60e0bccf..5b8921ace 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -803,9 +803,9 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, param_cv = shared_ptr(new Line(pos, dir)); param_cv->setParamBounds(t1, t2); - if (space_crv->isReversed()) { - param_cv->reverseParameterDirection(); - } + // if (space_crv->isReversed()) { + // param_cv->reverseParameterDirection(); + // } } else { @@ -862,7 +862,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, double ang_y_axis_par = y_axis_par.angle(y_axis_par_proj); bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); - bool reversed = (space_crv->isReversed()); + bool reversed = false;//(space_crv->isReversed()); //double sign = (reversed) ? -1.0 : 1.0; Point param_cv_axis(0.0, 0.0); @@ -873,10 +873,10 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, if (y_axis_reversed)// && (sign > 0)) { - std::cout << "Reversing y axis." << std::endl; + //std::cout << "Reversing y axis." << std::endl; //((Circle*)(param_cv.get()))->setYAxis(sign*y_axis_par_proj); ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); - std::cout << "Done reversing y axis." << std::endl; + //std::cout << "Done reversing y axis." << std::endl; } //std::cout << "Setting the Circle!" << std::endl; @@ -904,8 +904,10 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, std::cout << "d1: " << d1 << std::endl; } - std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; - std::cout << "tpar_dist: " << tpar_dist << std::endl; + if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { + std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; + std::cout << "tpar_dist: " << tpar_dist << std::endl; + } if (std::max(dist_par1, dist_par2) > epspar) { From 4f1d816502a162627a6d6493a8acce22e9aab70c Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Fri, 1 Nov 2024 20:42:37 +0100 Subject: [PATCH 13/24] Finally fixed issue with 2d circle and direction. --- .../include/GoTools/geometry/Circle.h | 6 +++ .../include/GoTools/geometry/GeometryTools.h | 4 +- gotools-core/src/geometry/BoundedUtils.C | 8 +++- gotools-core/src/geometry/CurveOnSurface.C | 9 ++-- gotools-core/src/geometry/Plane.C | 47 ++++++++++++++----- gotools-core/src/geometry/SplineDebugUtils.C | 9 +++- 6 files changed, 63 insertions(+), 20 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/Circle.h b/gotools-core/include/GoTools/geometry/Circle.h index b17999d66..16f2ef049 100644 --- a/gotools-core/include/GoTools/geometry/Circle.h +++ b/gotools-core/include/GoTools/geometry/Circle.h @@ -183,6 +183,12 @@ class Circle : public ElementaryCurve return vec2_; } + void setXAxis(const Point& vec1) + { + vec1_ = vec1; + setSpanningVectors(); + } + void setYAxis(const Point& vec2) { vec2_ = vec2; diff --git a/gotools-core/include/GoTools/geometry/GeometryTools.h b/gotools-core/include/GoTools/geometry/GeometryTools.h index 2592bf2b7..fd29290bb 100644 --- a/gotools-core/include/GoTools/geometry/GeometryTools.h +++ b/gotools-core/include/GoTools/geometry/GeometryTools.h @@ -405,7 +405,7 @@ namespace GeometryTools /// \param lc reference to the LineCloud that is to be rotated void GO_API rotateLineCloud(Point rot_axis, double alpha, LineCloud& lc); - /// Rotate the given 3D point a certain angle around a certain axis. + /// Rotate the given 2D or 3D point a certain angle around a certain axis. /// \param rot_axis the axis of rotation. It does not have to be normalized, but /// must of course be nonzero. /// \param alpha the angle of rotation, given in radians @@ -413,7 +413,7 @@ namespace GeometryTools /// the point are stored. These will be overwritten with the /// rotated coordinates. void GO_API rotatePoint(Point rot_axis, double alpha, double* space_pt); - /// Rotate the given 3D point a certain angle around a certain axis. + /// Rotate the given 2D or 3D point a certain angle around a certain axis. /// \param rot_axis the axis of rotation. It does not have to be normalized, but /// must of course be nonzero. /// \param alpha the angle of rotation, given in radians diff --git a/gotools-core/src/geometry/BoundedUtils.C b/gotools-core/src/geometry/BoundedUtils.C index bc3fc45f0..74eb3853a 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -39,7 +39,7 @@ #define DEBUG1 -//#define SBR_DBG +#define SBR_DBG #include "GoTools/geometry/BoundedUtils.h" #include "GoTools/utils/Logger.h" @@ -5002,6 +5002,12 @@ void BoundedUtils::fixInvalidBoundedSurface(shared_ptr& bd_sf, if ((pos_state%8 > 1) || (pos_state%16 > 1)) { // Loops not closed or in wrong order/orientation. +#ifdef SBR_DBG + { + std::ofstream outfile_failures("tmp/bd_sf_failures.g2"); + SplineDebugUtils::writeTrimmedInfo(*bd_sf, outfile_failures, 0.0); + } +#endif double max_gap = -1.0; bool success = bd_sf->fixInvalidSurface(max_gap); if (!success) diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index eeacf3b4d..9d5138faf 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -37,7 +37,7 @@ * written agreement between you and SINTEF ICT. */ -//#define SBR_DBG +#define SBR_DBG //#define DEBUG #include "GoTools/geometry/CurveOnSurface.h" @@ -1516,10 +1516,11 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, bool debug_same_trace = sameTrace(epspar); if (!debug_same_trace) { LOG_DEBUG("debug_same_trace: " + debug_same_trace); - std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl;; - } else { - std::cout << "DEBUG: OK: debug_same_trace: " << debug_same_trace << std::endl; + std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; } + // else { + // std::cout << "DEBUG: OK: debug_same_trace: " << debug_same_trace << std::endl; + // } } } } diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 5b8921ace..dcfdcf004 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -41,6 +41,7 @@ #include "GoTools/geometry/Line.h" #include "GoTools/geometry/Circle.h" #include "GoTools/geometry/SplineSurface.h" +#include "GoTools/geometry/GeometryTools.h" #include #include @@ -802,10 +803,16 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, param_cv = shared_ptr(new Line(pos, dir)); param_cv->setParamBounds(t1, t2); - - // if (space_crv->isReversed()) { - // param_cv->reverseParameterDirection(); - // } + // We do not reverse the parameter direction of the param cv even if space curve is reversed. Already handled by + // the evaluation. + + Point par_cv_1 = param_cv->point(t1); + double d1 = par1.dist(par_cv_1); + Point par_cv_2 = param_cv->point(t2); + double d2 = par2.dist(par_cv_2); + if (std::max(d1, d2) > epspar) { + std::cout << "DEBUG: Line: d1: " << d1 << ", d2: " << d2 << std::endl; + } } else { @@ -857,12 +864,12 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point y_axis_par_proj = y_axis_par_end - centre_2d; y_axis_par_proj.normalize(); - Point y_axis_par(-x_axis_par[1], x_axis_par[0]); + Point y_axis_par(-x_axis_par[1], x_axis_par[0]); // This is the hardcoded assignment of vec2_ in Circle. //= ((Circle*)(param_cv.get()))->getYAxis(); double ang_y_axis_par = y_axis_par.angle(y_axis_par_proj); bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); - bool reversed = false;//(space_crv->isReversed()); + bool reversed = (space_crv->isReversed()); //double sign = (reversed) ? -1.0 : 1.0; Point param_cv_axis(0.0, 0.0); @@ -875,11 +882,27 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, { //std::cout << "Reversing y axis." << std::endl; //((Circle*)(param_cv.get()))->setYAxis(sign*y_axis_par_proj); - ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); //std::cout << "Done reversing y axis." << std::endl; + +#if 0 + std::cout << "WARNING! This approach will fail as the vec2_ i2 forgotten when cloning or writing to file!" + << std::endl; + ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); +#else + //std::cout << "WARNING! Trying to reverse parameter direction to fix vec2_ issue!" << std::endl; + param_cv->reverseParameterDirection(); + double vec1_rot_ang = 2*M_PI - t2 - t1; // Rotate ccw. + Point new_x_axis_par = x_axis_par; + Point rot_axis_not_used(0.0, 0.0); // Rot axis is not used for 2d case. Always ccw in the plane. + GeometryTools::rotatePoint(rot_axis_not_used, vec1_rot_ang, &new_x_axis_par[0]); + param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, new_x_axis_par, !reversed)); + param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); + //((Circle*)param_cv.get() )->setXAxis(y_axis_par_proj); +#endif + + } - //std::cout << "Setting the Circle!" << std::endl; // We calculate the dist from the lifted point to the space cv point. Point par_start = param_cv->point(param_cv->startparam()); double dist_par1 = par1.dist(par_start); @@ -892,6 +915,10 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point tpar_2d = param_cv->point(tpar); Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); double tpar_dist = tpar_3d.dist(tpar_sf_3d); + if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { + std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; + std::cout << "tpar_dist: " << tpar_dist << std::endl; + } if (start_par_pt != NULL) { double d0 = start_par_pt->dist(par_start); @@ -904,10 +931,6 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, std::cout << "d1: " << d1 << std::endl; } - if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { - std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; - std::cout << "tpar_dist: " << tpar_dist << std::endl; - } if (std::max(dist_par1, dist_par2) > epspar) { diff --git a/gotools-core/src/geometry/SplineDebugUtils.C b/gotools-core/src/geometry/SplineDebugUtils.C index a171b222b..c7ad8e251 100644 --- a/gotools-core/src/geometry/SplineDebugUtils.C +++ b/gotools-core/src/geometry/SplineDebugUtils.C @@ -202,9 +202,16 @@ void SplineDebugUtils::writeTrimmedInfo(BoundedSurface& bd_sf, writeSpaceParamCurve(*line_cv, os); } + else if (par_cv->instanceType() == Class_Circle) { + shared_ptr circle_cv = + dynamic_pointer_cast + (par_cv); + writeSpaceParamCurve(*circle_cv, + os); + } else { - MESSAGE("Curve type not supported!"); + MESSAGE("Curve type not supported! par_cv->instanceType(): " + par_cv->instanceType()); } } shared_ptr space_cv = From 9e202ffcd0b8a562c83ba58c4308f82e97b00751 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 08:42:06 +0100 Subject: [PATCH 14/24] Removed debug file write. --- gotools-core/src/geometry/BoundedSurface.C | 2 +- gotools-core/src/geometry/BoundedUtils.C | 6 +++--- gotools-core/src/geometry/CurveOnSurface.C | 2 +- gotools-core/src/geometry/Plane.C | 6 +++--- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/gotools-core/src/geometry/BoundedSurface.C b/gotools-core/src/geometry/BoundedSurface.C index b086b40d7..b1a3732c9 100644 --- a/gotools-core/src/geometry/BoundedSurface.C +++ b/gotools-core/src/geometry/BoundedSurface.C @@ -73,7 +73,7 @@ using std::endl; //#define CHECK_PARAM_LOOP_ORIENTATION #ifndef NDEBUG -#define SBR_DBG +//#define SBR_DBG #include "GoTools/geometry/SplineDebugUtils.h" #endif diff --git a/gotools-core/src/geometry/BoundedUtils.C b/gotools-core/src/geometry/BoundedUtils.C index 74eb3853a..064f92ad0 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -37,9 +37,9 @@ * written agreement between you and SINTEF ICT. */ -#define DEBUG1 +//#define DEBUG1 -#define SBR_DBG +//#define SBR_DBG #include "GoTools/geometry/BoundedUtils.h" #include "GoTools/utils/Logger.h" @@ -5012,7 +5012,7 @@ void BoundedUtils::fixInvalidBoundedSurface(shared_ptr& bd_sf, bool success = bd_sf->fixInvalidSurface(max_gap); if (!success) { - LOG_INFO("max_gap = " + std::to_string(max_gap)); + LOG_INFO("Loop not closed or in wrong order/orientation, failed to fix. max_gap = " + std::to_string(max_gap)); } } diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 9d5138faf..3962fb505 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -37,7 +37,7 @@ * written agreement between you and SINTEF ICT. */ -#define SBR_DBG +//#define SBR_DBG //#define DEBUG #include "GoTools/geometry/CurveOnSurface.h" diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index dcfdcf004..bd4ba82d8 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -52,7 +52,7 @@ using std::numeric_limits; using std::streamsize; using std::swap; -#define DEBUG +//#define DEBUG namespace Go { @@ -916,8 +916,8 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); double tpar_dist = tpar_3d.dist(tpar_sf_3d); if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { - std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; - std::cout << "tpar_dist: " << tpar_dist << std::endl; + std::cout << "ERROR: dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << "tpar_dist: " << + tpar_dist << std::endl; } if (start_par_pt != NULL) { From 4f542b439ce62ea4d99ee346316b4e6fd59afde7 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 09:48:13 +0100 Subject: [PATCH 15/24] Increased max log size to 100 MB (before rotation). --- gotools-core/include/GoTools/utils/Logger.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gotools-core/include/GoTools/utils/Logger.h b/gotools-core/include/GoTools/utils/Logger.h index 711da3123..4435dd446 100644 --- a/gotools-core/include/GoTools/utils/Logger.h +++ b/gotools-core/include/GoTools/utils/Logger.h @@ -67,7 +67,7 @@ class Logger { try { std::cout << "Setting logfile: " << logfile_name << std::endl; // Open the file in truncate mode to clear its contents - file_logger = spdlog::rotating_logger_mt("file_logger", logfile_name, 1048576 * 5, 3); // 5 MB size, 3 files + file_logger = spdlog::rotating_logger_mt("file_logger", logfile_name, 1048576 * 100, 3); // 100 MB size, 3 files std::ofstream(logfile_name, std::ios::trunc).close(); // Truncate the file spdlog::set_default_logger(file_logger); spdlog::set_pattern("[%Y-%m-%d %H:%M:%S.%e] [%^%l%$] %v"); @@ -173,4 +173,4 @@ class Logger { #endif -} // namespace Go \ No newline at end of file +} // namespace Go From 98cbeb748b5e81361a018a8a75914471c78e9caa Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 09:51:46 +0100 Subject: [PATCH 16/24] Added logging of fixes. --- gotools-core/src/geometry/CurveLoop.C | 4 +++- gotools-core/src/geometry/Plane.C | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/gotools-core/src/geometry/CurveLoop.C b/gotools-core/src/geometry/CurveLoop.C index 75bf4a68b..b4aaa36cc 100644 --- a/gotools-core/src/geometry/CurveLoop.C +++ b/gotools-core/src/geometry/CurveLoop.C @@ -42,6 +42,7 @@ #include "GoTools/geometry/CurveOnSurface.h" #include "GoTools/geometry/orientCurves.h" #include "GoTools/geometry/SplineCurve.h" +#include "GoTools/utils/Logger.h" #include //#define DEBUG @@ -592,6 +593,7 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) if (max_gap2 < max_gap) { std::cout << "WARN: Replacing the loop curves." << std::endl; + LOG_WARN("Replacing the loop curves."); curves_ = curves; max_gap = max_gap2; } @@ -599,7 +601,7 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) #if 1 if (max_gap > space_epsilon_) { - MESSAGE("WARNING: Loop max_gap = " << max_gap << ", space_epsilon_ = " << space_epsilon_ << "!"); + LOG_WARN("Loop max_gap = " + std::to_string(max_gap) + ", space_epsilon_ = " + std::to_string(space_epsilon_) + "!"); } return true; diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index bd4ba82d8..7baeb2422 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -42,6 +42,7 @@ #include "GoTools/geometry/Circle.h" #include "GoTools/geometry/SplineSurface.h" #include "GoTools/geometry/GeometryTools.h" +#include "GoTools/utils/Logger.h" #include #include @@ -890,6 +891,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); #else //std::cout << "WARNING! Trying to reverse parameter direction to fix vec2_ issue!" << std::endl; + LOG_WARN("Fixing 2d circle with flipped normal."); param_cv->reverseParameterDirection(); double vec1_rot_ang = 2*M_PI - t2 - t1; // Rotate ccw. Point new_x_axis_par = x_axis_par; From 38feca293cd50fdd38ab760616433d476ecbe550 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 10:25:07 +0100 Subject: [PATCH 17/24] Removed debug cout, moved some of them to logger. --- .../include/GoTools/geometry/CurveLoop.h | 8 +-- gotools-core/src/geometry/CurveLoop.C | 14 +++--- gotools-core/src/geometry/Plane.C | 50 +------------------ 3 files changed, 12 insertions(+), 60 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/CurveLoop.h b/gotools-core/include/GoTools/geometry/CurveLoop.h index 784561d16..7a40fbef6 100644 --- a/gotools-core/include/GoTools/geometry/CurveLoop.h +++ b/gotools-core/include/GoTools/geometry/CurveLoop.h @@ -101,8 +101,8 @@ inline double computeLoopGap(const std::vector< PtrToCurveType >& curves) curves[i]->point(startp, startpar); dist = endp.dist(startp); if (dist > maxdist) { - if (dist > 0.1) - std::cout << "i: " << i << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; + // if (dist > 0.1) + // std::cout << "i: " << i << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; maxdist = dist; } } @@ -123,8 +123,8 @@ inline double computeLoopGap(const std::vector< PtrToCurveType >& curves) curves[0]->point(startp, startpar); dist = endp.dist(startp); if (dist > maxdist) { - if (dist > 0.1) - std::cout << "i: " << 0 << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; + // if (dist > 0.1) + // std::cout << "i: " << 0 << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; maxdist = dist; } //std::cout << "maxdist: " << maxdist << std::endl; diff --git a/gotools-core/src/geometry/CurveLoop.C b/gotools-core/src/geometry/CurveLoop.C index b4aaa36cc..29324d2a5 100644 --- a/gotools-core/src/geometry/CurveLoop.C +++ b/gotools-core/src/geometry/CurveLoop.C @@ -425,7 +425,8 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) computeParLoopSpaceGap(cvs_on_sf) : -1.0; // The distance is lifted to space. if ((max_gap < space_epsilon_) && (par_loop_gap > 10.0*space_epsilon_)) { - std::cout << "max_gap: " << max_gap << ", space_epsilon_: " << space_epsilon_ << ", par_loop_gap: " << par_loop_gap << std::endl; + LOG_WARN("max_gap: " + std::to_string(max_gap) + ", space_epsilon_: " + std::to_string(space_epsilon_) + + ", par_loop_gap: " + std::to_string(par_loop_gap)); return false; } @@ -551,10 +552,10 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) double maxgap_par = computeLoopGap(par_cvs); double maxgap_space = computeLoopGap(space_cvs); if (maxgap_space > space_epsilon_) { - std::cout << "WARN: maxgap_space: " << maxgap_space << std::endl; + LOG_WARN("maxgap_space: " + std::to_string(maxgap_space)); } if (maxgap_par > space_epsilon_) { - std::cout << "WARN: maxgap_par: " << maxgap_par << std::endl; + LOG_WARN("maxgap_par: " + std::to_string(maxgap_par)); } if ((maxgap_space >= 0.0) && (maxgap_space < space_epsilon_)) { @@ -591,12 +592,11 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) // If we got closer we replace the loop curves. max_gap2 = computeLoopGap(curves); if (max_gap2 < max_gap) - { - std::cout << "WARN: Replacing the loop curves." << std::endl; - LOG_WARN("Replacing the loop curves."); + { + LOG_WARN("Replacing the loop curves."); curves_ = curves; max_gap = max_gap2; - } + } #if 1 if (max_gap > space_epsilon_) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 7baeb2422..fe9810890 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -911,7 +911,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point par_end = param_cv->point(param_cv->endparam()); double dist_par2 = par2.dist(par_end); - // Calculate intermediate parameter. + // Testing dist for intermediate parameter. double tpar = (1.0*param_cv->startparam() + 2.0*param_cv->endparam())/3.0; Point tpar_3d = space_crv->point(tpar); Point tpar_2d = param_cv->point(tpar); @@ -921,54 +921,6 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, std::cout << "ERROR: dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << "tpar_dist: " << tpar_dist << std::endl; } - - if (start_par_pt != NULL) { - double d0 = start_par_pt->dist(par_start); - if (d0 > epspar) - std::cout << "d0: " << d0 << std::endl; - } - if (end_par_pt != NULL) { - double d1 = end_par_pt->dist(par_end); - if (d1 > epspar) - std::cout << "d1: " << d1 << std::endl; - } - - if (std::max(dist_par1, dist_par2) > epspar) - { - - double parb1 = ((Circle*)(space_crv))->param_bounds_1(); - double parb2 = ((Circle*)(space_crv))->param_bounds_2(); - std::cout << "\n\nparb1: " << parb1 << ", parb2: " << parb2 << std::endl; - std::cout << "Something wrong with the calculated circle! Mis-match at end point." << std::endl; - std::cout << "dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << std::endl; - std::cout << "par1: " << par1 << ", par2: " << par2 << std::endl; - std::cout << "par_start: " << par_start << ", par_end: " << par_end << std::endl; - std::cout << "normal_: " << normal_ << std::endl; - std::cout << "normal_flipped: " << normal_flipped << std::endl; - std::cout << "param_cv->startparam(): " << param_cv->startparam() << std::endl; - std::cout << "param_cv->endparam(): " << param_cv->endparam() << std::endl; - std::cout << "swapped: " << isSwapped() << std::endl; - if (start_par_pt != NULL) { - std::cout << "start_par_pt: " << *start_par_pt << std::endl; - } - if (end_par_pt != NULL) { - std::cout << "end_par_pt: " << *end_par_pt << std::endl; - } - // std::cout << "xvec: " << xvec << std::endl; - // std::cout << "xvec2: " << xvec2 << std::endl; - std::cout << "x_axis_par: " << x_axis_par << std::endl; - std::cout << "space_cv reversed: " << reversed << std::endl; - std::cout << "space_cv normal: " << ((Circle*)(space_crv))->getNormal() << std::endl; - std::cout << "normal_ang: " << normal_ang << std::endl; - std::cout << "y_axis_par: " << y_axis_par << std::endl; - std::cout << "y_axis_par_proj: " << y_axis_par_proj << std::endl; - if (ang_y_axis_par > 0.5*M_PI) { - std::cout << "Flipped: ang_y_axis_par: " << ang_y_axis_par << std::endl; - } else { - std::cout << "Not flipped: ang_y_axis_par: " << ang_y_axis_par << std::endl; - } - } - } #ifdef DEBUG // TEST From 3e439acad4ac2966461f07ebc97a3e26067a3f5e Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 10:38:19 +0100 Subject: [PATCH 18/24] Cleaned up the code. --- gotools-core/src/geometry/Plane.C | 44 ++++++++----------------------- 1 file changed, 11 insertions(+), 33 deletions(-) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index fe9810890..5fca5eb40 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -817,21 +817,19 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, } else { - //return param_cv; - // For the circle we must check: // 1) Space circle normal vs plane normal (may be flipped). - // 2) Plane swapped_. + // 2) Plane swapped_: Handled by the evaluation. // 3) Space curve: Is reversed. + // The space circle normal should be close to the plane normal (or flipped). double ang_tol = 1e-04; Point space_cv_normal = ((Circle*)(space_crv))->getNormal(); double normal_ang = normal_.angle(space_cv_normal); if ((normal_ang > ang_tol) && (normal_ang < M_PI - ang_tol)) { + LOG_WARN("Unexpected normal angle: " + std::to_string(normal_ang)); return param_cv; } - bool normal_flipped = (fabs(M_PI - normal_ang) < ang_tol); - //std::cout << "reversed: " << space_crv->isReversed() << std::endl; Point centre_2d(2); Point centre_3d = ((Circle*)(space_crv))->getCentre(); @@ -841,9 +839,8 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, if (d4 > epspar) return param_cv; - double rad_par = centre_2d.dist(par1); + double rad_par = centre_2d.dist(par1); // The radius in the parametric space. - //std::cout << "rad_par: " << rad_par << ", radius: " << radius << std::endl; Point x_axis = ((Circle*)(space_crv))->getXAxis(); Point x_axis_end = centre_3d + x_axis; double d5; @@ -866,7 +863,6 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, y_axis_par_proj.normalize(); Point y_axis_par(-x_axis_par[1], x_axis_par[0]); // This is the hardcoded assignment of vec2_ in Circle. - //= ((Circle*)(param_cv.get()))->getYAxis(); double ang_y_axis_par = y_axis_par.angle(y_axis_par_proj); bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); @@ -874,36 +870,18 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, //double sign = (reversed) ? -1.0 : 1.0; Point param_cv_axis(0.0, 0.0); - param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, x_axis_par, reversed)); - - param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); - //param_cv->setParameterInterval(space_crv->startparam(), space_crv->endparam()); - - if (y_axis_reversed)// && (sign > 0)) - { - //std::cout << "Reversing y axis." << std::endl; - //((Circle*)(param_cv.get()))->setYAxis(sign*y_axis_par_proj); - //std::cout << "Done reversing y axis." << std::endl; - -#if 0 - std::cout << "WARNING! This approach will fail as the vec2_ i2 forgotten when cloning or writing to file!" - << std::endl; - ((Circle*)(param_cv.get()))->setYAxis(y_axis_par_proj); -#else - //std::cout << "WARNING! Trying to reverse parameter direction to fix vec2_ issue!" << std::endl; + if (!y_axis_reversed) { + param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, x_axis_par, reversed)); + } else { LOG_WARN("Fixing 2d circle with flipped normal."); - param_cv->reverseParameterDirection(); + //param_cv->reverseParameterDirection(); double vec1_rot_ang = 2*M_PI - t2 - t1; // Rotate ccw. Point new_x_axis_par = x_axis_par; Point rot_axis_not_used(0.0, 0.0); // Rot axis is not used for 2d case. Always ccw in the plane. GeometryTools::rotatePoint(rot_axis_not_used, vec1_rot_ang, &new_x_axis_par[0]); param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, new_x_axis_par, !reversed)); - param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); - //((Circle*)param_cv.get() )->setXAxis(y_axis_par_proj); -#endif - - } + param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); // We calculate the dist from the lifted point to the space cv point. Point par_start = param_cv->point(param_cv->startparam()); @@ -918,8 +896,8 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); double tpar_dist = tpar_3d.dist(tpar_sf_3d); if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { - std::cout << "ERROR: dist_par1: " << dist_par1 << ", dist_par2: " << dist_par2 << "tpar_dist: " << - tpar_dist << std::endl; + LOG_WARN("dist_par1: " + std::to_string(dist_par1) + ", dist_par2: " + std::to_string(dist_par2) + + "tpar_dist: " + std::to_string(tpar_dist)); } } #ifdef DEBUG From ba56d25d1dcd26ec38ee8fe220bab07c91bdd51d Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 10:47:30 +0100 Subject: [PATCH 19/24] Cleaned up the code. --- .../include/GoTools/geometry/Circle.h | 21 ------------------- gotools-core/src/geometry/Circle.C | 1 - gotools-core/src/geometry/Plane.C | 1 - 3 files changed, 23 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/Circle.h b/gotools-core/include/GoTools/geometry/Circle.h index 16f2ef049..02ba8f29a 100644 --- a/gotools-core/include/GoTools/geometry/Circle.h +++ b/gotools-core/include/GoTools/geometry/Circle.h @@ -108,15 +108,6 @@ class Circle : public ElementaryCurve virtual double startparam() const; virtual double endparam() const; - double param_bounds_1() const - { - return parbound1_; - } - double param_bounds_2() const - { - return parbound2_; - } - //virtual void reverseParameterDirection(bool switchparam = false); /// The full circle is always parametrized on [0, 2*M_PI), it does @@ -183,18 +174,6 @@ class Circle : public ElementaryCurve return vec2_; } - void setXAxis(const Point& vec1) - { - vec1_ = vec1; - setSpanningVectors(); - } - - void setYAxis(const Point& vec2) - { - vec2_ = vec2; - vec2_.normalize(); - } - double getRadius() const { return radius_; diff --git a/gotools-core/src/geometry/Circle.C b/gotools-core/src/geometry/Circle.C index 8bdafa512..42b89c32a 100644 --- a/gotools-core/src/geometry/Circle.C +++ b/gotools-core/src/geometry/Circle.C @@ -254,7 +254,6 @@ Circle* Circle::clone() const isReversed_); circle->setParamBounds(parbound1_, parbound2_); circle->setParameterInterval(startparam_, endparam_); - circle->setYAxis(circle->getYAxis()); return circle; } diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 5fca5eb40..e699c6cd7 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -867,7 +867,6 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); bool reversed = (space_crv->isReversed()); - //double sign = (reversed) ? -1.0 : 1.0; Point param_cv_axis(0.0, 0.0); if (!y_axis_reversed) { From f8efe49f0a5fe30211c400010de519760663cb47 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 11:10:08 +0100 Subject: [PATCH 20/24] Removed cout. --- .../include/GoTools/geometry/CurveLoop.h | 48 +++---------------- gotools-core/src/geometry/CurveOnSurface.C | 26 +++++----- 2 files changed, 17 insertions(+), 57 deletions(-) diff --git a/gotools-core/include/GoTools/geometry/CurveLoop.h b/gotools-core/include/GoTools/geometry/CurveLoop.h index 7a40fbef6..69fc9fb0d 100644 --- a/gotools-core/include/GoTools/geometry/CurveLoop.h +++ b/gotools-core/include/GoTools/geometry/CurveLoop.h @@ -45,7 +45,6 @@ #include #include "GoTools/utils/config.h" #include "GoTools/geometry/ParamCurve.h" -#include "GoTools/geometry/ElementaryCurve.h" namespace Go @@ -84,50 +83,15 @@ inline double computeLoopGap(const std::vector< PtrToCurveType >& curves) double maxdist = -1.0; double dist; for (i = 1; i < n; ++i) { - double endpar = curves[i-1]->endparam(); - shared_ptr elem_cv_prev = - dynamic_pointer_cast(curves[i-1]); - // if (elem_cv_prev && elem_cv_prev->isReversed()) { - // endpar = curves[i-1]->startparam(); - // } - curves[i-1]->point(endp, endpar); - - double startpar = curves[i]->startparam(); - shared_ptr elem_cv_curr = - dynamic_pointer_cast(curves[i]); - // if (elem_cv_curr && elem_cv_curr->isReversed()) { - // startpar = curves[i]->endparam(); - // } - curves[i]->point(startp, startpar); + curves[i-1]->point(endp, curves[i-1]->endparam()); + curves[i]->point(startp, curves[i]->startparam()); dist = endp.dist(startp); - if (dist > maxdist) { - // if (dist > 0.1) - // std::cout << "i: " << i << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; - maxdist = dist; - } + if (dist > maxdist) maxdist = dist; } - double endpar = curves[n-1]->endparam(); - shared_ptr elem_cv_prev = - dynamic_pointer_cast(curves[n-1]); - // if (elem_cv_prev && elem_cv_prev->isReversed()) { - // endpar = curves[n-1]->startparam(); - // } - curves[n-1]->point(endp, endpar); - - double startpar = curves[0]->startparam(); - shared_ptr elem_cv_curr = - dynamic_pointer_cast(curves[0]); - // if (elem_cv_curr && elem_cv_curr->isReversed()) { - // startpar = curves[0]->endparam(); - // } - curves[0]->point(startp, startpar); + curves[n-1]->point(endp, curves[n-1]->endparam()); + curves[0]->point(startp, curves[0]->startparam()); dist = endp.dist(startp); - if (dist > maxdist) { - // if (dist > 0.1) - // std::cout << "i: " << 0 << ", dist: " << dist << ", elem_cv_prev: " << elem_cv_prev << ", elem_cv_curr: " << elem_cv_curr << std::endl; - maxdist = dist; - } - //std::cout << "maxdist: " << maxdist << std::endl; + if (dist > maxdist) maxdist = dist; return maxdist; } diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 3962fb505..9e87aecf6 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -1507,22 +1507,18 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, } } - if (elem_sf.get() && elem_cv.get()) - { + if (elem_sf.get() && elem_cv.get()) { // The function returns a curve only if the configuration is simple - pcurve_ = elem_sf->getElementaryParamCurve(elem_cv.get(), epspar, - start_par_pt, end_par_pt); - if (pcurve_) { - bool debug_same_trace = sameTrace(epspar); - if (!debug_same_trace) { - LOG_DEBUG("debug_same_trace: " + debug_same_trace); - std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; - } - // else { - // std::cout << "DEBUG: OK: debug_same_trace: " << debug_same_trace << std::endl; - // } - } - } + pcurve_ = elem_sf->getElementaryParamCurve(elem_cv.get(), epspar, + start_par_pt, end_par_pt); + if (pcurve_) { + bool same_trace = sameTrace(epspar); + if (!same_trace) { + LOG_WARN("Projected elementary curve: same_trace: " + same_trace); + //std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; + } + } + } } // If the space curve and surface are not elementary geometry or the parameter curve From 4f981c54d418f59313f19420feab4df0279e4121 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 11:50:13 +0100 Subject: [PATCH 21/24] Cleaned up the code. --- gotools-core/src/geometry/Plane.C | 147 +++++++++++++----------------- 1 file changed, 61 insertions(+), 86 deletions(-) diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index e699c6cd7..271ac8385 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -751,71 +751,51 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, const Point* start_par_pt, const Point* end_par_pt) const //=========================================================================== { - // Default is not simple elementary parameter curve exists - shared_ptr param_cv; - //std::cout << "epspar: " << epspar << std::endl; - // Bookkeeping related to swapped parameters - int ind1 = 0; - int ind2 = 1; - if (isSwapped()) - swap(ind1, ind2); - - double t1, t2; - int idx; - bool closed = false; - if (space_crv->instanceType() == Class_Line) - { - if (!((Line*)(space_crv))->isBounded()) - return param_cv; // Project endpoints onto the surface - t1 = space_crv->startparam(); - t2 = space_crv->endparam(); - idx = ind1; // 0 - } - else if (space_crv->instanceType() == Class_Circle) - { - t1 = space_crv->startparam(); - closed = ((Circle*)(space_crv))->isClosed(); - //std::cout << "closed: " << closed << std::endl; - // t2 = (closed) ? 0.5*(t1 + space_crv->endparam()) : - // space_crv->endparam(); - t2 = space_crv->endparam(); - idx = ind2; // 1 + // Default is not simple elementary parameter curve exists + shared_ptr param_cv; + + if (space_crv->instanceType() == Class_Line) { + if (!((Line*)(space_crv))->isBounded()) + return param_cv; // Project endpoints onto the surface + } else if (space_crv->instanceType() != Class_Circle) { + return param_cv; } - else - return param_cv; - - double fac = (parbound_.umax()-parbound_.umin())/(domain_.umax()-domain_.umin()); - double parval1[2], parval2[2]; - double d1, d2; - Point close1, close2; - Point pos1 = space_crv->ParamCurve::point(t1); - Point pos2 = space_crv->ParamCurve::point(t2); - closestPoint(pos1, parval1[0], parval1[1], close1, d1, epspar); - closestPoint(pos2, parval2[0], parval2[1], close2, d2, epspar); - if (d1 > epspar || d2 > epspar) - return param_cv; - Point par1(parval1[0], parval1[1]); - Point par2(parval2[0], parval2[1]); - if (space_crv->instanceType() == Class_Line) { - // We want L(t1) = par1 && L(t2) = par2. - Point pos = (t2*par1 - t1*par2)/(t2 - t1); - Point dir = (par2 - par1)/(t2 - t1); - - param_cv = shared_ptr(new Line(pos, dir)); - param_cv->setParamBounds(t1, t2); - // We do not reverse the parameter direction of the param cv even if space curve is reversed. Already handled by - // the evaluation. - - Point par_cv_1 = param_cv->point(t1); - double d1 = par1.dist(par_cv_1); - Point par_cv_2 = param_cv->point(t2); - double d2 = par2.dist(par_cv_2); - if (std::max(d1, d2) > epspar) { - std::cout << "DEBUG: Line: d1: " << d1 << ", d2: " << d2 << std::endl; - } + double t1 = space_crv->startparam(); + double t2 = space_crv->endparam(); + + double parval1[2], parval2[2]; + double d1, d2; + Point close1, close2; + Point pos1 = space_crv->ParamCurve::point(t1); + Point pos2 = space_crv->ParamCurve::point(t2); + closestPoint(pos1, parval1[0], parval1[1], close1, d1, epspar); + closestPoint(pos2, parval2[0], parval2[1], close2, d2, epspar); + if (d1 > epspar || d2 > epspar) + return param_cv; + + Point par1(parval1[0], parval1[1]); + Point par2(parval2[0], parval2[1]); + if (space_crv->instanceType() == Class_Line) { + // We want L(t1) = par1 && L(t2) = par2. + Point pos = (t2*par1 - t1*par2)/(t2 - t1); + Point dir = (par2 - par1)/(t2 - t1); + + param_cv = shared_ptr(new Line(pos, dir)); + param_cv->setParamBounds(t1, t2); + // We do not reverse the parameter direction of the param cv even if space curve is reversed. Already handled by + // the evaluation. + + // Issue warning if outside epspar. + Point par_cv_1 = param_cv->point(t1); + double d1 = par1.dist(par_cv_1); + Point par_cv_2 = param_cv->point(t2); + double d2 = par2.dist(par_cv_2); + if (std::max(d1, d2) > epspar) { + LOG_WARN("Line: d1: " + std::to_string(d1) + ", d2: " + std::to_string(d2)); + } } - else + else { // For the circle we must check: // 1) Space circle normal vs plane normal (may be flipped). @@ -833,31 +813,28 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, Point centre_2d(2); Point centre_3d = ((Circle*)(space_crv))->getCentre(); - double d4; - Point close4; - closestPoint(centre_3d, centre_2d[0], centre_2d[1], close4, d4, epspar); - if (d4 > epspar) + double clo_dist; + Point clo_pt; + closestPoint(centre_3d, centre_2d[0], centre_2d[1], clo_pt, clo_dist, epspar); + if (clo_dist > epspar) return param_cv; double rad_par = centre_2d.dist(par1); // The radius in the parametric space. Point x_axis = ((Circle*)(space_crv))->getXAxis(); Point x_axis_end = centre_3d + x_axis; - double d5; - Point close5; Point x_axis_par_end(2); - closestPoint(x_axis_end, x_axis_par_end[0], x_axis_par_end[1], close5, d5, epspar); - if (d5 > epspar) + closestPoint(x_axis_end, x_axis_par_end[0], x_axis_par_end[1], clo_pt, clo_dist, epspar); + if (clo_dist > epspar) return param_cv; Point x_axis_par = x_axis_par_end - centre_2d; + // The 2d circle does not support flipping the normal, using hardcoded assignment of y_axis based on x_axis. Point y_axis = ((Circle*)(space_crv))->getYAxis(); Point y_axis_end = centre_3d + y_axis; - double d6; - Point close6; Point y_axis_par_end(2); - closestPoint(y_axis_end, y_axis_par_end[0], y_axis_par_end[1], close6, d6, epspar); - if (d6 > epspar) + closestPoint(y_axis_end, y_axis_par_end[0], y_axis_par_end[1], clo_pt, clo_dist, epspar); + if (clo_dist > epspar) return param_cv; Point y_axis_par_proj = y_axis_par_end - centre_2d; y_axis_par_proj.normalize(); @@ -867,12 +844,11 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, bool y_axis_reversed = (ang_y_axis_par > 0.5*M_PI); bool reversed = (space_crv->isReversed()); - Point param_cv_axis(0.0, 0.0); if (!y_axis_reversed) { param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, x_axis_par, reversed)); } else { - LOG_WARN("Fixing 2d circle with flipped normal."); + //LOG_WARN("Fixing 2d circle with flipped normal."); //param_cv->reverseParameterDirection(); double vec1_rot_ang = 2*M_PI - t2 - t1; // Rotate ccw. Point new_x_axis_par = x_axis_par; @@ -882,31 +858,30 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, } param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); - // We calculate the dist from the lifted point to the space cv point. + // Add warning to the log if the distance is too large. Point par_start = param_cv->point(param_cv->startparam()); double dist_par1 = par1.dist(par_start); Point par_end = param_cv->point(param_cv->endparam()); double dist_par2 = par2.dist(par_end); - - // Testing dist for intermediate parameter. + // Also testing dist for intermediate parameter. double tpar = (1.0*param_cv->startparam() + 2.0*param_cv->endparam())/3.0; Point tpar_3d = space_crv->point(tpar); Point tpar_2d = param_cv->point(tpar); Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); double tpar_dist = tpar_3d.dist(tpar_sf_3d); if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { - LOG_WARN("dist_par1: " + std::to_string(dist_par1) + ", dist_par2: " + std::to_string(dist_par2) + - "tpar_dist: " + std::to_string(tpar_dist)); + LOG_WARN("Circle: dist_par1: " + std::to_string(dist_par1) + ", dist_par2: " + std::to_string(dist_par2) + + "tpar_dist: " + std::to_string(tpar_dist)); } } #ifdef DEBUG - // TEST - Point p1 = param_cv->ParamCurve::point(param_cv->startparam()); - Point p2 = param_cv->ParamCurve::point(param_cv->endparam()); - int stop_break = 1; + // TEST + Point p1 = param_cv->ParamCurve::point(param_cv->startparam()); + Point p2 = param_cv->ParamCurve::point(param_cv->endparam()); + int stop_break = 1; #endif - return param_cv; + return param_cv; } //=========================================================================== From 05916df08f102223fe2213acb111df139ebf2160 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 12:37:16 +0100 Subject: [PATCH 22/24] getElementaryParamCurve(): Added elem type for cv and sf when projection from space cv failed (trace mismatch). --- gotools-core/src/geometry/CurveOnSurface.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 9e87aecf6..7276ce839 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -1514,7 +1514,9 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, if (pcurve_) { bool same_trace = sameTrace(epspar); if (!same_trace) { - LOG_WARN("Projected elementary curve: same_trace: " + same_trace); + LOG_WARN("Projected elementary curve: same_trace: " + std::to_string(same_trace) + + ", elem_sf->instanceType(): " + std::to_string(elem_sf->instanceType()) + + ", elem_cv->instanceType(): " + std::to_string(elem_cv->instanceType())); //std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; } } From deea0c7fde999beab837d7978ac1a3d162473b4f Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 12:58:49 +0100 Subject: [PATCH 23/24] Cleaned up the code. --- gotools-core/src/geometry/CurveOnSurface.C | 4 +++ gotools-core/src/geometry/Plane.C | 40 ++-------------------- 2 files changed, 7 insertions(+), 37 deletions(-) diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 7276ce839..44e3fe824 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -1519,6 +1519,10 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, ", elem_cv->instanceType(): " + std::to_string(elem_cv->instanceType())); //std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; } + } else { + LOG_WARN("Failed to project elementary curve: elem_sf->instanceType(): " + + std::to_string(elem_sf->instanceType()) + ", elem_cv->instanceType(): " + + std::to_string(elem_cv->instanceType())); } } } diff --git a/gotools-core/src/geometry/Plane.C b/gotools-core/src/geometry/Plane.C index 271ac8385..f530783c0 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -785,22 +785,12 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, param_cv->setParamBounds(t1, t2); // We do not reverse the parameter direction of the param cv even if space curve is reversed. Already handled by // the evaluation. - - // Issue warning if outside epspar. - Point par_cv_1 = param_cv->point(t1); - double d1 = par1.dist(par_cv_1); - Point par_cv_2 = param_cv->point(t2); - double d2 = par2.dist(par_cv_2); - if (std::max(d1, d2) > epspar) { - LOG_WARN("Line: d1: " + std::to_string(d1) + ", d2: " + std::to_string(d2)); - } } else { - // For the circle we must check: + // For the circle we must check (swapped sf is handled by evaluation): // 1) Space circle normal vs plane normal (may be flipped). - // 2) Plane swapped_: Handled by the evaluation. - // 3) Space curve: Is reversed. + // 2) Space curve: Is reversed. // The space circle normal should be close to the plane normal (or flipped). double ang_tol = 1e-04; @@ -848,8 +838,6 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, if (!y_axis_reversed) { param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, x_axis_par, reversed)); } else { - //LOG_WARN("Fixing 2d circle with flipped normal."); - //param_cv->reverseParameterDirection(); double vec1_rot_ang = 2*M_PI - t2 - t1; // Rotate ccw. Point new_x_axis_par = x_axis_par; Point rot_axis_not_used(0.0, 0.0); // Rot axis is not used for 2d case. Always ccw in the plane. @@ -857,29 +845,7 @@ Plane::getElementaryParamCurve(ElementaryCurve* space_crv, double epspar, param_cv = shared_ptr(new Circle(rad_par, centre_2d, param_cv_axis, new_x_axis_par, !reversed)); } param_cv->setParamBounds(space_crv->startparam(), space_crv->endparam()); - - // Add warning to the log if the distance is too large. - Point par_start = param_cv->point(param_cv->startparam()); - double dist_par1 = par1.dist(par_start); - Point par_end = param_cv->point(param_cv->endparam()); - double dist_par2 = par2.dist(par_end); - // Also testing dist for intermediate parameter. - double tpar = (1.0*param_cv->startparam() + 2.0*param_cv->endparam())/3.0; - Point tpar_3d = space_crv->point(tpar); - Point tpar_2d = param_cv->point(tpar); - Point tpar_sf_3d = ParamSurface::point(tpar_2d[0], tpar_2d[1]); - double tpar_dist = tpar_3d.dist(tpar_sf_3d); - if ((dist_par1 > epspar) || (dist_par2 > epspar) || (tpar_dist > epspar)) { - LOG_WARN("Circle: dist_par1: " + std::to_string(dist_par1) + ", dist_par2: " + std::to_string(dist_par2) + - "tpar_dist: " + std::to_string(tpar_dist)); - } - } -#ifdef DEBUG - // TEST - Point p1 = param_cv->ParamCurve::point(param_cv->startparam()); - Point p2 = param_cv->ParamCurve::point(param_cv->endparam()); - int stop_break = 1; -#endif + } return param_cv; } From 7430afab5aa0adfa7efc96312326ca51a22d4872 Mon Sep 17 00:00:00 2001 From: Sverre Briseid Date: Sat, 2 Nov 2024 13:05:37 +0100 Subject: [PATCH 24/24] Resetting param cv if mismatch between trace of space and param cv. --- gotools-core/src/geometry/CurveOnSurface.C | 1 + 1 file changed, 1 insertion(+) diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index 44e3fe824..9adf92d2b 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -1518,6 +1518,7 @@ bool CurveOnSurface::ensureParCrvExistence(double epsgeo, ", elem_sf->instanceType(): " + std::to_string(elem_sf->instanceType()) + ", elem_cv->instanceType(): " + std::to_string(elem_cv->instanceType())); //std::cout << "DEBUG: debug_same_trace: " << debug_same_trace << std::endl; + pcurve_ = shared_ptr(); } } else { LOG_WARN("Failed to project elementary curve: elem_sf->instanceType(): " +