41label NURBS3DSurface::sgn(
const scalar val)
const
43 return val >= scalar(0) ? 1 : -1;
47scalar NURBS3DSurface::abs(
const scalar val)
const
49 return (sgn(val) == 1)? val: -val;
53label NURBS3DSurface::mod(
const label
x,
const label interval)
const
55 label ratio(
x%interval);
56 return ratio < 0 ? ratio+interval : ratio;
60void NURBS3DSurface::setCPUVLinking()
62 const label uNCPs(uBasis_.
nCPs());
63 const label vNCPs(vBasis_.
nCPs());
65 CPsUCPIs_.
setSize(uNCPs*vNCPs, -1);
66 CPsVCPIs_.
setSize(uNCPs*vNCPs, -1);
68 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
70 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
72 const label CPI(vCPI*uNCPs + uCPI);
73 CPsUCPIs_[CPI] = uCPI;
74 CPsVCPIs_[CPI] = vCPI;
80void NURBS3DSurface::setUniformUV
89 v.setSize(nUPts*nVPts,
Zero);
91 for (label uI = 0; uI<nUPts; uI++)
93 scalar uValue = scalar(uI)/scalar(nUPts - 1);
94 for (label vI = 0; vI<nVPts; vI++)
96 const label ptI(uI*nVPts + vI);
98 v[ptI] = scalar(vI)/scalar(nVPts - 1);
104void NURBS3DSurface::setUniformUV()
106 setUniformUV(u_, v_, nUPts_, nVPts_);
119 boundDirection(u, minVal, maxVal)
120 || boundDirection(v, minVal, maxVal);
126bool NURBS3DSurface::boundDirection
133 bool boundPoint(
false);
140 else if (u > scalar(1))
150void NURBS3DSurface::setEquidistantR
155 const label lenAcc = 25,
156 const label maxIter = 10,
157 const label spacingCorrInterval = -1,
158 const scalar tolerance = 1.e-5
161 const label nPts(
R.size());
163 scalar xLength(
Zero);
164 const scalar rLength(scalar(1) / scalar(nPts - 1));
166 if (paramR == PARAMU)
168 xLength =
lengthU(SHeld) / (nPts - 1);
172 xLength =
lengthV(SHeld) / (nPts - 1);
176 R[nPts-1] = scalar(1);
178 for (label ptI=1; ptI<(nPts - 1); ptI++)
180 const scalar& RPrev(
R[ptI - 1]);
181 scalar& RCurr(
R[ptI]);
185 bool overReached(
false);
187 RCurr = RPrev + rLength;
194 if (paramR == PARAMU)
196 bounded = bound(RCurr, SNull);
201 bounded = bound(SNull, RCurr);
205 xDiff = xLength -
delta;
208 if (abs(xDiff) < tolerance)
216 if (bounded && (direc == 1))
225 else if (direc == scalar(1))
240 while (iter < maxIter)
244 RCurr += direc * rLength;
246 if (paramR == PARAMU)
258 (spacingCorrInterval != -1)
259 && (mod(ptI, spacingCorrInterval) == 0)
262 if (paramR == PARAMU)
271 xDiff = (ptI * xLength) -
delta;
275 if (paramR == PARAMU)
284 xDiff = xLength -
delta;
288 if (abs(xDiff) < tolerance)
294 const scalar oldDirec(direc);
295 direc = sgn(xDiff) * abs(oldDirec);
320 u_(nUPts*nVPts,
Zero),
321 v_(nUPts*nVPts,
Zero),
322 weights_(CPs.size(), scalar(1)),
334 nrmOrientation_(ALIGNED),
336 boundaryCPIDs_(nullptr)
358 u_(nUPts*nVPts,
Zero),
359 v_(nUPts*nVPts,
Zero),
372 nrmOrientation_(ALIGNED),
374 boundaryCPIDs_(nullptr)
397 u_(nUPts*nVPts,
Zero),
398 v_(nUPts*nVPts,
Zero),
399 weights_(CPs.size(), scalar(1)),
403 uBasis_(nCPsU, uDegree),
404 vBasis_(nCPsV, vDegree),
411 nrmOrientation_(ALIGNED),
413 boundaryCPIDs_(nullptr)
416 if (nCPsU*nCPsV != CPs_.
size())
419 <<
"nCPsU*nCPsV " << nCPsU*nCPsV
420 <<
" not equal to size of CPs " << CPs_.
size()
447 u_(nUPts*nVPts,
Zero),
448 v_(nUPts*nVPts,
Zero),
449 weights_(CPs.size(), scalar(1)),
453 uBasis_(nCPsU, uDegree, knotsU),
454 vBasis_(nCPsV, vDegree, knotsV),
461 nrmOrientation_(ALIGNED),
463 boundaryCPIDs_(nullptr)
466 if (nCPsU*nCPsV != CPs_.
size())
469 <<
"nCPsU*nCPsV " << nCPsU*nCPsV
470 <<
" not equal to size of CPs " << CPs_.
size()
496 u_(nUPts*nVPts,
Zero),
497 v_(nUPts*nVPts,
Zero),
502 uBasis_(nCPsU, uDegree),
503 vBasis_(nCPsV, vDegree),
510 nrmOrientation_(ALIGNED),
512 boundaryCPIDs_(nullptr)
515 if (nCPsU*nCPsV != CPs_.
size())
518 <<
"nCPsU*nCPsV " << nCPsU*nCPsV
519 <<
" not equal to size of CPs " << CPs_.
size()
548 u_(nUPts*nVPts,
Zero),
549 v_(nUPts*nVPts,
Zero),
554 uBasis_(nCPsU, uDegree, knotsU),
555 vBasis_(nCPsV, vDegree, knotsV),
562 nrmOrientation_(ALIGNED),
564 boundaryCPIDs_(nullptr)
567 if (nCPsU*nCPsV != CPs_.
size())
570 <<
"nCPsU*nCPsV " << nCPsU*nCPsV
571 <<
" not equal to size of CPs " << CPs_.
size()
593 givenInitNrm_ = givenNrm;
594 surfaceNrm /=
mag(surfaceNrm);
596 const scalar relation(givenNrm & surfaceNrm);
600 nrmOrientation_ = ALIGNED;
604 nrmOrientation_ = OPPOSED;
607 Info<<
"Initial nrmOrientation after comparison to NURBS u="
608 << u <<
",v=" << v <<
" nrm: " << nrmOrientation_
615 if (nrmOrientation_ == ALIGNED)
617 nrmOrientation_ = OPPOSED;
621 nrmOrientation_ = ALIGNED;
646 const label uDegree(uBasis_.
degree());
647 const label vDegree(vBasis_.
degree());
648 const label uNCPs(uBasis_.
nCPs());
649 const label vNCPs(vBasis_.
nCPs());
660 for (label uI = 0; uI<nUPts_; uI++)
662 for (label vI = 0; vI<nVPts_; vI++)
664 const label ptI(uI*nVPts_ + vI);
665 const scalar& u(u_[ptI]);
666 const scalar& v(v_[ptI]);
670 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
672 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
674 const label CPI(vCPI*uNCPs + uCPI);
684 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
686 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
688 const label CPI(vCPI*uNCPs + uCPI);
705 Info<<
"Inverting NURBS surface " << name_ <<
" in u." <<
endl;
707 const label uNCPs(uBasis_.
nCPs());
708 const label vNCPs(vBasis_.
nCPs());
712 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
714 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
716 const label CPI(vCPI*uNCPs + uCPI);
717 const label invUCPI(uNCPs-1-uCPI);
718 const label uInvCPI(vCPI*uNCPs + invUCPI);
720 invertedCPs[CPI] = CPs_[uInvCPI];
721 invertedWeights[CPI] = weights_[uInvCPI];
726 weights_ = invertedWeights;
734 Info<<
"Inverting NURBS surface " << name_ <<
" in v." <<
endl;
736 const label uNCPs(uBasis_.
nCPs());
737 const label vNCPs(vBasis_.
nCPs());
741 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
743 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
745 const label CPI(vCPI*uNCPs + uCPI);
746 const label invVCPI(vNCPs-1-vCPI);
747 const label vInvCPI(invVCPI*uNCPs + uCPI);
749 invertedCPs[CPI] = CPs_[vInvCPI];
750 invertedWeights[CPI] = weights_[vInvCPI];
755 weights_ = invertedWeights;
763 Info<<
"Inverting NURBS surface " << name_ <<
" in u and v." <<
endl;
765 const label uNCPs(uBasis_.
nCPs());
766 const label vNCPs(vBasis_.
nCPs());
770 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
772 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
774 const label CPI(vCPI*uNCPs + uCPI);
775 const label invUCPI(uNCPs - 1 - uCPI);
776 const label invVCPI(vNCPs - 1 - vCPI);
777 const label uvInvCPI(invVCPI*uNCPs + invUCPI);
779 invertedCPs[CPI] = CPs_[uvInvCPI];
780 invertedWeights[CPI] = weights_[uvInvCPI];
785 weights_ = invertedWeights;
795 const label spacingCorrInterval,
796 const scalar tolerance
805 for (label vI = 0; vI<nVPts_; vI++)
808 const scalar VHeld(v_[vI]);
814 const label ptI(uI*nVPts_ + vI);
815 uAddressing[uI] = ptI;
832 const label& uAddress(uAddressing[uI]);
833 u_[uAddress] = UI[uI];
838 for (label uI = 0; uI<nUPts_; uI++)
841 const scalar UHeld(u_[uI*nVPts_]);
847 const label ptI(uI*nVPts_ + vI);
848 vAddressing[vI] = ptI;
866 const label& vAddress(vAddressing[vI]);
867 v_[vAddress] = VI[vI];
883 const label uDegree(uBasis_.
degree());
884 const label vDegree(vBasis_.
degree());
885 const label uNCPs(uBasis_.
nCPs());
886 const label vNCPs(vBasis_.
nCPs());
890 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
892 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
894 const label CPI(vCPI*uNCPs + uCPI);
906 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
908 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
910 const label CPI(vCPI*uNCPs + uCPI);
926 const vector& targetPoint,
928 const scalar tolerance
938 const scalar distLoc(
mag(targetPoint-surface[ptI]));
948 scalar u(u_[closePtI]);
949 scalar v(v_[closePtI]);
952 scalar resOld(GREAT);
953 scalar resDeriv(GREAT);
987 const scalar a((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
988 const scalar
b((dxdu&dxdv) + ((xuv-targetPoint) & d2xduv));
990 const scalar d((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
991 const scalar invDenom = 1./(a*d-
b*c);
993 const scalar uRHS(-((xuv-targetPoint) & dxdu));
994 const scalar vRHS(-((xuv-targetPoint) & dxdv));
999 u += ( d*uRHS-
b*vRHS)*invDenom;
1000 v += (-c*uRHS+a*vRHS)*invDenom;
1002 if (boundDirection(u))
1006 if (boundDirection(v))
1016 resDeriv =
mag(res-resOld)/resOld;
1021 else if (nBoundsV >= 5)
1024 resDeriv =
mag(res-resOld)/resOld;
1028 else if (nBoundsU <= 5 && nBoundsV <= 5)
1032 resDeriv =
mag(res-resOld)/resOld;
1038 <<
"More than 5 bounds in both the u and v directions!"
1039 <<
"Something seems weird" << nBoundsU <<
" " << nBoundsV
1043 resDeriv =
mag(res-resOld)/resOld;
1047 while ((iter++ < maxIter) && (res > tolerance));
1050 closestParameters[1] = v;
1055 <<
"Finding surface point closest to " << targetPoint
1056 <<
" for surface " << name_ <<
" failed \n"
1057 <<
" Number of bounding operations in u,v "
1058 << nBoundsU <<
" " << nBoundsV <<
endl
1059 <<
" Residual value and derivative " << res <<
" " << resDeriv
1061 closestParameters = -1;
1064 return closestParameters;
1071 const label maxIter,
1072 const scalar tolerance
1076 auto& paramCoors = tparamCoors.ref();
1079 label nBoundedPoints(0);
1080 scalar maxResidual(0);
1081 scalar maxResidualDeriv(0);
1084 const vector& targetPoint(targetPoints[pI]);
1094 const scalar distLoc(
mag(targetPoint - surface[ptI]));
1104 scalar u(u_[closePtI]);
1105 scalar v(v_[closePtI]);
1108 scalar resOld(GREAT);
1109 scalar resDeriv(GREAT);
1120 const scalar a((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
1121 const scalar
b((dxdu&dxdv) + ((xuv-targetPoint) & d2xduv));
1123 const scalar d((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
1124 const scalar invDenom = 1./(a*d-
b*c);
1126 const scalar uRHS(-((xuv-targetPoint) & dxdu));
1127 const scalar vRHS(-((xuv-targetPoint) & dxdv));
1132 u += ( d*uRHS-
b*vRHS)*invDenom;
1133 v += (-c*uRHS+a*vRHS)*invDenom;
1135 if (boundDirection(u))
1139 if (boundDirection(v))
1149 resDeriv =
mag(res-resOld)/resOld;
1153 else if (nBoundsV >= 5)
1156 resDeriv =
mag(res-resOld)/resOld;
1159 else if (nBoundsU<=5 && nBoundsV<=5)
1163 resDeriv =
mag(res-resOld)/resOld;
1169 <<
"More than 5 bounding operations in both the u and v directions!"
1170 <<
"u direction " << nBoundsU <<
endl
1171 <<
"v direction " << nBoundsV <<
endl
1176 resDeriv =
mag(res-resOld)/resOld;
1180 while ((iter++ < maxIter) && (res > tolerance));
1186 maxResidual =
max(maxResidual, res);
1187 maxResidualDeriv =
max(maxResidualDeriv, resDeriv);
1190 paramCoors[pI].x() = u;
1191 paramCoors[pI].y() = v;
1196 Info<<
"Points that couldn't reach the residual limit:: "
1197 << nBoundedPoints <<
endl
1198 <<
"Max residual after reaching iterations limit:: "
1199 << maxResidual <<
endl
1200 <<
"Max residual derivative after reaching iterations limit:: "
1201 << maxResidualDeriv <<
endl
1210 const vector& targetPoint,
1211 const scalar& uInitGuess,
1212 const scalar& vInitGuess,
1213 const label maxIter,
1214 const scalar tolerance
1219 scalar u(uInitGuess);
1220 scalar v(vInitGuess);
1230 const scalar uLHS((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
1231 const scalar uRHS(-((xuv-targetPoint) & dxdu));
1232 const scalar vLHS((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
1233 const scalar vRHS(-((xuv-targetPoint) & dxdv));
1238 u += uRHS/(uLHS+SMALL);
1239 v += vRHS/(vLHS+SMALL);
1248 while ((iter++ < maxIter) && (res > tolerance));
1254 <<
"Finding surface point closest to " << targetPoint <<
" failed."
1259 closestParameters[1] = v;
1261 return closestParameters;
1269 if (nrmOrientation_ == ALIGNED)
1278 surfaceNrm /=
mag(surfaceNrm);
1289 const label maxIter,
1290 const label spacingCorrInterval,
1291 const scalar tolerance
1305 setUniformUV(
U, V, nUPts, nVPts);
1308 for (label vI = 0; vI<nVPts; vI++)
1311 const scalar VHeld(V[vI]);
1317 const label ptI(uI*nVPts + vI);
1318 uAddressing[uI] = ptI;
1328 spacingCorrInterval,
1335 const label& uAddress(uAddressing[uI]);
1336 U[uAddress] = UI[uI];
1341 for (label uI = 0; uI<nUPts; uI++)
1344 const scalar UHeld(
U[uI*nVPts]);
1350 const label ptI(uI*nVPts + vI);
1351 vAddressing[vI] = ptI;
1362 spacingCorrInterval,
1369 const label& vAddress(vAddressing[vI]);
1370 V[vAddress] = VI[vI];
1387 const label uCPI(CPsUCPIs_[CPI]);
1399 const label uDegree(uBasis_.
degree());
1412 const label vCPI(CPsVCPIs_[CPI]);
1424 const label vDegree(vBasis_.
degree());
1435 const label uDegree,
1455 const label uDegree(uBasis_.
degree());
1456 const label vDegree(vBasis_.
degree());
1464 const label vIConst,
1465 const label uIStart,
1470 const label uLenSize(uIEnd - uIStart + 1);
1472 scalar uLength(
Zero);
1476 const label ptI((uIStart+uI)*nVPts_ + vIConst);
1477 const label& u(u_[ptI]);
1478 const label& v(v_[ptI]);
1484 for(label uI = 0; uI<(uLenSize - 1); uI++)
1486 const label ptI((uIStart+uI)*nVPts_ + vIConst);
1489 0.5*(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))*(u_[ptI + 1]-u_[ptI]);
1498 const scalar vConst,
1499 const scalar uStart,
1507 scalar uLength(
Zero);
1511 scalar& uLocal(localU[uI]);
1512 uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1517 for(label uI = 0; uI<(nPts - 1); uI++)
1520 0.5*(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))*(localU[uI + 1]-localU[uI]);
1529 return lengthU(vIConst, 0, (nUPts_ - 1));
1535 return lengthU(vConst, scalar(0), scalar(1), 100);
1541 const label uIConst,
1542 const label vIStart,
1547 const label vLenSize(vIEnd - vIStart + 1);
1549 scalar vLength(
Zero);
1553 const label ptI((uIConst)*nVPts_ + (vIStart+vI));
1554 const label& u(u_[ptI]);
1555 const label& v(v_[ptI]);
1561 for(label vI = 0; vI<(vLenSize - 1); vI++)
1563 const label ptI((uIConst)*nVPts_ + (vIStart + vI));
1566 0.5*(
mag(dxdv[vI + 1]) +
mag(dxdv[vI]))*(v_[ptI + 1] - v_[ptI]);
1575 const scalar uConst,
1576 const scalar vStart,
1584 scalar vLength(
Zero);
1588 scalar& vLocal(localV[vI]);
1589 vLocal = vStart + scalar(vI)/scalar(nPts - 1)*(vEnd - vStart);
1594 for(label vI = 0; vI<(nPts - 1); vI++)
1597 0.5*(
mag(dxdv[vI + 1]) +
mag(dxdv[vI]))*(localV[vI + 1]-localV[vI]);
1606 return lengthV(uIConst, 0, (nVPts_ - 1));
1612 return lengthV(uConst, scalar(0), scalar(1), 100);
1624 const label uDegree(uBasis_.
degree());
1625 const label vDegree(vBasis_.
degree());
1626 const label uNCPs(uBasis_.
nCPs());
1627 const label vNCPs(vBasis_.
nCPs());
1631 scalar dNduMW(
Zero);
1637 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1639 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1641 const label CPI(vCPI*uNCPs + uCPI);
1642 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1643 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1644 const scalar uBasisDeriv
1647 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1648 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1649 NMW += uBasisValue * vBasisValue * weights_[CPI];
1650 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1654 const vector uDerivative((dNduMWP - dNduMW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1666 const label uDegree(uBasis_.
degree());
1667 const label vDegree(vBasis_.
degree());
1668 const label uNCPs(uBasis_.
nCPs());
1669 const label vNCPs(vBasis_.
nCPs());
1673 scalar dMdvNW(
Zero);
1679 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1681 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1683 const label CPI(vCPI*uNCPs + uCPI);
1684 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1685 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1686 const scalar vBasisDeriv
1689 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1690 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1691 NMW += uBasisValue * vBasisValue * weights_[CPI];
1692 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1696 const vector vDerivative((dMdvNWP - dMdvNW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1708 const label uDegree(uBasis_.
degree());
1709 const label vDegree(vBasis_.
degree());
1710 const label uNCPs(uBasis_.
nCPs());
1711 const label vNCPs(vBasis_.
nCPs());
1717 scalar dNduMW(
Zero);
1718 scalar dMdvNW(
Zero);
1719 scalar dNMduvW(
Zero);
1725 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1727 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1729 const label CPI(vCPI*uNCPs + uCPI);
1730 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1731 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1732 const scalar uBasisDeriv
1734 const scalar vBasisDeriv
1740 NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1741 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1742 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1743 dNMduvWP += uBasisDeriv * vBasisDeriv * weights_[CPI] * CPs_[CPI];
1744 NMW += vBasisValue * uBasisValue * weights_[CPI];
1745 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1746 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1747 dNMduvW += uBasisDeriv * vBasisDeriv * weights_[CPI];
1751 const vector uvDerivative
1755 - (dNMduvW*NMWP + dMdvNW*dNduMWP + dNduMW*dMdvNWP)/(NMW+SMALL)
1756 + scalar(2)*dNduMW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1760 return uvDerivative;
1770 const label uDegree(uBasis_.
degree());
1771 const label vDegree(vBasis_.
degree());
1772 const label uNCPs(uBasis_.
nCPs());
1773 const label vNCPs(vBasis_.
nCPs());
1778 scalar dNduMW(
Zero);
1779 scalar d2Ndu2MW(
Zero);
1785 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1787 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1789 const label CPI(vCPI*uNCPs + uCPI);
1790 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1791 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1792 const scalar uBasisDeriv
1794 const scalar uBasis2Deriv
1797 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1798 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1799 d2Ndu2MWP += uBasis2Deriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1800 NMW += uBasisValue * vBasisValue * weights_[CPI];
1801 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1802 d2Ndu2MW += uBasis2Deriv * vBasisValue * weights_[CPI];
1806 const vector uuDerivative
1810 - scalar(2)*dNduMW*dNduMWP/(NMW+SMALL)
1811 - d2Ndu2MW*NMWP/(NMW+SMALL)
1812 + scalar(2)*dNduMW*dNduMW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1816 return uuDerivative;
1826 const label uDegree(uBasis_.
degree());
1827 const label vDegree(vBasis_.
degree());
1828 const label uNCPs(uBasis_.
nCPs());
1829 const label vNCPs(vBasis_.
nCPs());
1834 scalar dMdvNW(
Zero);
1835 scalar d2Mdv2NW(
Zero);
1841 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1843 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1845 const label CPI(vCPI*uNCPs + uCPI);
1846 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1847 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1848 const scalar vBasisDeriv
1850 const scalar vBasis2Deriv
1853 NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1854 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1855 d2Mdv2NWP += vBasis2Deriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1856 NMW += vBasisValue * uBasisValue * weights_[CPI];
1857 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1858 d2Mdv2NW += vBasis2Deriv * uBasisValue * weights_[CPI];
1862 const vector vvDerivative
1866 - scalar(2)*dMdvNW*dMdvNWP/(NMW+SMALL)
1867 - d2Mdv2NW*NMWP/(NMW+SMALL)
1868 + scalar(2)*dMdvNW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1872 return vvDerivative;
1885 const label uDegree(uBasis_.
degree());
1886 const label vDegree(vBasis_.
degree());
1887 const label uNCPs(uBasis_.
nCPs());
1888 const label vNCPs(vBasis_.
nCPs());
1889 const label uCPI(CPsUCPIs_[CPI]);
1890 const label vCPI(CPsVCPIs_[CPI]);
1893 for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1895 for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1897 const label CPJ(vCPJ*uNCPs + uCPJ);
1898 const scalar uBasisValue(uBasis_.
basisValue(uCPJ, uDegree, u));
1899 const scalar vBasisValue(vBasis_.
basisValue(vCPJ, vDegree, v));
1901 NMW += vBasisValue * uBasisValue * weights_[CPJ];
1906 const scalar CPDerivative
1914 return CPDerivative;
1926 const label uDegree(uBasis_.
degree());
1927 const label vDegree(vBasis_.
degree());
1928 const label uNCPs(uBasis_.
nCPs());
1929 const label vNCPs(vBasis_.
nCPs());
1930 const label uCPI(CPsUCPIs_[CPI]);
1931 const label vCPI(CPsVCPIs_[CPI]);
1935 for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1937 for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1939 const label CPJ(vCPJ*uNCPs + uCPJ);
1940 const scalar uBasisValue(uBasis_.
basisValue(uCPJ, uDegree, u));
1941 const scalar vBasisValue(vBasis_.
basisValue(vCPJ, vDegree, v));
1943 NMWP += uBasisValue * vBasisValue * weights_[CPJ] * CPs_[CPJ];
1944 NMW += uBasisValue * vBasisValue * weights_[CPJ];
1953 * (CPs_[CPI] - (NMWP / (NMW+SMALL)))
1963 const scalar vConst,
1964 const scalar uStart,
1973 scalar ulDerivative(
Zero);
1977 scalar& uLocal(localU[uI]);
1978 uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1984 for(label uI=0; uI<(nPts-1); uI++)
1989 (dxdu[uI+1]&d2xdu2[uI+1])/(
mag(dxdu[uI+1])+SMALL)
1990 + (dxdu[uI ]&d2xdu2[uI ])/(
mag(dxdu[uI ])+SMALL)
1992 * (localU[uI+1]-localU[uI]);
1995 return ulDerivative;
2001 const scalar uConst,
2002 const scalar vStart,
2011 scalar vlDerivative(
Zero);
2015 scalar& vLocal(localV[vI]);
2016 vLocal = vStart + scalar(vI)/scalar(nPts-1)*(vEnd-vStart);
2022 for(label vI=0; vI<(nPts-1); vI++)
2027 (dxdv[vI+1]&d2xdv2[vI+1])/(
mag(dxdv[vI+1])+SMALL)
2028 + (dxdv[vI ]&d2xdv2[vI ])/(
mag(dxdv[vI ])+SMALL)
2030 * (localV[vI+1]-localV[vI]);
2033 return vlDerivative;
2041 if (!boundaryCPIDs_)
2043 const label uNCPs(uBasis_.
nCPs());
2044 const label vNCPs(vBasis_.
nCPs());
2045 const label nBoundCPs(2*uNCPs+2*vNCPs-4);
2046 boundaryCPIDs_.reset(
new labelList(nBoundCPs, -1));
2047 whichBoundaryCPID_.reset(
new labelList(uNCPs*vNCPs, -1));
2051 for(label vI=0; vI<vNCPs; vI+=vNCPs-1)
2053 for(label uI=0; uI<uNCPs; uI++)
2055 const label CPI(vI*uNCPs + uI);
2056 whichBoundaryCPID_()[CPI] = bID;
2057 boundaryCPIDs_()[bID++] = CPI;
2061 for(label uI=0; uI<uNCPs; uI+=uNCPs-1)
2064 for(label vI=1; vI<vNCPs-1; vI++)
2066 const label CPI(vI*uNCPs + uI);
2067 whichBoundaryCPID_()[CPI] = bID;
2068 boundaryCPIDs_()[bID++] = CPI;
2073 return boundaryCPIDs_();
2085 if (!whichBoundaryCPID_)
2090 return whichBoundaryCPID_()[globalCPI];
2113 << surface[ptI].x() <<
" "
2114 << surface[ptI].y() <<
" "
2122 << CPs_[CPI].x() <<
" "
2123 << CPs_[CPI].y() <<
" "
2170 << surface[ptI].x() <<
" "
2171 << surface[ptI].y() <<
" "
2179 << CPs_[CPI].x() <<
" "
2180 << CPs_[CPI].y() <<
" "
2234 << surface[ptI].x() <<
" "
2235 << surface[ptI].y() <<
" "
2236 << surface[ptI].z() <<
")"
2244 << CPs_[CPI].x() <<
" "
2245 << CPs_[CPI].y() <<
" "
2246 << CPs_[CPI].z() <<
")"
2297 << surface[ptI].x() <<
" "
2298 << surface[ptI].y() <<
" "
2299 << surface[ptI].z() <<
")"
2307 << CPs_[CPI].x() <<
" "
2308 << CPs_[CPI].y() <<
" "
2309 << CPs_[CPI].z() <<
")"
2355 <<
"Do not supply a file extension."
2365 faceList surfaceFaces((nUPts_ - 1)*(nUPts_ - 1),
face(4));
2367 for (label fuI = 0; fuI < (nUPts_ - 1); fuI++)
2369 for (label fvI = 0; fvI < (nVPts_ - 1); fvI++)
2371 const label fI(fuI*(nUPts_ - 1) + fvI);
2372 face& surfaceFace(surfaceFaces[fI]);
2374 surfaceFace[0] = (fuI)*nVPts_ + (fvI);
2375 surfaceFace[1] = (fuI + 1)*nVPts_ + (fvI);
2376 surfaceFace[2] = (fuI + 1)*nVPts_ + (fvI + 1);
2377 surfaceFace[3] = (fuI)*nVPts_ + (fvI + 1);
2387 vtkDirName/vtkFileName,
2394 fileName vtkCPFileName(vtkFileName+
"CPs");
2396 const label uNCPs(uBasis_.
nCPs());
2397 const label vNCPs(vBasis_.
nCPs());
2400 for (label fvCPI=0; fvCPI<(vNCPs-1); fvCPI++)
2402 for (label fuCPI=0; fuCPI<(uNCPs-1); fuCPI++)
2404 const label fCPI(fvCPI*(uNCPs-1) + fuCPI);
2405 face& surfaceCPFace(surfaceCPFaces[fCPI]);
2407 surfaceCPFace[0] = (fvCPI)*uNCPs + (fuCPI);
2408 surfaceCPFace[1] = (fvCPI + 1)*uNCPs + (fuCPI);
2409 surfaceCPFace[2] = (fvCPI + 1)*uNCPs + (fuCPI + 1);
2410 surfaceCPFace[3] = (fvCPI)*uNCPs + (fuCPI + 1);
2418 vtkDirName/vtkCPFileName,
#define R(A, B, C, D, E, F, K, M)
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
bool surfacePoint() const
Part of a surface point pair.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void setSize(const label n)
Alias for resize()
scalar surfaceDerivativeCP(const scalar u, const scalar v, const label CPI) const
Surface derivative wrt the weight of CPI at point u,v.
List< scalarList > genEquidistant(const label nUPts=100, const label nVPts=100, const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
Generate points on the surface which are evenly spaced in cartesian.
const labelList & getBoundaryCPIs()
vector surfaceDerivativeUV(const scalar u, const scalar v) const
Surface second derivative wrt u and v at point u,v.
bool checkRangeU(const scalar u, const label CPI, const label uDegree) const
void setCPs(const List< vector > &CPs)
scalar lengthU(const label vIConst, const label uIStart, const label uIEnd) const
void makeEquidistant(const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
scalarList findClosestSurfacePoint(const vector &targetPoint, const label maxIter=100, const scalar tolerance=1.e-6)
vector surfacePoint(const scalar &u, const scalar &v)
vector surfaceDerivativeV(const scalar u, const scalar v) const
Surface derivative wrt v at point u,v.
void setName(const word &name)
scalar lengthDerivativeV(const scalar uConst, const scalar vStart, const scalar vEnd, const label nPts) const
Surface derivative wrt v length along u=const contour range.
void setNrmOrientation(const vector &givenNrm, const scalar u, const scalar v)
scalar lengthV(const label uIConst, const label vIStart, const label vIEnd) const
const labelList & getBoundaryCPIDs()
Get IDs of boundary control points.
void write()
Write curve to file.
vector surfaceDerivativeUU(const scalar u, const scalar v) const
Surface second derivative wrt u at point u,v.
vector surfaceDerivativeVV(const scalar u, const scalar v) const
Surface second derivative wrt v at point u,v.
const vector nrm(scalar u, scalar v)
bool checkRangeV(const scalar v, const label CPI, const label vDegree) const
vector surfaceDerivativeW(const scalar u, const scalar v, const label CPI) const
Surface derivative wrt WI at point u,v.
vector surfaceDerivativeU(const scalar u, const scalar v) const
Surface derivative wrt u at point u,v.
scalar lengthDerivativeU(const scalar vConst, const scalar uStart, const scalar uEnd, const label nPts) const
Surface derivative wrt u length along v=const contour range.
void setWeights(const scalarList &weights)
void flipNrmOrientation()
Flip the orientation of the nrm.
bool checkRangeUV(const scalar v, const scalar u, const label CPI, const label uDegree, const label vDegree) const
const label & whichBoundaryCPI(const label &globalCPI)
Get the boundary CP ID based on the global CP ID.
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.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
T & operator[](const label i)
Return element of UList.
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
word ext() const
Return file name extension (part after last .)
void writeVTK() const
Write VTK freeSurface mesh.
splitCell * master() const
A surfaceWriter for VTK legacy (.vtk) or XML (.vtp) format.
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
List< label > labelList
A List of labels.
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)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.