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..91c92fc5e 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 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. /// \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..02ba8f29a 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, @@ -168,6 +169,11 @@ class Circle : public ElementaryCurve return vec1_; } + Point getYAxis() const + { + return vec2_; + } + double getRadius() const { return radius_; 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/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/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..2366774f3 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 @@ -179,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_; @@ -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..a7fa2733c 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. @@ -75,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/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/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 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/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 4baede2fb..064f92ad0 100644 --- a/gotools-core/src/geometry/BoundedUtils.C +++ b/gotools-core/src/geometry/BoundedUtils.C @@ -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 radius1, + double radius2, 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(), 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, + "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(), radius1, radius2, 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, @@ -4367,11 +5002,17 @@ 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) { - 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)); } } @@ -4737,7 +5378,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/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..a371f1e7e 100644 --- a/gotools-core/src/geometry/Cone.C +++ b/gotools-core/src/geometry/Cone.C @@ -696,10 +696,18 @@ 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 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()); @@ -713,10 +721,18 @@ 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 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; @@ -890,7 +906,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 +1356,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 +1423,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 +1465,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/CurveLoop.C b/gotools-core/src/geometry/CurveLoop.C index 70b8ce7dc..29324d2a5 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 @@ -424,6 +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_)) { + 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; } @@ -548,6 +551,12 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) } double maxgap_par = computeLoopGap(par_cvs); double maxgap_space = computeLoopGap(space_cvs); + if (maxgap_space > space_epsilon_) { + LOG_WARN("maxgap_space: " + std::to_string(maxgap_space)); + } + if (maxgap_par > space_epsilon_) { + LOG_WARN("maxgap_par: " + std::to_string(maxgap_par)); + } if ((maxgap_space >= 0.0) && (maxgap_space < space_epsilon_)) { for (int k = 0; k < (int)curves.size(); ++k) @@ -583,15 +592,16 @@ bool CurveLoop::fixInvalidLoop(double& max_gap) // If we got closer we replace the loop curves. max_gap2 = computeLoopGap(curves); if (max_gap2 < max_gap) - { + { + LOG_WARN("Replacing the loop curves."); curves_ = curves; max_gap = max_gap2; - } + } #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/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index a3b7637fe..9adf92d2b 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -55,9 +55,11 @@ #include "GoTools/creators/HermiteAppS.h" #include "GoTools/creators/CurveCreators.h" #include "GoTools/creators/CoonsPatchGen.h" +#include "GoTools/utils/Logger.h" #include #include + using namespace Go; using std::vector; using std::max; @@ -149,7 +151,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 +815,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 +872,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 +896,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,14 +1498,34 @@ 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()) - { + 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); - } + 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: " + 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; + pcurve_ = shared_ptr(); + } + } 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())); + } + } } // If the space curve and surface are not elementary geometry or the parameter curve @@ -1508,47 +1537,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 +2394,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..f530783c0 100644 --- a/gotools-core/src/geometry/Plane.C +++ b/gotools-core/src/geometry/Plane.C @@ -39,7 +39,10 @@ #include "GoTools/geometry/Plane.h" #include "GoTools/geometry/Line.h" +#include "GoTools/geometry/Circle.h" #include "GoTools/geometry/SplineSurface.h" +#include "GoTools/geometry/GeometryTools.h" +#include "GoTools/utils/Logger.h" #include #include @@ -50,6 +53,7 @@ using std::numeric_limits; using std::streamsize; using std::swap; +//#define DEBUG namespace Go { @@ -741,6 +745,110 @@ bool Plane::isDegenerate(bool& b, bool& r, return false; } +//=========================================================================== +shared_ptr +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; + + 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; + } + + 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. + } + else + { + // For the circle we must check (swapped sf is handled by evaluation): + // 1) Space circle normal vs plane normal (may be flipped). + // 2) 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; + } + + Point centre_2d(2); + Point centre_3d = ((Circle*)(space_crv))->getCentre(); + 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; + Point x_axis_par_end(2); + 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; + Point y_axis_par_end(2); + 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(); + + Point y_axis_par(-x_axis_par[1], x_axis_par[0]); // This is the hardcoded assignment of vec2_ in Circle. + 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()); + 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 { + 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()); + } + + return param_cv; +} //=========================================================================== void Plane::getDegenerateCorners(vector& deg_corners, double tol) const @@ -980,6 +1088,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..c7ad8e251 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) @@ -132,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 = @@ -229,11 +306,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); diff --git a/gotools-core/src/geometry/SurfaceTools.C b/gotools-core/src/geometry/SurfaceTools.C index 7c6bcfb48..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; @@ -186,7 +187,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; } } @@ -816,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 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/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);