41 label NURBS3DSurface::sgn(
const scalar val)
const 43 return val >= scalar(0) ? 1 : -1;
47 scalar NURBS3DSurface::abs(
const scalar val)
const 49 return (sgn(val) == 1)? val: -val;
53 label NURBS3DSurface::mod(
const label
x,
const label interval)
const 55 label ratio(
x%interval);
56 return ratio < 0 ? ratio+interval : ratio;
60 void 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;
80 void 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);
104 void NURBS3DSurface::setUniformUV()
106 setUniformUV(u_, v_, nUPts_, nVPts_);
110 bool NURBS3DSurface::bound
119 boundDirection(u, minVal, maxVal)
120 || boundDirection(v, minVal, maxVal);
126 bool NURBS3DSurface::boundDirection
133 bool boundPoint(
false);
140 else if (u > scalar(1))
150 void 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);
309 const List<vector>& CPs,
312 const NURBSbasis& uBasis,
313 const NURBSbasis& vBasis,
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());
658 field = vector::zero;
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());
709 List<vector> invertedCPs(CPs_.
size(),
Zero);
710 List<scalar> invertedWeights(CPs_.
size(),
Zero);
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;
1194 Info<<
"Points that couldn't reach the residual limit:: " 1196 <<
"Max residual after reaching iterations limit:: " 1198 <<
"Max residual derivative after reaching iterations limit:: " 1208 const vector& targetPoint,
1209 const scalar& uInitGuess,
1210 const scalar& vInitGuess,
1211 const label maxIter,
1212 const scalar tolerance
1217 scalar u(uInitGuess);
1218 scalar v(vInitGuess);
1228 const scalar uLHS((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
1229 const scalar uRHS(-((xuv-targetPoint) & dxdu));
1230 const scalar vLHS((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
1231 const scalar vRHS(-((xuv-targetPoint) & dxdv));
1236 u += uRHS/(uLHS+SMALL);
1237 v += vRHS/(vLHS+SMALL);
1246 while ((iter++ < maxIter) && (res > tolerance));
1252 <<
"Finding surface point closest to " << targetPoint <<
" failed." 1257 closestParameters[1] = v;
1259 return closestParameters;
1267 if (nrmOrientation_ == ALIGNED)
1276 surfaceNrm /=
mag(surfaceNrm);
1287 const label maxIter,
1288 const label spacingCorrInterval,
1289 const scalar tolerance
1303 setUniformUV(
U, V, nUPts, nVPts);
1306 for (label vI = 0; vI<nVPts; vI++)
1309 const scalar VHeld(V[vI]);
1315 const label ptI(uI*nVPts + vI);
1316 uAddressing[uI] = ptI;
1326 spacingCorrInterval,
1333 const label& uAddress(uAddressing[uI]);
1334 U[uAddress] = UI[uI];
1339 for (label uI = 0; uI<nUPts; uI++)
1342 const scalar UHeld(
U[uI*nVPts]);
1348 const label ptI(uI*nVPts + vI);
1349 vAddressing[vI] = ptI;
1360 spacingCorrInterval,
1367 const label& vAddress(vAddressing[vI]);
1368 V[vAddress] = VI[vI];
1385 const label uCPI(CPsUCPIs_[CPI]);
1410 const label vCPI(CPsVCPIs_[CPI]);
1433 const label uDegree,
1453 const label uDegree(uBasis_.
degree());
1462 const label vIConst,
1463 const label uIStart,
1468 const label uLenSize(uIEnd - uIStart + 1);
1470 scalar uLength(
Zero);
1474 const label ptI((uIStart+uI)*nVPts_ + vIConst);
1475 const label& u(u_[ptI]);
1476 const label& v(v_[ptI]);
1482 for (label uI = 0; uI<(uLenSize - 1); uI++)
1484 const label ptI((uIStart+uI)*nVPts_ + vIConst);
1487 0.5*(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))*(u_[ptI + 1]-u_[ptI]);
1496 const scalar vConst,
1497 const scalar uStart,
1505 scalar uLength(
Zero);
1509 scalar& uLocal(localU[uI]);
1510 uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1515 for (label uI = 0; uI<(nPts - 1); uI++)
1518 0.5*(
mag(dxdu[uI + 1]) +
mag(dxdu[uI]))*(localU[uI + 1]-localU[uI]);
1527 return lengthU(vIConst, 0, (nUPts_ - 1));
1533 return lengthU(vConst, scalar(0), scalar(1), 100);
1539 const label uIConst,
1540 const label vIStart,
1545 const label vLenSize(vIEnd - vIStart + 1);
1547 scalar vLength(
Zero);
1551 const label ptI((uIConst)*nVPts_ + (vIStart+vI));
1552 const label& u(u_[ptI]);
1553 const label& v(v_[ptI]);
1559 for (label vI = 0; vI<(vLenSize - 1); vI++)
1561 const label ptI((uIConst)*nVPts_ + (vIStart + vI));
1564 0.5*(
mag(dxdv[vI + 1]) +
mag(dxdv[vI]))*(v_[ptI + 1] - v_[ptI]);
1573 const scalar uConst,
1574 const scalar vStart,
1582 scalar vLength(
Zero);
1586 scalar& vLocal(localV[vI]);
1587 vLocal = vStart + scalar(vI)/scalar(nPts - 1)*(vEnd - vStart);
1592 for (label vI = 0; vI<(nPts - 1); vI++)
1595 0.5*(
mag(dxdv[vI + 1]) +
mag(dxdv[vI]))*(localV[vI + 1]-localV[vI]);
1604 return lengthV(uIConst, 0, (nVPts_ - 1));
1610 return lengthV(uConst, scalar(0), scalar(1), 100);
1622 const label uDegree(uBasis_.
degree());
1623 const label vDegree(vBasis_.
degree());
1624 const label uNCPs(uBasis_.
nCPs());
1625 const label vNCPs(vBasis_.
nCPs());
1629 scalar dNduMW(
Zero);
1635 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1637 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1639 const label CPI(vCPI*uNCPs + uCPI);
1640 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1641 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1642 const scalar uBasisDeriv
1645 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1646 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1647 NMW += uBasisValue * vBasisValue * weights_[CPI];
1648 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1652 const vector uDerivative((dNduMWP - dNduMW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1664 const label uDegree(uBasis_.
degree());
1665 const label vDegree(vBasis_.
degree());
1666 const label uNCPs(uBasis_.
nCPs());
1667 const label vNCPs(vBasis_.
nCPs());
1671 scalar dMdvNW(
Zero);
1677 for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1679 for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1681 const label CPI(vCPI*uNCPs + uCPI);
1682 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1683 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1684 const scalar vBasisDeriv
1687 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1688 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1689 NMW += uBasisValue * vBasisValue * weights_[CPI];
1690 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1694 const vector vDerivative((dMdvNWP - dMdvNW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1706 const label uDegree(uBasis_.
degree());
1707 const label vDegree(vBasis_.
degree());
1708 const label uNCPs(uBasis_.
nCPs());
1709 const label vNCPs(vBasis_.
nCPs());
1715 scalar dNduMW(
Zero);
1716 scalar dMdvNW(
Zero);
1717 scalar dNMduvW(
Zero);
1723 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1725 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1727 const label CPI(vCPI*uNCPs + uCPI);
1728 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1729 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1730 const scalar uBasisDeriv
1732 const scalar vBasisDeriv
1738 NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1739 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1740 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1741 dNMduvWP += uBasisDeriv * vBasisDeriv * weights_[CPI] * CPs_[CPI];
1742 NMW += vBasisValue * uBasisValue * weights_[CPI];
1743 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1744 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1745 dNMduvW += uBasisDeriv * vBasisDeriv * weights_[CPI];
1749 const vector uvDerivative
1753 - (dNMduvW*NMWP + dMdvNW*dNduMWP + dNduMW*dMdvNWP)/(NMW+SMALL)
1754 + scalar(2)*dNduMW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1758 return uvDerivative;
1768 const label uDegree(uBasis_.
degree());
1769 const label vDegree(vBasis_.
degree());
1770 const label uNCPs(uBasis_.
nCPs());
1771 const label vNCPs(vBasis_.
nCPs());
1776 scalar dNduMW(
Zero);
1777 scalar d2Ndu2MW(
Zero);
1783 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1785 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1787 const label CPI(vCPI*uNCPs + uCPI);
1788 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1789 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1790 const scalar uBasisDeriv
1792 const scalar uBasis2Deriv
1795 NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1796 dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1797 d2Ndu2MWP += uBasis2Deriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1798 NMW += uBasisValue * vBasisValue * weights_[CPI];
1799 dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1800 d2Ndu2MW += uBasis2Deriv * vBasisValue * weights_[CPI];
1804 const vector uuDerivative
1808 - scalar(2)*dNduMW*dNduMWP/(NMW+SMALL)
1809 - d2Ndu2MW*NMWP/(NMW+SMALL)
1810 + scalar(2)*dNduMW*dNduMW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1814 return uuDerivative;
1824 const label uDegree(uBasis_.
degree());
1825 const label vDegree(vBasis_.
degree());
1826 const label uNCPs(uBasis_.
nCPs());
1827 const label vNCPs(vBasis_.
nCPs());
1832 scalar dMdvNW(
Zero);
1833 scalar d2Mdv2NW(
Zero);
1839 for (label vCPI=0; vCPI<vNCPs; vCPI++)
1841 for (label uCPI=0; uCPI<uNCPs; uCPI++)
1843 const label CPI(vCPI*uNCPs + uCPI);
1844 const scalar uBasisValue(uBasis_.
basisValue(uCPI, uDegree, u));
1845 const scalar vBasisValue(vBasis_.
basisValue(vCPI, vDegree, v));
1846 const scalar vBasisDeriv
1848 const scalar vBasis2Deriv
1851 NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1852 dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1853 d2Mdv2NWP += vBasis2Deriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1854 NMW += vBasisValue * uBasisValue * weights_[CPI];
1855 dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1856 d2Mdv2NW += vBasis2Deriv * uBasisValue * weights_[CPI];
1860 const vector vvDerivative
1864 - scalar(2)*dMdvNW*dMdvNWP/(NMW+SMALL)
1865 - d2Mdv2NW*NMWP/(NMW+SMALL)
1866 + scalar(2)*dMdvNW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1870 return vvDerivative;
1883 const label uDegree(uBasis_.
degree());
1884 const label vDegree(vBasis_.
degree());
1885 const label uNCPs(uBasis_.
nCPs());
1886 const label vNCPs(vBasis_.
nCPs());
1887 const label uCPI(CPsUCPIs_[CPI]);
1888 const label vCPI(CPsVCPIs_[CPI]);
1891 for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1893 for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1895 const label CPJ(vCPJ*uNCPs + uCPJ);
1896 const scalar uBasisValue(uBasis_.
basisValue(uCPJ, uDegree, u));
1897 const scalar vBasisValue(vBasis_.
basisValue(vCPJ, vDegree, v));
1899 NMW += vBasisValue * uBasisValue * weights_[CPJ];
1904 const scalar CPDerivative
1912 return CPDerivative;
1924 const label uDegree(uBasis_.
degree());
1925 const label vDegree(vBasis_.
degree());
1926 const label uNCPs(uBasis_.
nCPs());
1927 const label vNCPs(vBasis_.
nCPs());
1928 const label uCPI(CPsUCPIs_[CPI]);
1929 const label vCPI(CPsVCPIs_[CPI]);
1933 for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1935 for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1937 const label CPJ(vCPJ*uNCPs + uCPJ);
1938 const scalar uBasisValue(uBasis_.
basisValue(uCPJ, uDegree, u));
1939 const scalar vBasisValue(vBasis_.
basisValue(vCPJ, vDegree, v));
1941 NMWP += uBasisValue * vBasisValue * weights_[CPJ] * CPs_[CPJ];
1942 NMW += uBasisValue * vBasisValue * weights_[CPJ];
1951 * (CPs_[CPI] - (NMWP / (NMW+SMALL)))
1961 const scalar vConst,
1962 const scalar uStart,
1971 scalar ulDerivative(
Zero);
1975 scalar& uLocal(localU[uI]);
1976 uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1982 for (label uI=0; uI<(nPts-1); uI++)
1987 (dxdu[uI+1]&d2xdu2[uI+1])/(
mag(dxdu[uI+1])+SMALL)
1988 + (dxdu[uI ]&d2xdu2[uI ])/(
mag(dxdu[uI ])+SMALL)
1990 * (localU[uI+1]-localU[uI]);
1993 return ulDerivative;
1999 const scalar uConst,
2000 const scalar vStart,
2009 scalar vlDerivative(
Zero);
2013 scalar& vLocal(localV[vI]);
2014 vLocal = vStart + scalar(vI)/scalar(nPts-1)*(vEnd-vStart);
2020 for (label vI=0; vI<(nPts-1); vI++)
2025 (dxdv[vI+1]&d2xdv2[vI+1])/(
mag(dxdv[vI+1])+SMALL)
2026 + (dxdv[vI ]&d2xdv2[vI ])/(
mag(dxdv[vI ])+SMALL)
2028 * (localV[vI+1]-localV[vI]);
2031 return vlDerivative;
2039 if (!boundaryCPIDs_)
2041 const label uNCPs(uBasis_.
nCPs());
2042 const label vNCPs(vBasis_.
nCPs());
2043 const label nBoundCPs(2*uNCPs+2*vNCPs-4);
2044 boundaryCPIDs_.reset(
new labelList(nBoundCPs, -1));
2045 whichBoundaryCPID_.reset(
new labelList(uNCPs*vNCPs, -1));
2049 for (label vI=0; vI<vNCPs; vI+=vNCPs-1)
2051 for (label uI=0; uI<uNCPs; uI++)
2053 const label CPI(vI*uNCPs + uI);
2054 whichBoundaryCPID_()[CPI] = bID;
2055 boundaryCPIDs_()[bID++] = CPI;
2059 for (label uI=0; uI<uNCPs; uI+=uNCPs-1)
2062 for (label vI=1; vI<vNCPs-1; vI++)
2064 const label CPI(vI*uNCPs + uI);
2065 whichBoundaryCPID_()[CPI] = bID;
2066 boundaryCPIDs_()[bID++] = CPI;
2071 return boundaryCPIDs_();
2083 if (!whichBoundaryCPID_)
2088 return whichBoundaryCPID_()[globalCPI];
2120 << CPs_[CPI].x() <<
" " 2121 << CPs_[CPI].y() <<
" " 2161 OFstream surfaceFile(dirName/fileName);
2162 OFstream surfaceFileCPs(dirName/fileName+
"CPs");
2177 << CPs_[CPI].x() <<
" " 2178 << CPs_[CPI].y() <<
" " 2242 << CPs_[CPI].x() <<
" " 2243 << CPs_[CPI].y() <<
" " 2244 << CPs_[CPI].z() <<
")" 2281 const fileName dirName,
2282 const fileName fileName
2287 OFstream surfaceFile(dirName/fileName);
2288 OFstream surfaceFileCPs(dirName/fileName+
"CPs");
2305 << CPs_[CPI].x() <<
" " 2306 << CPs_[CPI].y() <<
" " 2307 << CPs_[CPI].z() <<
")" 2344 const fileName vtkDirName,
2345 const fileName vtkFileName
2350 if (vtkFileName.has_ext())
2353 <<
"Do not supply a file extension." 2361 OFstream surfaceFile(vtkFileName);
2363 faceList surfaceFaces((nUPts_ - 1)*(nUPts_ - 1), face(4));
2365 for (label fuI = 0; fuI < (nUPts_ - 1); fuI++)
2367 for (label fvI = 0; fvI < (nVPts_ - 1); fvI++)
2369 const label fI(fuI*(nUPts_ - 1) + fvI);
2370 face& surfaceFace(surfaceFaces[fI]);
2372 surfaceFace[0] = (fuI)*nVPts_ + (fvI);
2373 surfaceFace[1] = (fuI + 1)*nVPts_ + (fvI);
2374 surfaceFace[2] = (fuI + 1)*nVPts_ + (fvI + 1);
2375 surfaceFace[3] = (fuI)*nVPts_ + (fvI + 1);
2379 surfaceWriters::vtkWriter
writer;
2385 vtkDirName/vtkFileName,
2392 fileName vtkCPFileName(vtkFileName+
"CPs");
2394 const label uNCPs(uBasis_.
nCPs());
2395 const label vNCPs(vBasis_.
nCPs());
2396 faceList surfaceCPFaces((uNCPs-1)*(vNCPs-1), face(4));
2398 for (label fvCPI=0; fvCPI<(vNCPs-1); fvCPI++)
2400 for (label fuCPI=0; fuCPI<(uNCPs-1); fuCPI++)
2402 const label fCPI(fvCPI*(uNCPs-1) + fuCPI);
2403 face& surfaceCPFace(surfaceCPFaces[fCPI]);
2405 surfaceCPFace[0] = (fvCPI)*uNCPs + (fuCPI);
2406 surfaceCPFace[1] = (fvCPI + 1)*uNCPs + (fuCPI);
2407 surfaceCPFace[2] = (fvCPI + 1)*uNCPs + (fuCPI + 1);
2408 surfaceCPFace[3] = (fvCPI)*uNCPs + (fuCPI + 1);
2416 vtkDirName/vtkCPFileName,
scalar lengthU(const label vIConst, const label uIStart, const label uIEnd) const
List< scalar > scalarList
List of scalar.
void size(const label n)
Older name for setAddressableSize.
vector surfaceDerivativeUV(const scalar u, const scalar v) const
Surface second derivative wrt u and v at point u,v.
A class for handling file names.
vector surfaceDerivativeVV(const scalar u, const scalar v) const
Surface second derivative wrt v at point u,v.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void makeEquidistant(const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Output to file stream as an OSstream, normally using std::ofstream for the actual output...
constexpr char nl
The newline '\n' character (0x0a)
const labelList & getBoundaryCPIDs()
Get IDs of boundary control points.
Ostream & endl(Ostream &os)
Add newline and flush stream.
scalar lengthV(const label uIConst, const label vIStart, const label vIEnd) const
scalar surfaceDerivativeCP(const scalar u, const scalar v, const label CPI) const
Surface derivative wrt the weight of CPI at point u,v.
void setCPs(const List< vector > &CPs)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
#define forAll(list, i)
Loop across all elements in list.
T & operator[](const label i)
Return element of UList.
vector surfacePoint(const scalar &u, const scalar &v)
void write()
Write curve to file.
List< face > faceList
List of faces.
const label & nCPs() const
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const vector nrm(scalar u, scalar v)
NURBSbasis function. Used to construct NURBS curves, surfaces and volumes.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
NURBS3DSurface(const List< vector > &CPs, const label nPointsU, const label nPointsV, const NURBSbasis &uBasis, const NURBSbasis &vBasis, const word name="NURBS3DSurface")
Construct from number of control points and basis functions.
const label & degree() const
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
vector surfaceDerivativeW(const scalar u, const scalar v, const label CPI) const
Surface derivative wrt WI at point u,v.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
bool checkRangeUV(const scalar v, const scalar u, const label CPI, const label uDegree, const label vDegree) const
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.
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.
bool checkRangeU(const scalar u, const label CPI, const label uDegree) const
bool checkRangeV(const scalar v, const label CPI, const label vDegree) const
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.
vector point
Point is a vector.
scalarList findClosestSurfacePoint(const vector &targetPoint, const label maxIter=100, const scalar tolerance=1.e-6)
#define R(A, B, C, D, E, F, K, M)
#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.
const dimensionedScalar c
Speed of light in a vacuum.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
vector surfaceDerivativeU(const scalar u, const scalar v) const
Surface derivative wrt u at point u,v.
vector surfaceDerivativeUU(const scalar u, const scalar v) const
Surface second derivative wrt u at point u,v.
const label & whichBoundaryCPI(const label &globalCPI)
Get the boundary CP ID based on the global CP ID.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
void writeVTK(const fileName vtkDirName, const fileName vtkFileName)
void setWeights(const scalarList &weights)
vector surfaceDerivativeV(const scalar u, const scalar v) const
Surface derivative wrt v at point u,v.
List< label > labelList
A List of labels.
void setNrmOrientation(const vector &givenNrm, const scalar u, const scalar v)
const labelList & getBoundaryCPIs()
void flipNrmOrientation()
Flip the orientation of the nrm.
void setName(const word &name)
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.
static constexpr const zero Zero
Global zero (0)