40label NURBS3DCurve::sgn(
const scalar val)
const
42 return val>=scalar(0) ? 1 : -1;
46scalar NURBS3DCurve::abs(
const scalar val)
const
48 return (sgn(val) == 1) ? val : -val;
52label NURBS3DCurve::mod(
const label
x,
const label interval)
const
54 label ratio(
x%interval);
55 return ratio < 0 ? ratio + interval : ratio;
59void NURBS3DCurve::setUniformU()
62 label nPts(curve.size());
66 u_[ptI] = scalar(ptI)/scalar(nPts - 1);
74 const scalar minVal = 1
e-7,
75 const scalar maxVal = 0.999999
98void NURBS3DCurve::setEquidistantU
101 const label lenAcc = 25,
102 const label maxIter = 10,
103 const label spacingCorrInterval=-1,
104 const scalar tolerance = 1.e-5
107 const label nPts(
U.size());
108 const scalar xLength(
length() /(nPts - 1));
109 const scalar uLength(scalar(1) / scalar(nPts - 1));
112 U[nPts - 1] = scalar(1);
114 for (label ptI = 1; ptI<(nPts - 1); ptI++)
116 const scalar UPrev(
U[ptI - 1]);
117 scalar& UCurr(
U[ptI]);
118 scalar direc(scalar(1));
119 scalar xDiff(scalar(0));
120 scalar
delta(scalar(0));
121 bool overReached(
false);
123 UCurr = UPrev + uLength;
128 bool bounded(bound(UCurr));
131 xDiff = xLength -
delta;
134 if (abs(xDiff) < tolerance)
142 if (bounded && (direc == 1))
151 else if (direc == scalar(1))
166 while (iter < maxIter)
170 UCurr += direc * uLength;
176 (spacingCorrInterval != -1)
177 && (mod(ptI, spacingCorrInterval) == 0)
181 xDiff = (ptI * xLength) -
delta;
186 xDiff = xLength -
delta;
190 if (abs(xDiff) < tolerance)
196 const scalar oldDirec(direc);
197 direc = sgn(xDiff) * abs(oldDirec);
228 nrmOrientation_(ALIGNED)
253 nrmOrientation_(ALIGNED)
271 weights_(CPs.size(), scalar(1)),
296 givenInitNrm_ = givenNrm;
301 if ((givenNrm & curveNrm) >= scalar(0))
303 nrmOrientation_ = ALIGNED;
307 nrmOrientation_ = OPPOSED;
310 Info<<
"Initial nrmOrientation after comparison to NURBS u = 0 nrm : "
322 givenInitNrm_ = givenNrm;
327 curveNrm.
x() = -
tan.y();
328 curveNrm.
y() =
tan.x();
331 if ((givenNrm & curveNrm) >= scalar(0))
333 nrmOrientation_ = ALIGNED;
337 nrmOrientation_ = OPPOSED;
340 Info<<
"Initial nrmOrientation after comparison to NURBS u = 0 nrm : "
348 if (nrmOrientation_ == ALIGNED)
350 nrmOrientation_ = OPPOSED;
354 nrmOrientation_ = ALIGNED;
379 const label degree(basis_.
degree());
386 const scalar u(u_[ptI]);
393 denom += basis_.
basisValue(CPJ, degree, u)*weights_[CPJ];
402 * weights_[CPI]/denom;
410 Info<<
"Inverting NURBS curve " << name_ <<
endl;
412 const label nCPs(CPs_.
size());
416 for (label CPI = 0; CPI<nCPs; CPI++)
418 invertedCPs[CPI] = CPs_[nCPs - 1 - CPI];
419 invertedWeights[CPI] = weights_[nCPs - 1 - CPI];
423 weights_ = invertedWeights;
433 const label spacingCorrInterval,
434 const scalar tolerance
460 const label degree(basis_.
degree());
465 NW += basis_.
basisValue(CPI, degree, u)*weights_[CPI];
485 const vector& targetPoint,
487 const scalar tolerance
497 const scalar distLoc(
mag(targetPoint -
curve[ptI]));
507 scalar u(scalar(closeI)/scalar(this->
size()-1));
515 scalar lhs((dxdu&dxdu) + ((xu - targetPoint) & d2xdu2));
516 scalar rhs(-((xu - targetPoint) & dxdu));
532 while ((iter++< maxIter) && (res > tolerance));
545 <<
"Finding curve point closest to " << targetPoint <<
" failed."
555 const vector& targetPoint,
556 const scalar initGuess,
558 const scalar tolerance
571 scalar lhs((dxdu&dxdu) + ((xu - targetPoint) & d2xdu2));
572 scalar rhs(-((xu - targetPoint) & dxdu));
590 while ((iter++< maxIter) && (res > tolerance));
603 <<
"Finding curve point closest to " << targetPoint
616 if (nrmOrientation_ == ALIGNED)
636 curveNrm.
x() = -nrmOrientation_*
tan.y();
637 curveNrm.
y() = nrmOrientation_*
tan.x();
639 curveNrm /=
mag(curveNrm);
654 const label degree(basis_.
degree());
655 const label nCPs(basis_.
nCPs());
660 for (label CPI = 0; CPI < (kInsert - degree + 1); CPI++)
662 newCPs[CPI] = CPs_[CPI];
665 for (label CPI = (kInsert - degree + 1); CPI < (kInsert + 1); CPI++)
667 const scalar uIOld(oldKnots[CPI]);
668 const scalar uIDOld(oldKnots[CPI + degree]);
669 const scalar ratio((uBar - uIOld) /(uIDOld - uIOld));
671 newCPs[CPI] = (ratio*CPs_[CPI] + (1 - ratio)*CPs_[CPI - 1]);
674 for (label CPI= (kInsert + 1); CPI<newCPs.
size(); CPI++)
676 newCPs[CPI] = CPs_[CPI - 1];
681 weights_ = newWeights;
692 const label spacingCorrInterval,
693 const scalar tolerance
731 const label degree(basis_.
degree());
739 const label lenSize(uIEnd - uIStart + 1);
749 for (label uI = 0; uI < (lenSize - 1); uI++)
753 *(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))
754 *(u_[uIStart + uI + 1]-u_[uIStart + uI]);
775 localU[uI] = uStart + scalar(uI)/scalar(nPts - 1)*(uEnd - uStart);
780 for (label uI = 0; uI < (nPts - 1); uI++)
784 *(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))
785 *(localU[uI + 1]-localU[uI]);
802 const label degree(basis_.
degree());
810 const scalar basisValue(basis_.
basisValue(CPI, degree, u));
813 NWP += basisValue * weights_[CPI] * CPs_[CPI];
814 dNduWP += basisDeriv * weights_[CPI] * CPs_[CPI];
815 NW += basisValue * weights_[CPI];
816 dNduW += basisDeriv * weights_[CPI];
819 const vector uDerivative((dNduWP - NWP*dNduW/NW)/NW);
827 const label degree(basis_.
degree());
833 scalar d2Ndu2W(
Zero);
837 const scalar basisValue(basis_.
basisValue(CPI, degree, u));
841 NWP += basisValue * weights_[CPI] * CPs_[CPI];
842 dNduWP += basisDeriv * weights_[CPI] * CPs_[CPI];
843 d2Ndu2WP += basis2Deriv * weights_[CPI] * CPs_[CPI];
844 NW += basisValue * weights_[CPI];
845 dNduW += basisDeriv * weights_[CPI];
846 d2Ndu2W += basis2Deriv * weights_[CPI];
853 - scalar(2)*dNduWP*dNduW/NW
855 + scalar(2)*NWP*dNduW*dNduW/NW/NW
870 const label degree(basis_.
degree());
875 NW += basis_.
basisValue(CPJ, degree, u) * weights_[CPJ];
878 const scalar basisValueI(basis_.
basisValue(CPI, degree, u));
879 const scalar CPDerivative(basisValueI * weights_[CPI] / NW);
892 const label degree(basis_.
degree());
898 const scalar basisValue(basis_.
basisValue(CPJ, degree, u));
899 NWP += basisValue * weights_[CPJ] * CPs_[CPJ];
900 NW += basisValue * weights_[CPJ];
904 const scalar basisValueI(basis_.
basisValue(CPI, degree, u));
905 const vector WDerivative(basisValueI/NW * (CPs_[CPI] - NWP/NW));
922 scalar lDerivative(
Zero);
926 scalar& uLocal(localU[uI]);
927 uLocal = uStart + scalar(uI)/scalar(nPts - 1)*(uEnd - uStart);
933 for (label uI = 0; uI<(nPts - 1); uI++)
938 (dxdu[uI + 1]&d2xdu2[uI + 1])/
mag(dxdu[uI + 1])
939 + (dxdu[uI]&d2xdu2[uI])/
mag(dxdu[uI])
941 * (localU[uI + 1]-localU[uI]);
970 label degree(basis_.
degree());
976 curveFile <<
field[uI].x() <<
" "
977 <<
field[uI].y() <<
" "
984 curveFileCPs << CPs_[CPI].x() <<
" "
985 << CPs_[CPI].y() <<
" "
992 const scalar u(u_[uI]);
995 curveFileBases << u <<
" ";
999 basesValues[CPI] = basis_.
basisValue(CPI, degree, u);
1000 curveFileBases << basesValues[CPI] <<
" ";
1003 curveFileBases <<
endl;
1020 label degree(basis_.
degree());
1026 curveFile <<
field[uI].x() <<
" "
1027 <<
field[uI].y() <<
" "
1034 curveFileCPs << CPs_[CPI].x() <<
" "
1035 << CPs_[CPI].y() <<
" "
1042 const scalar u(u_[uI]);
1045 curveFileBases << u <<
" ";
1049 basesValues[CPI] = basis_.
basisValue(CPI, degree, u);
1050 curveFileBases << basesValues[CPI] <<
" ";
1053 curveFileBases <<
endl;
1075 label degree(basis_.
degree());
1082 <<
field[uI].x() <<
" "
1083 <<
field[uI].y() <<
" "
1084 <<
field[uI].z() <<
")"
1091 << CPs_[CPI].x() <<
" "
1092 << CPs_[CPI].y() <<
" "
1093 << CPs_[CPI].z() <<
")"
1099 const scalar u(u_[uI]);
1102 curveFileBases << u <<
" ";
1106 basesValues[CPI] = basis_.
basisValue(CPI, degree, u);
1107 curveFileBases << basesValues[CPI] <<
" ";
1110 curveFileBases <<
endl;
1127 label degree(basis_.
degree());
1134 <<
field[uI].x() <<
" "
1135 <<
field[uI].y() <<
" "
1136 <<
field[uI].z() <<
")"
1143 << CPs_[CPI].x() <<
" "
1144 << CPs_[CPI].y() <<
" "
1145 << CPs_[CPI].z() <<
")"
1151 const scalar u(u_[uI]);
1154 curveFileBases << u <<
" ";
1158 basesValues[CPI] = basis_.
basisValue(CPI, degree, u);
1159 curveFileBases << basesValues[CPI] <<
" ";
1162 curveFileBases <<
endl;
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
const vector nrm3D(const vector &refTan, const scalar u) const
Find the normal to the curve, with the option of forcing a z-plane.
void insertKnot(const scalarField &oldKnots, const scalar uBar, const label kInsert)
Insert a knot by re-computing the control points.
void setNrm3DOrientation(const vector &givenNrm, const vector &givenTan)
vector curveDerivativeUU(const scalar u) const
Curve second derivative wrt u at point ui.
void setCPs(const List< vector > &CPs)
Set CPs.
void makeEquidistant(const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
Make curve points equidistant in cartesian space.
vector curveDerivativeWeight(const scalar u, const label CPI)
Curve derivative wrt CPII at point u.
scalar length() const
Calculate length for the entire curve length.
void setName(const word &name)
Set name.
void invert()
Invert CPs order and re-build curve.
scalar findClosestCurvePoint(const vector &targetPoint, const label maxIter=1000, const scalar tolerance=1.e-13)
void buildCurve()
Build curve.
scalarList genEquidistant(const label nPts=100, const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
const vector nrm2D(const scalar zVal, const scalar u) const
scalar lengthDerivativeU(const scalar uStart, const scalar uEnd, const label nPts)
void write()
Write curve to file.
bool checkRange(const scalar u, const label CPI, const label degree) const
Check if given parametric coordinate u and CP are linked.
vector curveDerivativeU(const scalar u) const
Curve derivative wrt u at point ui.
void setNrm2DOrientation(const vector &givenNrm, const scalar zVal)
void flipNrmOrientation()
Flip the orientation of the nrm.
void setWeights(const List< scalar > &weights)
Set weights.
vector curvePoint(const scalar u) const
Curve point cartesian coordinates at ui.
scalar curveDerivativeCP(const scalar u, const label CPI)
NURBSbasis function. Used to construct NURBS curves, surfaces and volumes.
scalar basisDerivativeUU(const label iCP, const label degree, const scalar u) const
Basis second derivative w.r.t u.
scalar basisDerivativeU(const label iCP, const label degree, const scalar u) const
Basis derivative w.r.t u.
const label & nCPs() const
const label & degree() const
scalar basisValue(const label iCP, const label degree, const scalar u) const
Basis value.
bool checkRange(const scalar u, const label CPI, const label degree) const
Checks to see if given u is affected by given CP.
Output to file stream, using an OSstream.
label size() const noexcept
The number of elements in the UList.
void size(const label n)
Older name for setAddressableSize.
T & operator[](const label i)
Return element of UList.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
A single curve in a graph.
A class for handling file names.
splitCell * master() const
A class for handling words, derived from Foam::string.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensionedScalar tan(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
List< scalar > scalarList
A List of scalars.
vector point
Point is a vector.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.