From 7ecc97f0d68671aa062ebb01a0e3e91727798fac Mon Sep 17 00:00:00 2001 From: vsk Date: Wed, 23 Sep 2026 09:51:28 +0200 Subject: [PATCH] Update with accumulated changes --- .../GoTools/compositemodel/CurveModel.h | 12 ++ .../GoTools/compositemodel/SurfaceModel.h | 23 ++- .../GoTools/compositemodel/ftPointSet.h | 9 + compositemodel/src/AdaptSurface.C | 4 +- compositemodel/src/CurveModel.C | 52 ++++-- compositemodel/src/SurfaceModel.C | 161 +++++++++++++++++- compositemodel/src/Vertex.C | 3 + .../include/GoTools/creators/ApproxCurve.h | 11 +- .../include/GoTools/geometry/BoundedCurve.h | 6 + .../GoTools/geometry/CurveBoundedDomain.h | 2 +- .../include/GoTools/geometry/CurveOnSurface.h | 8 +- .../include/GoTools/geometry/Domain.h | 4 + .../include/GoTools/geometry/ParamCurve.h | 6 + .../include/GoTools/geometry/RectDomain.h | 4 +- .../include/GoTools/geometry/SplineCurve.h | 2 +- .../GoTools/geometry/SplineDebugUtils.h | 5 +- .../include/GoTools/utils/QRFactorization.h | 58 +++++++ gotools-core/src/creators/ApproxCurve.C | 6 +- gotools-core/src/geometry/CurveOnSurface.C | 2 +- gotools-core/src/geometry/ParamCurve.C | 24 +++ gotools-core/src/geometry/SplineDebugUtils.C | 25 +++ gotools-core/src/utils/BoundingBox.C | 2 +- gotools-core/src/utils/QRFactorization.C | 144 ++++++++++++++++ .../GoTools/trivariate/CurveOnVolume.h | 6 + 24 files changed, 539 insertions(+), 40 deletions(-) create mode 100644 gotools-core/include/GoTools/utils/QRFactorization.h create mode 100644 gotools-core/src/utils/QRFactorization.C diff --git a/compositemodel/include/GoTools/compositemodel/CurveModel.h b/compositemodel/include/GoTools/compositemodel/CurveModel.h index 74fd89474..4fc2109b9 100644 --- a/compositemodel/include/GoTools/compositemodel/CurveModel.h +++ b/compositemodel/include/GoTools/compositemodel/CurveModel.h @@ -53,6 +53,7 @@ namespace Go { class CompositeCurve; + class Vertex; //=========================================================================== /** A curve model including topological information @@ -107,6 +108,12 @@ class CurveModel : public CompositeModel /// \return Index to curve int getIndex(ParamCurve* curve) const; + /// Fetch all edges + std::vector > allEdges() + { + return edges_; + } + /// Evaluate position /// \param idx Index of curve /// \param par[] Parameter value @@ -225,6 +232,11 @@ class CurveModel : public CompositeModel /// \return Vector of pointers to the composite curves std::vector > fetchCompositeCurves() const; + /// Return all vertices associated with this surface model + /// \retval vertices Vector of pointers to all vertices. + void getVertices(std::vector >& vertices) const; + + private: std::vector > edges_; diff --git a/compositemodel/include/GoTools/compositemodel/SurfaceModel.h b/compositemodel/include/GoTools/compositemodel/SurfaceModel.h index 632f1e37b..104fae2df 100644 --- a/compositemodel/include/GoTools/compositemodel/SurfaceModel.h +++ b/compositemodel/include/GoTools/compositemodel/SurfaceModel.h @@ -285,6 +285,21 @@ class GO_API SurfaceModel : public CompositeModel std::vector& der) const; // Result + /// Closest point between a given point and the outer boundary/boundaries + /// of this surface model + /// Returns one point + /// \param pnt Input point + /// \param clo_pnt Found closest point + /// \param idx Index of surface where the closest point is found + /// \param clo_par[] Parameter value corresponding to the closest point + /// \param dist Distance between input point and found closest point + void + closestBoundaryPoint(Point& pnt, // Input point + Point& clo_pnt, // Found closest point + int& idx, // Index of surface where the closest point is found + double clo_par[], // Parameter value corresponding to the closest point + double& dist); // Distance between input point and found closest point + /// Closest point between a given point and this surface model /// Returns one point /// \param pnt Input point @@ -304,7 +319,11 @@ class GO_API SurfaceModel : public CompositeModel /// \return Closest point ftPoint closestPoint(const Point& point); - /// Closest point between a given point and this surface model + void closestPoint(Point& point, int seed_ix, double seed[], + Point& clo_pt, int& idx, double clo_par[], + double& dist); + + /// Closest point between a given point and this surface model /// \param point Input point /// \return Closest point ftPoint closestPoint(const ftPoint& point) { return closestPoint(point.position()); } @@ -913,7 +932,7 @@ class GO_API SurfaceModel : public CompositeModel std::vector >& crv_bound, bool compute_curves=true) const; - ftPoint closestPointLocal(const ftPoint& point) const; + ftPoint closestPointLocal(const ftPoint& point, bool use_seed=false) const; void localExtreme(ftSurface *face, Point& dir, Point& ext_pnt, int& ext_id, diff --git a/compositemodel/include/GoTools/compositemodel/ftPointSet.h b/compositemodel/include/GoTools/compositemodel/ftPointSet.h index fb9b3e2ed..608239b4b 100644 --- a/compositemodel/include/GoTools/compositemodel/ftPointSet.h +++ b/compositemodel/include/GoTools/compositemodel/ftPointSet.h @@ -118,6 +118,15 @@ class ftSamplePoint /// Am I on sub surface boundary? bool isOnSubSurfaceBoundary() const { return ((at_boundary_ == 2) || (at_boundary_ == 1)); } + /// Am I on an inner boundary? + bool isOnInnerBoundary() const + { return at_boundary_ == 2; } + + virtual bool isCorner() const + { + return false; + } + /// Get number of neighbours. int getNmbNeighbour() const { return (int)next_.size();} diff --git a/compositemodel/src/AdaptSurface.C b/compositemodel/src/AdaptSurface.C index 5ed637e30..a32752a72 100644 --- a/compositemodel/src/AdaptSurface.C +++ b/compositemodel/src/AdaptSurface.C @@ -1010,7 +1010,7 @@ namespace Go // Fetch constant parameter curves in the selected curve direction double u1 = (cv_dir == 0) ? dom.umin() : dom.vmin(); double u2 = (cv_dir == 0) ? dom.umax() : dom.vmax(); - double udel = (u2 - u1)/(double)(nmb_cv+1); + double udel = (u2 - u1)/(double)(nmb_cv+2); double tol1 = std::max(1.0e-7, 1.0e-5*(u2-u1)); double par[2]; //int ki, kj; @@ -1071,7 +1071,7 @@ namespace Go nmb = std::max(nmb, min_samples); double v1 = crvs[kr]->startparam(); double v2 = crvs[kr]->endparam(); - double vdel = (v2 - v1)/(double)(nmb+1); + double vdel = (v2 - v1)/(double)(nmb+2); double tol2 = std::max(1.0e-7, 1.0e-5*(v2-v1)); if (consider_joint) par[pt_dir] = std::min(v1+vdel, surf->nextSegmentVal(pt_dir, v1, true, tol2)); diff --git a/compositemodel/src/CurveModel.C b/compositemodel/src/CurveModel.C index 0eacd8abd..37b1a25a6 100644 --- a/compositemodel/src/CurveModel.C +++ b/compositemodel/src/CurveModel.C @@ -47,8 +47,8 @@ using std::vector; -namespace Go -{ +using namespace Go; + //=========================================================================== CurveModel::CurveModel(double gap, // Gap between adjacent curves double neighbour, // Threshold for whether curves are adjacent @@ -391,20 +391,40 @@ void CurveModel::tesselate(vector >& meshes) const vector > vertices(all_vertices.begin(), all_vertices.end()); for (ki=0; kigetVertexPoint().dist(vertices[kj]->getVertexPoint()); - if (dist < toptol_.neighbour) - { - vector connected_edges = vertices[kj]->allEdges(); - vertices[ki]->joinVertex(vertices[kj]); - for (size_t kr=0; krreplaceVertex(vertices[kj], vertices[ki]); - vertices.erase(vertices.begin() + kj); - break; - } + { + for (kj=ki+1; kjgetVertexPoint().dist(vertices[kj]->getVertexPoint()); + if (dist < toptol_.neighbour) + { + vector connected_edges = vertices[kj]->allEdges(); + vertices[ki]->joinVertex(vertices[kj]); + for (size_t kr=0; krreplaceVertex(vertices[kj], vertices[ki]); + vertices.erase(vertices.begin() + kj); + } + else + ++kj; + } } } -} // namespace Go + + //=========================================================================== + void CurveModel::getVertices(vector >& vertices) const + //=========================================================================== + { + std::set > all_vertices; // All vertices in the model represented once + for (size_t ki=0; ki vert = edges_[ki]->getVertex(true); + all_vertices.insert(vert); + vert = edges_[ki]->getVertex(false); + all_vertices.insert(vert); + } + vertices.clear(); + vertices.insert(vertices.end(), all_vertices.begin(), all_vertices.end()); + } + + diff --git a/compositemodel/src/SurfaceModel.C b/compositemodel/src/SurfaceModel.C index cd4852f13..160daf25e 100644 --- a/compositemodel/src/SurfaceModel.C +++ b/compositemodel/src/SurfaceModel.C @@ -746,6 +746,53 @@ namespace Go #endif } + //=========================================================================== + void SurfaceModel::closestBoundaryPoint(Point& pnt, // Input point + Point& clo_pnt, // Found closest point + int& idx, // Index of surface where the closest point is found + double clo_par[], // Parameter value corresponding to the closest point + double& dist) // Distance between input point and found closest point + //=========================================================================== + { + idx = -1; + int loop_ix1 = -1, loop_ix2 = -1; + size_t cv_ix; + dist = std::numeric_limits::max(); + double cv_par; + for (size_t ki=0; kiclosestPoint(pnt, cv_ind, tpar, + cv_close, tdist); + if (tdist < dist) + { + cv_ix = cv_ind; + loop_ix1 = (int)ki; + loop_ix2 = (int)kj; + cv_par = tpar; + dist = tdist; + clo_pnt = cv_close; + } + } + } + + if (loop_ix1 >= 0) + { + shared_ptr edg0 = + boundary_curves_[loop_ix1][loop_ix2]->getEdge(cv_ix); + ftEdge* edg = edg0->geomEdge(); + ftFaceBase *face0 = edg->face(); + ftSurface *face = face0->asFtSurface(); + idx = getIndex(face); + Point fpar = edg->faceParameter(cv_par); + clo_par[0] = fpar[0]; + clo_par[1] = fpar[1]; + } + } @@ -870,7 +917,34 @@ namespace Go return ftPoint(bestcp, bestface, bestu, bestv); } - + //=========================================================================== + void SurfaceModel::closestPoint(Point& point, int seed_ix, double seed[], + Point& clo_pt, int& idx, double clo_par[], + double& dist) + //=========================================================================== + { + fill(face_checked_.begin(), face_checked_.end(), false); + ftPoint inpt; + bool use_seed = (seed_ix >= 0); + int ix = (seed_ix >= 0) ? seed_ix : 0; + ftSurface *first = static_pointer_cast(faces_[ix]).get(); + if (seed_ix >= 0) + inpt = ftPoint(point, first, seed[0], seed[1]); + else + inpt = ftPoint(point, first); + ftPoint close = closestPointLocal(inpt, use_seed); + if (close.face() != 0) + { + clo_pt = close.position(); + idx = getIndex(close.face()); + clo_par[0] = close.u(); + clo_par[1] = close.v(); + dist = point.dist(clo_pt); + } + else + closestPoint(point, clo_pt, idx, clo_par, dist); + } + //=========================================================================== int SurfaceModel::nmbEntities() const @@ -1397,9 +1471,10 @@ void SurfaceModel::swapFaces(int idx1, int idx2) } //=========================================================================== - ftPoint SurfaceModel::closestPointLocal(const ftPoint& point) const + ftPoint SurfaceModel::closestPointLocal(const ftPoint& point, bool use_seed) const //=========================================================================== { + double tol = toptol_.gap; // Should be a parameter tolerance const Point& pt = point.position(); int id = getIndex(point.face()); ftSurface* curface = 0; @@ -1413,13 +1488,33 @@ void SurfaceModel::swapFaces(int idx1, int idx2) double closestpt_epsilon = toptol_.neighbour; // Maybe gap instead? ftSurface* bestface = 0; int nmb_checked = 0; + double seedval[2]; + double *seed = use_seed ? seedval : 0; + seedval[0] = point.u(); + seedval[1] = point.v(); while (!finished) { if (!face_checked_[id]) { curface = dynamic_cast(faces_[id].get()); face_checked_[id] = true; // cout << "Face: " << id << endl; ASSERT(curface != 0); - curface->closestPoint(pt, u, v, cp, dist, closestpt_epsilon); + const Domain& dom = curface->surface()->parameterDomain(); + if (seed) + { + bool in_domain = false; + try { + Vector2D seed2(seed[0], seed[1]); + in_domain = dom.isInDomain(seed2, tol); + } + catch (...) + { + in_domain = false; + } + if (!in_domain) + seed = 0; + } + curface->closestPoint(pt, u, v, cp, dist, closestpt_epsilon, NULL, + seed); nmb_checked++; if (dist < bestdist) { bestdist = dist; @@ -1429,7 +1524,6 @@ void SurfaceModel::swapFaces(int idx1, int idx2) bestface = curface; } // Check if the point was on the boundary - const Domain& dom = curface->surface()->parameterDomain(); bool on_boundary = dom.isOnBoundary(Vector2D(u, v), toptol_.neighbour); if (on_boundary) { @@ -1438,11 +1532,59 @@ void SurfaceModel::swapFaces(int idx1, int idx2) = curface->edgeClosestToPoint(u, v); ftEdgeBase* twin = boundary_edge->twin(); if (!twin) // That is, there is no neighbour - finished = true; - else { + { + bool corner = dom.isOnCorner(Vector2D(u, v), + toptol_.neighbour); + if (corner) + { + shared_ptr vx1, vx2; + boundary_edge->geomEdge()->getVertices(vx1, vx2); + Point pt1 = vx1->getVertexPoint(); + Point pt2 = vx2->getVertexPoint(); + shared_ptr vx = (pt1.dist(pt) <= pt2.dist(pt)) ? vx1 : vx2; + vector edgs = vx->getEdges(curface); + for (size_t kr=0; krtwin(); + if (twin) + break; + } + } + } + if (!twin) + finished = true; + } + if (twin) + { // cout << "We're crossing a boundary!" << endl; - id = getIndex(twin->face()->asFtSurface()); - } + int id2 = getIndex(twin->face()->asFtSurface()); + if (id2 == id) + { + shared_ptr surf = getSurface(id); + SurfaceTools::surface_seedfind(pt, *surf, 0, + seedval[0], seedval[1]); + } + else + { + id = id2; + if (use_seed) + { + ftEdge *twin2 = twin->geomEdge(); + if (twin2) + { + double tpar, tdist; + Point close; + twin2->closestPoint(bestcp, tpar, close, tdist); + Point fpar = twin2->faceParameter(tpar); + seedval[0] = fpar[0]; + seedval[1] = fpar[1]; + seed = seedval; + } + } + } + } } else // point was in the interior finished = true; } else { // if face_checked_[id] @@ -2149,7 +2291,8 @@ void SurfaceModel::swapFaces(int idx1, int idx2) } } -//=========================================================================== + +///=========================================================================== vector > SurfaceModel::getBoundaryEdges() const //=========================================================================== { diff --git a/compositemodel/src/Vertex.C b/compositemodel/src/Vertex.C index e98cba76e..f718fa8b1 100644 --- a/compositemodel/src/Vertex.C +++ b/compositemodel/src/Vertex.C @@ -841,6 +841,9 @@ namespace Go for (size_t ki=0; kiface()->asFtSurface(); + shared_ptr surf = curr_face->surface(); + if (!surf.get()) + continue; size_t kj; for (kj=0; kj& points, const std::vector& parvals, @@ -106,7 +106,7 @@ class ApproxCurve /// approximate /// \param dim the spatial dimension of the points (usually 3) /// \param aepsge the geometric tolerance to work with - /// \param in the number of control points of the resulting spline curve + /// \param in the number of control points of the initial spline curve /// \param ik the order of the resulting spline curve (pol. degree + 1) /// \param knots specifies the knotvector of the resulting spline curve. ApproxCurve(const std::vector& points, @@ -142,6 +142,12 @@ class ApproxCurve /// Approximate C1 continuity with a given impartance 0 <= fac < 1 void setC1Approx(double fac); + /// Flag for parameter iteration. Default: true for dimension > 1 + void setParameterIteration(bool repar) + { + repar_ = repar; + } + /// When everything else is set, this function can be used to fetch the /// approximating curve /// \retval maxdist report the maximum distance between the generated curve and @@ -166,6 +172,7 @@ class ApproxCurve double smoothweight_; double smoothfac_; double c1fac_; + bool repar_; int dim_; std::vector points_; diff --git a/gotools-core/include/GoTools/geometry/BoundedCurve.h b/gotools-core/include/GoTools/geometry/BoundedCurve.h index 6799c0d34..bd6abd79e 100644 --- a/gotools-core/include/GoTools/geometry/BoundedCurve.h +++ b/gotools-core/include/GoTools/geometry/BoundedCurve.h @@ -121,6 +121,12 @@ class GO_API BoundedCurve : public ParamCurve /// \param endpar end parameter virtual void setParameterInterval(double t1, double t2); + // Translate the curve along a given vector + virtual void translateCurve(const Point& dir) + { + curve_->translateCurve(dir); + } + virtual SplineCurve* geometryCurve(); virtual bool isDegenerate(double degenerate_epsilon); diff --git a/gotools-core/include/GoTools/geometry/CurveBoundedDomain.h b/gotools-core/include/GoTools/geometry/CurveBoundedDomain.h index 8d975318f..4b0808550 100644 --- a/gotools-core/include/GoTools/geometry/CurveBoundedDomain.h +++ b/gotools-core/include/GoTools/geometry/CurveBoundedDomain.h @@ -125,7 +125,7 @@ class GO_API CurveBoundedDomain : public Domain /// Check if the given parameter pair is located on the endpoint of some /// curve in the curve loop - bool isOnCorner(const Array& point, + virtual bool isOnCorner(const Array& point, double tolerance) const; /// Find the parameter pair contained in the domain that is closest (using diff --git a/gotools-core/include/GoTools/geometry/CurveOnSurface.h b/gotools-core/include/GoTools/geometry/CurveOnSurface.h index 60a0aa637..c376337d6 100644 --- a/gotools-core/include/GoTools/geometry/CurveOnSurface.h +++ b/gotools-core/include/GoTools/geometry/CurveOnSurface.h @@ -273,7 +273,13 @@ class GO_API CurveOnSurface : public ParamCurve spacecurve_ = spacecurve; } - /// Replace the parameter curve corresponding to this curve on surface curve. + // Translate the curve along a given vector + virtual void translateCurve(const Point& dir) + { + spacecurve_->translateCurve(dir); + } + + /// Replace the parameter curve corresponding to this curve on surface curve. /// Used for instance in relation to reparameterizations of the related surface. /// Use with care! void setParameterCurve(shared_ptr parametercurve) diff --git a/gotools-core/include/GoTools/geometry/Domain.h b/gotools-core/include/GoTools/geometry/Domain.h index f0d13fa2f..834bde3b7 100644 --- a/gotools-core/include/GoTools/geometry/Domain.h +++ b/gotools-core/include/GoTools/geometry/Domain.h @@ -87,6 +87,10 @@ class GO_API Domain virtual bool isOnBoundary(const Array& point, double tolerance) const = 0; + /// Check if the given parameter pair is located at a domain corner + virtual bool isOnCorner(const Array& point, + double tolerance) const = 0; + /// Find the (u, v) point in the Domain that is closest (using Euclidean distance /// in R^2) to a given (u, v) point. If the given point is in the domain, then /// the answer is obviously the same point. diff --git a/gotools-core/include/GoTools/geometry/ParamCurve.h b/gotools-core/include/GoTools/geometry/ParamCurve.h index db5fd109e..44efc9b51 100644 --- a/gotools-core/include/GoTools/geometry/ParamCurve.h +++ b/gotools-core/include/GoTools/geometry/ParamCurve.h @@ -113,6 +113,9 @@ class GO_API ParamCurve : public GeomObject /// put a 'using ParamCurve::point' in the class definition. std::vector point(double tpar, int derivs, bool from_right = true) const; + /// Curvature radius of curve in a given point + double curvatureRadius(double tpar, bool from_right = true) const; + /// Evaluate points on a regular set of parameter values /// \param num number of values to evaluate /// \param points upon function return, this vector holds all the evaluated points @@ -139,6 +142,9 @@ class GO_API ParamCurve : public GeomObject /// Linear reparametrization. The meaning is changed for elementary curves virtual void setParameterInterval(double t1, double t2) = 0; + // Translate the curve along a given vector + virtual void translateCurve(const Point& dir) = 0; + /// If the definition of this ParamCurve contains a SplineCurve describing its /// spatial shape, then this function will return a pointer to this SplineCurve. /// Otherwise it will return a null pointer. diff --git a/gotools-core/include/GoTools/geometry/RectDomain.h b/gotools-core/include/GoTools/geometry/RectDomain.h index 2d6b65a9e..6c267bfc1 100644 --- a/gotools-core/include/GoTools/geometry/RectDomain.h +++ b/gotools-core/include/GoTools/geometry/RectDomain.h @@ -77,7 +77,7 @@ class GO_API RectDomain : public Domain virtual int isInDomain2(const Array& point, double tolerance) const; - // check whether a gien parameter pair is located on the Domain boundary. + // check whether a given parameter pair is located on the Domain boundary. // DOXYGEN documentation can be found in the base class header Domain.h virtual bool isOnBoundary(const Array& point, double tolerance) const; @@ -94,7 +94,7 @@ class GO_API RectDomain : public Domain /// Check if a given parameter pair lies on a corner in the domain within /// the given tolerance - bool isOnCorner(const Array& point, + virtual bool isOnCorner(const Array& point, double tolerance) const; /// Given two parameter pairs, check if they specify a domain boundary diff --git a/gotools-core/include/GoTools/geometry/SplineCurve.h b/gotools-core/include/GoTools/geometry/SplineCurve.h index 769cbabb7..73ee5ea30 100644 --- a/gotools-core/include/GoTools/geometry/SplineCurve.h +++ b/gotools-core/include/GoTools/geometry/SplineCurve.h @@ -546,7 +546,7 @@ class GO_API SplineCurve : public ParamCurve void replaceEndPoint(Point pnt, bool at_start); // Translate the curve along a given vector - void translateCurve(const Point& dir); + virtual void translateCurve(const Point& dir); // Translate the curve along a given vector and swap sign if specified /// Modify in 1. (pdir == 1), 2. (pdir == 2) or both (pdir == 3) diff --git a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h index a7fa2733c..0db28a721 100644 --- a/gotools-core/include/GoTools/geometry/SplineDebugUtils.h +++ b/gotools-core/include/GoTools/geometry/SplineDebugUtils.h @@ -78,7 +78,10 @@ namespace SplineDebugUtils 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 + void GO_API writeSpaceParamSurf(const SplineSurface& psurf, + 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. void GO_API writeTrimmedInfo(BoundedSurface& bd_sf, diff --git a/gotools-core/include/GoTools/utils/QRFactorization.h b/gotools-core/include/GoTools/utils/QRFactorization.h new file mode 100644 index 000000000..b7d81d694 --- /dev/null +++ b/gotools-core/include/GoTools/utils/QRFactorization.h @@ -0,0 +1,58 @@ +/* + * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, + * Applied Mathematics, Norway. + * + * Contact information: E-mail: tor.dokken@sintef.no + * SINTEF ICT, Department of Applied Mathematics, + * P.O. Box 124 Blindern, + * 0314 Oslo, Norway. + * + * This file is part of GoTools. + * + * GoTools is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * GoTools is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public + * License along with GoTools. If not, see + * . + * + * In accordance with Section 7(b) of the GNU Affero General Public + * License, a covered work must retain the producer line in every data + * file that is created or manipulated using GoTools. + * + * Other Usage + * You can be released from the requirements of the license by purchasing + * a commercial license. Buying such a license is mandatory as soon as you + * develop commercial activities involving the GoTools library without + * disclosing the source code of your own applications. + * + * This file may be used in accordance with the terms contained in a + * written agreement between you and SINTEF ICT. + */ + +#ifndef _QRFACTORIZATION_H +#define _QRFACTORIZATION_H + +#include + +namespace Go +{ + + namespace QRFactorization + { + void QRDecomp(std::vector& A, int n, int m, + std::vector& Q, std::vector& R); + + void QRSolve(std::vector& Q, std::vector& R, int n, int m, + std::vector& b, int dim, std::vector& x); + }; +}; + +#endif diff --git a/gotools-core/src/creators/ApproxCurve.C b/gotools-core/src/creators/ApproxCurve.C index da19b00a2..2456310f5 100644 --- a/gotools-core/src/creators/ApproxCurve.C +++ b/gotools-core/src/creators/ApproxCurve.C @@ -73,6 +73,7 @@ ApproxCurve::ApproxCurve() smoothweight_ = 0.000000001; smoothfac_ = 1.0; c1fac_ = 0.0; + repar_ = true; } //*************************************************************************** @@ -96,6 +97,7 @@ ApproxCurve::ApproxCurve(const std::vector& points, aepsge_ = aepsge; smoothweight_ = 0.000000001; c1fac_ = 0.0; + repar_ = (dim > 1); points_.reserve(points.size()); parvals_.reserve(parvals.size()); @@ -133,6 +135,7 @@ ApproxCurve::ApproxCurve(const std::vector& points, aepsge_ = aepsge; smoothweight_ = 0.000000001; c1fac_ = 0.0; + repar_ = (dim > 1); points_.reserve(points.size()); parvals_.reserve(parvals.size()); @@ -168,6 +171,7 @@ ApproxCurve::ApproxCurve(const std::vector& points, aepsge_ = aepsge; smoothweight_ = 0.000000001; c1fac_ = 0.0; + repar_ = (dim > 1); points_.reserve(points.size()); parvals_.reserve(parvals.size()); @@ -423,7 +427,7 @@ void ApproxCurve::checkAccuracy(std::vector& newknots, int uniform) curr_crv_->writeStandardHeader(of); curr_crv_->write(of); #endif - bool reparam = (dim_ == 1) ? false : true; + bool reparam = (dim_ == 1) ? false : repar_; // double par_tol = 0.000000000001; maxdist_ = -10000.0; diff --git a/gotools-core/src/geometry/CurveOnSurface.C b/gotools-core/src/geometry/CurveOnSurface.C index a602f3c1b..4f04b817a 100644 --- a/gotools-core/src/geometry/CurveOnSurface.C +++ b/gotools-core/src/geometry/CurveOnSurface.C @@ -2298,6 +2298,7 @@ Point CurveOnSurface::faceParameter(double crv_par, if (pcurve_.get()) { param = pcurve_->ParamCurve::point(crv_par); + seed = param.begin(); // crv_par = param[2-constdir_]; // same = true; } @@ -2306,7 +2307,6 @@ Point CurveOnSurface::faceParameter(double crv_par, param[constdir_-1] = constval_; param[2-constdir_] = crv_par; } - seed = param.begin(); } else if (pcurve_.get()) { diff --git a/gotools-core/src/geometry/ParamCurve.C b/gotools-core/src/geometry/ParamCurve.C index d51c6b865..4ae215806 100644 --- a/gotools-core/src/geometry/ParamCurve.C +++ b/gotools-core/src/geometry/ParamCurve.C @@ -129,6 +129,30 @@ ParamCurve::point(double tpar, return pts; } +//=========================================================================== + double ParamCurve::curvatureRadius(double tpar, bool from_right) const +//=========================================================================== + { + double eps = 1.0e-10; + int dim = dimension(); + if (dim > 3) + return -1.0; + int derivs = 2; + std::vector der(derivs+1, Point()); + point(der, tpar, derivs, from_right); + double kappa; + if (dim == 1) + kappa = fabs(der[2][0])/pow(1.0+der[1].length(), 3); + else if (dim == 2) + kappa = pow(der[1].length(), 3)/fabs(der[1][0]*der[2][1]-der[1][1]*der[2][0]); + else + { + Point vec = der[1].cross(der[2]); + kappa = vec.length()/pow(der[1].length(), 3); + } + return (kappa < eps) ? -1 : 1.0/kappa; + } + //=========================================================================== bool ParamCurve::isClosed() diff --git a/gotools-core/src/geometry/SplineDebugUtils.C b/gotools-core/src/geometry/SplineDebugUtils.C index c7ad8e251..4fff7605b 100644 --- a/gotools-core/src/geometry/SplineDebugUtils.C +++ b/gotools-core/src/geometry/SplineDebugUtils.C @@ -170,6 +170,31 @@ void SplineDebugUtils::writeSpaceParamCurve(shared_ptr pcurve, } } +//=========================================================================== +void SplineDebugUtils::writeSpaceParamSurf(const SplineSurface& psurf, std::ostream& os, + double z) +//=========================================================================== +{ + ALWAYS_ERROR_IF(psurf.dimension() != 2, + "Expecting input of 2D-curve."); + + std::vector space_coefs; + for (int i = 0; i < psurf.numCoefs_u()*psurf.numCoefs_v(); ++i) { + space_coefs.insert(space_coefs.end(), + psurf.coefs_begin() + i*2, + psurf.coefs_begin() + (i + 1)*2); + space_coefs.push_back(z); // Make param_curve live in plane parallell to the xy-plane. + } + + SplineSurface space_psurf = + SplineSurface(psurf.numCoefs_u(), psurf.numCoefs_v(), + psurf.order_u(), psurf.order_v(), + psurf.basis_u().begin(), psurf.basis_v().begin(), + space_coefs.begin(), 3); + space_psurf.writeStandardHeader(os); + space_psurf.write(os); +} + //=========================================================================== void SplineDebugUtils::writeTrimmedInfo(BoundedSurface& bd_sf, std::ostream& os, double z) diff --git a/gotools-core/src/utils/BoundingBox.C b/gotools-core/src/utils/BoundingBox.C index e82294500..be7b96c83 100644 --- a/gotools-core/src/utils/BoundingBox.C +++ b/gotools-core/src/utils/BoundingBox.C @@ -148,7 +148,7 @@ bool BoundingBox::containsBox(const BoundingBox& box, double tol) const void BoundingBox::addUnionWith(const Point& pt) //=========================================================================== { - ALWAYS_ERROR_IF (low_.dimension() != pt.dimension(), + ALWAYS_ERROR_IF (valid_ && low_.dimension() != pt.dimension(), "Dimension mismatch."); if (!valid_) { diff --git a/gotools-core/src/utils/QRFactorization.C b/gotools-core/src/utils/QRFactorization.C new file mode 100644 index 000000000..1dd18fe97 --- /dev/null +++ b/gotools-core/src/utils/QRFactorization.C @@ -0,0 +1,144 @@ +/* + * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, + * Applied Mathematics, Norway. + * + * Contact information: E-mail: tor.dokken@sintef.no + * SINTEF ICT, Department of Applied Mathematics, + * P.O. Box 124 Blindern, + * 0314 Oslo, Norway. + * + * This file is part of GoTools. + * + * GoTools is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * GoTools is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public + * License along with GoTools. If not, see + * . + * + * In accordance with Section 7(b) of the GNU Affero General Public + * License, a covered work must retain the producer line in every data + * file that is created or manipulated using GoTools. + * + * Other Usage + * You can be released from the requirements of the license by purchasing + * a commercial license. Buying such a license is mandatory as soon as you + * develop commercial activities involving the GoTools library without + * disclosing the source code of your own applications. + * + * This file may be used in accordance with the terms contained in a + * written agreement between you and SINTEF ICT. + */ + +#include "GoTools/utils/QRFactorization.h" +#include "GoTools/utils/errormacros.h" +#include +#include + +using namespace Go; +using std::vector; + +//============================================================================== +void QRFactorization::QRDecomp(vector& A, int n, int m, + vector& Q, vector& R) +//============================================================================== +{ + if ((int)A.size() != n*m) + THROW("QRFactorization::QRDecomp: Expecting matrix size n*m"); + + double eps = 1.0e-12; + int kt = std::min(n,m); + + R = A; + Q.resize(m*m, 0.0); + for (int kr=0; kr v0(m-kr, 0.0); + for (int ki=kr; ki= 0.0) ? 1.0 : -1.0; + double sgn = (fabs(v0[0]-normv0) < 0.1*normv0) ? 1.0 : -1.0; + + vector v(v0.begin(), v0.end()); + v[0] += sgn*normv0; + double normv = 0.0; + for (size_t i=0; i eps) + { + for (size_t i=0; i& Q, vector& R, int n, int m, + vector& b, int dim, vector& x) +//============================================================================== +{ + if ((int)Q.size() != m*m) + THROW("QRFactorization::QRSolve: Expecting matrix size m*m"); + if ((int)b.size() != dim*m) + THROW("QRFactorization::QRSolve: Expecting left hand size dim*m"); + + x.resize(dim*n, 0.0); + for (int kb=0; kb b2(n, 0.0); + for (int kj=0; kj=0; --kr) + { + x[kb*n+kr] = b2[kr]; + double val = 0.0; + for (int ki=kr+1; kitranslateCurve(dir); + } + // inherited from ParamCurve virtual SplineCurve* geometryCurve();