40 label NURBS3DCurve::sgn(
const scalar val)
const 42 return val>=scalar(0) ? 1 : -1;
46 scalar NURBS3DCurve::abs(
const scalar val)
const 48 return (sgn(val) == 1) ? val : -val;
52 label NURBS3DCurve::mod(
const label
x,
const label interval)
const 54 label ratio(
x%interval);
55 return ratio < 0 ? ratio + interval : ratio;
59 void NURBS3DCurve::setUniformU()
62 label nPts(curve.size());
66 u_[ptI] = scalar(ptI)/scalar(nPts - 1);
71 bool NURBS3DCurve::bound
74 const scalar minVal = 1
e-7,
75 const scalar maxVal = 0.999999
98 void 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);
211 const NURBSbasis& basis,
212 const List<vector>& CPs,
213 const List<scalar>& weights,
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());
413 List<vector> invertedCPs(nCPs,
Zero);
414 List<scalar> invertedWeights(nCPs,
Zero);
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)
625 curveNrm.normalise();
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] <<
" ";
1011 const fileName dirName,
1012 const fileName fileName
1017 OFstream curveFile(dirName/fileName);
1018 OFstream curveFileCPs(dirName/fileName +
"CPs");
1019 OFstream curveFileBases(dirName/fileName +
"Bases");
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] <<
" ";
1118 const fileName dirName,
1119 const fileName fileName
1124 OFstream curveFile(dirName/fileName);
1125 OFstream curveFileCPs(dirName/fileName +
"CPs");
1126 OFstream curveFileBases(dirName/fileName +
"Bases");
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;
void size(const label n)
Older name for setAddressableSize.
A class for handling file names.
NURBS3DCurve(const NURBSbasis &basis, const List< vector > &CPs, const List< scalar > &weights, const scalarField &u, const label nPts, const word name="NURBS3DCurve")
Construct from control points, weights and basis order and parametric coordinates.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
vector curvePoint(const scalar u) const
Curve point cartesian coordinates at ui.
Output to file stream, using an OSstream.
Ostream & endl(Ostream &os)
Add newline and flush stream.
scalar findClosestCurvePoint(const vector &targetPoint, const label maxIter=1000, const scalar tolerance=1.e-13)
Find curve point which is closest to given point using Newton-Raphson. Returns the param coordinate...
scalar length() const
Calculate length for the entire curve length.
#define forAll(list, i)
Loop across all elements in list.
void buildCurve()
Build curve.
Type & operator[](const label i)
Return element of UList.
const label & nCPs() const
void flipNrmOrientation()
Flip the orientation of the nrm.
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
vector curveDerivativeU(const scalar u) const
Curve derivative wrt u at point ui.
vector curveDerivativeUU(const scalar u) const
Curve second derivative wrt u at point ui.
void setNrm3DOrientation(const vector &givenNrm, const vector &givenTan)
Take a given normal and use to determine if NURBS normals should be reversed. Computation taken from ...
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
bool checkRange(const scalar u, const label CPI, const label degree) const
Check if given parametric coordinate u and CP are linked.
const vector nrm2D(const scalar zVal, const scalar u) const
List< scalar > scalarList
A List of scalars.
void setCPs(const List< vector > &CPs)
Set CPs.
NURBSbasis function. Used to construct NURBS curves, surfaces and volumes.
scalarList genEquidistant(const label nPts=100, const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
Generate points on the curve which are evenly spaced in cartesian coordinate distances.
const label & degree() const
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 invert()
Invert CPs order and re-build curve.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
scalar basisDerivativeUU(const label iCP, const label degree, const scalar u) const
Basis second derivative w.r.t u.
scalar basisValue(const label iCP, const label degree, const scalar u) const
Basis value.
scalar basisDerivativeU(const label iCP, const label degree, const scalar u) const
Basis derivative w.r.t u.
vector curveDerivativeWeight(const scalar u, const label CPI)
Curve derivative wrt CPII at point u.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
bool checkRange(const scalar u, const label CPI, const label degree) const
Checks to see if given u is affected by given CP.
void setNrm2DOrientation(const vector &givenNrm, const scalar zVal)
scalar lengthDerivativeU(const scalar uStart, const scalar uEnd, const label nPts)
Calculate length from starting to ending indices via computational evaluation using trapezoid rule...
static bool master(const label communicator=worldComm)
Am I the master rank.
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.
messageStream Info
Information stream (stdout output on master, null elsewhere)
label size() const noexcept
The number of elements in the UList.
Field< vector > vectorField
Specialisation of Field<T> for vector.
void insertKnot(const scalarField &oldKnots, const scalar uBar, const label kInsert)
Insert a knot by re-computing the control points.
void setName(const word &name)
Set name.
dimensionedScalar tan(const dimensionedScalar &ds)
scalar curveDerivativeCP(const scalar u, const label CPI)
Curve derivative wrt b at point ui: scalar since dx/dX = dy/dY = dz/dZ.
void write()
Write curve to file.
void setWeights(const List< scalar > &weights)
Set weights.
static constexpr const zero Zero
Global zero (0)