56 <<
"Attempting to recompute points residing within control boxes"
73 Info<<
"Control Points bounds \n"
74 <<
"\tX1 : (" << lowerX <<
" " << upperX <<
")\n"
75 <<
"\tX2 : (" << lowerY <<
" " << upperY <<
")\n"
76 <<
"\tX3 : (" << lowerZ <<
" " << upperZ <<
")\n" <<
endl;
81 const vector& pointI = meshPoints[pI];
84 pointI.
x() >= lowerX && pointI.
x() <= upperX
85 && pointI.
y() >= lowerY && pointI.
y() <= upperY
86 && pointI.
z() >= lowerZ && pointI.
z() <= upperZ
90 reverseMap[pI] = count;
99 Info<<
"Initially found " << count <<
" points inside control boxes"
109 scalar timeBef = mesh_.time().elapsedCpuTime();
111 if (parametricCoordinatesPtr_)
114 <<
"Attempting to recompute parametric coordinates"
119 labelList& reverseMap = reverseMapPtr_();
121 parametricCoordinatesPtr_.reset
127 "parametricCoordinates" + name_,
128 mesh_.time().timeName(),
137 vectorField& paramCoors = parametricCoordinatesPtr_().primitiveFieldRef();
142 parametricCoordinatesPtr_().typeHeaderOk<pointVectorField>(
true)
150 Info<<
"Reading parametric coordinates from file" <<
endl;
151 IOobject& header = parametricCoordinatesPtr_().ref();
152 parametricCoordinatesPtr_() =
166 const label globalPointIndex = map[pI];
169 actualMap[curIndex] = map[pI];
170 reverseMap[globalPointIndex] = curIndex;
175 reverseMap[globalPointIndex] = -1;
183 Info<<
"Read non-zero parametric coordinates for " << curIndex
184 <<
" points" <<
endl;
193 scalar minX1 =
min(cps_.component(0));
194 scalar maxX1 =
max(cps_.component(0));
195 scalar minX2 =
min(cps_.component(1));
196 scalar maxX2 =
max(cps_.component(1));
197 scalar minX3 =
min(cps_.component(2));
198 scalar maxX3 =
max(cps_.component(2));
200 scalar oneOverDenomX(1./(maxX1 - minX1));
201 scalar oneOverDenomY(1./(maxX2 - minX2));
202 scalar oneOverDenomZ(1./(maxX3 - minX3));
206 const label globalPI = map[pI];
207 paramCoors[globalPI].x() = (
points[pI].x() - minX1)*oneOverDenomX;
208 paramCoors[globalPI].y() = (
points[pI].y() - minX2)*oneOverDenomY;
209 paramCoors[globalPI].z() = (
points[pI].z() - minX3)*oneOverDenomZ;
215 label nDropedPoints(0);
223 Info<<
"Mapping of mesh points to parametric space for box " << name_
226 label maxIterNeeded(0);
230 label nBoundIters(0);
231 vector res(GREAT, GREAT, GREAT);
234 const label globalPI = map[pI];
235 vector& uVec = paramCoors[globalPI];
236 vector& coorPointI = splinesBasedCoors[pI];
237 uVec += ((
inv(JacobianUVW(uVec))) & (
points[pI] - coorPointI));
245 if (nBoundIters > nMaxBound_)
247 dropOffPoints[pI] =
true;
268 <<
"Mapping to parametric space for point " << pI
269 <<
" failed." <<
endl
270 <<
"Residual after " << maxIter_ + 1 <<
" iterations : "
272 <<
"parametric coordinates " << paramCoors[map[pI]]
274 <<
"Local system coordinates " <<
points[pI] <<
endl
275 <<
"Threshold residual per direction : " << tolerance_
278 maxIterNeeded =
max(maxIterNeeded, iter);
282 label nParameterizedPoints = map.
size() - nDropedPoints;
288 map.
setSize(nParameterizedPoints);
293 if (!dropOffPoints[pI])
295 map[curIndex] = mapOld[pI];
296 reverseMap[mapOld[pI]] = curIndex;
302 reverseMap[mapOld[pI]] = -1;
308 Info<<
"Found " << nDropedPoints
309 <<
" to discard from morphing boxes" <<
endl;
310 Info<<
"Keeping " << nParameterizedPoints
311 <<
" parameterized points in boxes" <<
endl;
314 scalar maxDiff(-GREAT);
315 forAll(splinesBasedCoors, pI)
318 mag(splinesBasedCoors[pI] - localSystemCoordinates_[map[pI]]);
325 scalar timeAft = mesh_.time().elapsedCpuTime();
326 Info<<
"\tMapping completed in " << timeAft - timeBef <<
" seconds"
328 Info<<
"\tMax iterations per point needed to compute parametric "
330 << maxIterNeeded <<
endl;
331 Info<<
"\tMax difference between original mesh points and "
332 <<
"parameterized ones "
344 computeParametricCoordinates(
points);
355 bool boundPoint(
false);
357 if (vec.
x() < scalar(0))
362 if (vec.
y() < scalar(0))
367 if (vec.
z() < scalar(0))
396 mkDir(mesh_.time().globalPath()/
"optimisation"/cpsFolder_);
403 label nCPs = cps_.size();
404 activeControlPoints_ =
boolList(nCPs,
true);
405 activeDesignVariables_ =
boolList(3*nCPs,
true);
408 confineBoundaryControlPoints();
411 continuityRealatedConfinement();
414 confineControlPointsDirections();
418 forAll(activeControlPoints_, cpI)
422 !activeDesignVariables_[3*cpI]
423 && !activeDesignVariables_[3*cpI + 1]
424 && !activeDesignVariables_[3*cpI + 2]
427 activeControlPoints_[cpI] =
false;
435 const label nCPsU = basisU_.nCPs();
436 const label nCPsV = basisV_.nCPs();
437 const label nCPsW = basisW_.nCPs();
440 if (confineBoundaryControlPoints_)
443 for (label iCPw = 0; iCPw < nCPsW; iCPw += nCPsW - 1)
445 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
447 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
449 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
454 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
456 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
458 for (label iCPu = 0; iCPu < nCPsU; iCPu += nCPsU - 1)
460 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
465 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
467 for (label iCPv = 0; iCPv < nCPsV; iCPv += nCPsV - 1)
469 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
471 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
481 const label nCPsU = basisU_.nCPs();
482 const label nCPsV = basisV_.nCPs();
483 const label nCPsW = basisW_.nCPs();
487 forAll(confineUMinCPs_, iCPu)
489 const boolVector& confineSlice = confineUMinCPs_[iCPu];
491 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
493 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
495 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
500 forAll(confineUMaxCPs_, sliceI)
502 const boolVector& confineSlice = confineUMaxCPs_[sliceI];
503 label iCPu = nCPsU - 1 - sliceI;
505 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
507 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
509 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
516 forAll(confineVMinCPs_, iCPv)
518 const boolVector& confineSlice = confineVMinCPs_[iCPv];
520 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
522 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
524 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
529 forAll(confineVMaxCPs_, sliceI)
531 const boolVector& confineSlice = confineVMaxCPs_[sliceI];
532 label iCPv = nCPsV - 1 - sliceI;
534 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
536 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
538 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
545 forAll(confineWMinCPs_, iCPw)
547 const boolVector& confineSlice = confineWMinCPs_[iCPw];
549 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
551 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
553 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
558 forAll(confineWMaxCPs_, sliceI)
560 const boolVector& confineSlice = confineWMaxCPs_[sliceI];
561 label iCPw = nCPsW - 1 - sliceI;
563 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
565 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
567 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
576 for (label cpI = 0; cpI < cps_.size(); ++cpI)
578 if (confineUMovement_) activeDesignVariables_[3*cpI] =
false;
579 if (confineVMovement_) activeDesignVariables_[3*cpI + 1] =
false;
580 if (confineWMovement_) activeDesignVariables_[3*cpI + 2] =
false;
587 if (cpI < 0 || cpI > cps_.size() -1)
590 <<
"Attempted to confine control point movement for a control point "
591 <<
" ID which is out of bounds"
596 activeDesignVariables_[3*cpI] =
false;
597 activeDesignVariables_[3*cpI + 1] =
false;
598 activeDesignVariables_[3*cpI + 2] =
false;
609 if (cpI < 0 || cpI > cps_.size() -1)
612 <<
"Attempted to confine control point movement for a control point "
613 <<
" ID which is out of bounds"
618 if (confineDirections.
x()) activeDesignVariables_[3*cpI] =
false;
619 if (confineDirections.
y()) activeDesignVariables_[3*cpI + 1] =
false;
620 if (confineDirections.
z()) activeDesignVariables_[3*cpI + 2] =
false;
631 bool computeParamCoors
650 basisU_(
dict.get<label>(
"nCPsU"),
dict.get<label>(
"degreeU")),
651 basisV_(
dict.get<label>(
"nCPsV"),
dict.get<label>(
"degreeV")),
652 basisW_(
dict.get<label>(
"nCPsW"),
dict.get<label>(
"degreeW")),
653 maxIter_(
dict.getOrDefault<label>(
"maxIterations", 10)),
654 tolerance_(
dict.getOrDefault<scalar>(
"tolerance", 1.e-10)),
655 nMaxBound_(
dict.getOrDefault<scalar>(
"nMaxBoundIterations", 4)),
658 reverseMapPtr_(nullptr),
659 parametricCoordinatesPtr_(nullptr),
665 "confineUMovement", {{
"confineX1movement", 1912}}, false
672 "confineVMovement", {{
"confineX2movement", 1912}}, false
679 "confineWMovement", {{
"confineX3movement", 1912}}, false
682 confineBoundaryControlPoints_
690 "confineUMinCPs", {{
"boundUMinCPs", 1912}}, boolVectorList()
697 "confineUMaxCPs", {{
"boundUMaxCPs", 1912}}, boolVectorList()
704 "confineVMinCPs", {{
"boundVMinCPs", 1912}}, boolVectorList()
711 "confineVMaxCPs", {{
"boundVMaxCPs", 1912}}, boolVectorList()
718 "confineWMinCPs", {{
"boundWMinCPs", 1912}}, boolVectorList()
725 "confineWMaxCPs", {{
"boundWMaxCPs", 1912}}, boolVectorList()
728 activeControlPoints_(0),
729 activeDesignVariables_(0),
730 cpsFolder_(
"controlPoints"),
739 (confineUMinCPs_.size() + confineUMaxCPs_.size() >= basisU_.nCPs())
740 || (confineVMinCPs_.size() + confineVMaxCPs_.size() >= basisV_.nCPs())
741 || (confineWMinCPs_.size() + confineWMaxCPs_.size() >= basisW_.nCPs())
745 <<
"Number of control point slices to be kept frozen at "
746 <<
"the boundaries is invalid \n"
747 <<
"Number of control points in u " << basisU_.nCPs() <<
"\n"
748 <<
"Number of control points in v " << basisV_.nCPs() <<
"\n"
749 <<
"Number of control points in w " << basisW_.nCPs() <<
"\n"
754 if (
found(
"controlPoints"))
761 basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()
768 determineActiveDesignVariablesAndPoints();
778 bool computeParamCoors
783 Info<<
"NURBS3DVolume type : " << modelType <<
endl;
785 auto* ctorPtr = dictionaryConstructorTable(modelType);
794 *dictionaryConstructorTablePtr_
811 const label degreeU = basisU_.degree();
812 const label degreeV = basisV_.degree();
813 const label degreeW = basisW_.degree();
815 const label nCPsU = basisU_.nCPs();
816 const label nCPsV = basisV_.nCPs();
817 const label nCPsW = basisW_.nCPs();
821 for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
823 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
824 for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
826 const scalar basisVW = basisW*basisV_.basisValue(iCPv, degreeV, v);
827 for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
830 cps_[getCPID(iCPu, iCPv, iCPw)]
831 *basisU_.basisDerivativeU(iCPu, degreeU, u)
848 const label degreeU = basisU_.degree();
849 const label degreeV = basisV_.degree();
850 const label degreeW = basisW_.degree();
852 const label nCPsU = basisU_.nCPs();
853 const label nCPsV = basisV_.nCPs();
854 const label nCPsW = basisW_.nCPs();
858 for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
860 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
861 for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
863 const scalar basisWDeriV =
864 basisW*basisV_.basisDerivativeU(iCPv, degreeV, v);
865 for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
868 cps_[getCPID(iCPu, iCPv, iCPw)]
869 *basisU_.basisValue(iCPu, degreeU, u)
886 const label degreeU = basisU_.degree();
887 const label degreeV = basisV_.degree();
888 const label degreeW = basisW_.degree();
890 const label nCPsU = basisU_.nCPs();
891 const label nCPsV = basisV_.nCPs();
892 const label nCPsW = basisW_.nCPs();
896 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
898 const scalar derivW(basisW_.basisDerivativeU(iCPw, degreeW, w));
899 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
901 const scalar derivWBasisV =
902 derivW*basisV_.basisValue(iCPv, degreeV, v);
903 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
906 cps_[getCPID(iCPu, iCPv, iCPw)]
907 *basisU_.basisValue(iCPu, degreeU, u)
922 const scalar u = uVector.
x();
923 const scalar v = uVector.
y();
924 const scalar w = uVector.
z();
926 vector uDeriv = volumeDerivativeU(u, v, w);
927 vector vDeriv = volumeDerivativeV(u, v, w);
928 vector wDeriv = volumeDerivativeW(u, v, w);
930 tensor Jacobian(uDeriv, vDeriv, wDeriv,
true);
942 const scalar u = uVector.
x();
943 const scalar v = uVector.
y();
944 const scalar w = uVector.
z();
946 const label nCPsU = basisU_.nCPs();
947 const label nCPsV = basisV_.nCPs();
949 const label degreeU = basisU_.degree();
950 const label degreeV = basisV_.degree();
951 const label degreeW = basisW_.degree();
953 label iCPw = cpI/label(nCPsU*nCPsV);
954 label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
955 label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
961 basisU_.basisValue(iCPu, degreeU, u)
962 *basisV_.basisValue(iCPv, degreeV, v)
963 *basisW_.basisValue(iCPw, degreeW, w);
978 const vectorField& parametricCoordinates = getParametricCoordinates();
980 forAll(controlPointDerivs, cpI)
982 forAll(sensitivityPatchIDs, pI)
984 const label patchI = sensitivityPatchIDs[pI];
985 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
986 const labelList& meshPoints = patch.meshPoints();
991 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
995 if (whichPointInBox != -1)
997 controlPointDerivs[cpI] +=
1015 return controlPointDerivs;
1026 computeControlPointSensitivities
1044 const vectorField& parametricCoordinates = getParametricCoordinates();
1048 const labelList& reverseMap = reverseMapPtr_();
1050 forAll(controlPointDerivs, cpI)
1052 forAll(sensitivityPatchIDs, pI)
1054 const label patchI = sensitivityPatchIDs[pI];
1055 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1056 const label patchStart = patch.start();
1062 const face& fGlobal = mesh_.faces()[fI + patchStart];
1069 const label whichPointInBox = reverseMap[
globalIndex];
1072 if (whichPointInBox != -1)
1076 facePointDerivs[pI] =
1078 * volumeDerivativeCP
1093 controlPointDerivs[cpI] += patchSens[fI] & fCtrs_d;
1100 return controlPointDerivs;
1114 const vectorField& parametricCoordinates = getParametricCoordinates();
1118 const labelList& reverseMap = reverseMapPtr_();
1120 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1121 const label patchStart = patch.start();
1125 const face& fGlobal = mesh_.faces()[fI + patchStart];
1132 const label whichPointInBox = reverseMap[
globalIndex];
1135 if (whichPointInBox != -1)
1139 facePointDerivs[pI] =
1155 cpSens += faceSens[fI] & fCtrs_d;
1168 bool DimensionedNormalSens
1171 const fvPatch& patch = mesh_.boundary()[patchI];
1172 const polyPatch& ppatch = patch.patch();
1178 const label patchStart = ppatch.
start();
1179 const labelList& reverseMap = reverseMapPtr_();
1182 const vectorField& parametricCoordinates = getParametricCoordinates();
1187 const face& fGlobal = mesh_.faces()[fI + patchStart];
1194 const label whichPointInBox = reverseMap[
globalIndex];
1197 if (whichPointInBox != -1)
1201 facePointDerivs[pI] =
1219 if (DimensionedNormalSens)
1221 dndbSens[fI] = dNdbSens[1];
1225 dndbSens[fI] = dNdbSens[2];
1240 const vectorField& parametricCoordinates = getParametricCoordinates();
1243 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1244 const labelList& meshPoints = patch.meshPoints();
1248 auto& dxdb = tdxdb.ref();
1253 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
1256 if (whichPointInBox != -1)
1279 const vectorField& parametricCoordinates = getParametricCoordinates();
1282 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1283 const label patchStart = patch.start();
1287 auto& dxdb = tdxdb.ref();
1294 const face& fGlobal = mesh_.faces()[fI + patchStart];
1301 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
1304 if (whichPointInBox != -1)
1308 facePointDerivs[pI] =
1335 const label
nPoints = mapPtr_().size();
1337 auto&
points = tpoints.ref();
1341 const label globalPI = mapPtr_()[pI];
1354 const label degreeU = basisU_.degree();
1355 const label degreeV = basisV_.degree();
1356 const label degreeW = basisW_.degree();
1358 const label nCPsU = basisU_.nCPs();
1359 const label nCPsV = basisV_.nCPs();
1360 const label nCPsW = basisW_.nCPs();
1362 const scalar u = uVector.
x();
1363 const scalar v = uVector.
y();
1364 const scalar w = uVector.
z();
1367 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
1369 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
1370 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
1372 const scalar basisVW =
1373 basisW*basisV_.basisValue(iCPv, degreeV, v);
1374 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
1377 cps_[getCPID(iCPu, iCPv, iCPw)]
1378 *basisU_.basisValue(iCPu, degreeU, u)
1394 const vectorField& paramCoors = getParametricCoordinates();
1398 cps_ += controlPointsMovement;
1399 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1403 const vectorField& parameterizedPoints = tparameterizedPoints();
1410 forAll(parameterizedPoints, pI)
1412 newPoints[map[pI]] = transformPointToCartesian(parameterizedPoints[pI]);
1416 updateLocalCoordinateSystem(newPoints);
1418 <<
"Max mesh movement equal to "
1429 const bool updateCPs
1433 const vectorField& paramCoors = getParametricCoordinates();
1436 cps_ += controlPointsMovement;
1440 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1448 for (
const label patchI : patchesToBeMoved)
1450 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1451 const labelList& meshPoints = patch.meshPoints();
1455 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
1458 if (whichPointInBox != -1)
1461 transformPointToCartesian
1475 updateLocalCoordinateSystem(newPoints);
1480 cps_ -= controlPointsMovement;
1484 <<
"Max mesh movement equal to "
1498 const label nCPsU = basisU_.nCPs();
1499 const label nCPsV = basisV_.nCPs();
1501 return k*nCPsU*nCPsV + j*nCPsU + i;
1507 if (cps_.size() != newCps.
size())
1510 <<
"Attempting to replace control points with a set of "
1523 forAll(controlPointsMovement, cpI)
1525 if (!activeDesignVariables_[3*cpI])
1527 controlPointsMovement[cpI].x() =
Zero;
1529 if (!activeDesignVariables_[3*cpI + 1])
1531 controlPointsMovement[cpI].y() =
Zero;
1533 if (!activeDesignVariables_[3*cpI + 2])
1535 controlPointsMovement[cpI].z() =
Zero;
1550 const vectorField& paramCoors = getParametricCoordinates();
1552 cps_ += controlPointsMovement;
1554 scalar maxDisplacement(
Zero);
1555 for (
const label patchI : patchesToBeMoved)
1557 const polyPatch& patch = mesh_.boundaryMesh()[patchI];
1558 const labelList& meshPoints = patch.meshPoints();
1562 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
1565 if (whichPointInBox != -1)
1568 transformPointToCartesian
1587 return maxDisplacement;
1595 findPointsInBox(localSystemCoordinates_);
1599 new vectorField(localSystemCoordinates_, mapPtr_())
1610 findPointsInBox(localSystemCoordinates_);
1619 if (!reverseMapPtr_)
1621 findPointsInBox(localSystemCoordinates_);
1624 return reverseMapPtr_();
1631 if (!parametricCoordinatesPtr_)
1637 findPointsInBox(localSystemCoordinates_);
1639 computeParametricCoordinates(getPointsInBox()());
1642 return parametricCoordinatesPtr_();
1649 const vectorField& parametricCoordinates = getParametricCoordinates();
1659 mesh_.time().timeName(),
1695 const vectorField& parametricCoordinates = getParametricCoordinates();
1705 mesh_.time().timeName(),
1736 const label cellI = pointCellsI[cI];
1737 DxDb[cellI] += C_d[cI] & pointDxDb;
1742 forAll(mesh_.boundary(), pI)
1744 const fvPatch& patch = mesh_.boundary()[pI];
1745 if (!isA<coupledFvPatch>(patch))
1760 label nU(basisU_.nCPs());
1761 return label(nU % 2 == 0 ? 0.5*nU : 0.5*(nU - 1) + 1);
1767 label nV(basisV_.nCPs());
1768 return label(nV % 2 == 0 ? 0.5*nV : 0.5*(nV - 1) + 1);
1774 label nW(basisW_.nCPs());
1775 return label(nW % 2 == 0 ? 0.5*nW : 0.5*(nW - 1) + 1);
1781 return Vector<label>(nUSymmetry(), nVSymmetry(), nWSymmetry());
1791 const label nCPsU = basisU_.nCPs();
1792 const label nCPsV = basisV_.nCPs();
1797 forAll(cpsInCartesian, cpI)
1799 cpsInCartesian[cpI] = transformPointToCartesian(cps_[cpI]);
1803 Info<<
"Writing control point positions to file" <<
endl;
1807 OFstream cpsFile(
"optimisation"/cpsFolder_/name_ + baseName +
".csv");
1810 <<
"\"Points : 0\", \"Points : 1\", \"Points : 2\","
1811 <<
"\"i\", \"j\", \"k\","
1812 <<
"\"active : 0\", \"active : 1\", \"active : 2\"" <<
endl;
1814 forAll(cpsInCartesian, cpI)
1816 const label iCPw = cpI/label(nCPsU*nCPsV);
1817 const label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
1818 const label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
1821 << cpsInCartesian[cpI].x() <<
", "
1822 << cpsInCartesian[cpI].y() <<
", "
1823 << cpsInCartesian[cpI].z() <<
", "
1827 << activeDesignVariables_[3*cpI] <<
", "
1828 << activeDesignVariables_[3*cpI + 1] <<
", "
1829 << activeDesignVariables_[3*cpI + 2] <<
endl;
1837 parametricCoordinatesPtr_().write();
1843 cps_.writeEntry(
"controlPoints",
os);
Macros for easy insertion into run-time selection tables.
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Generic GeometricBoundaryField class.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void setSize(const label n)
Alias for resize()
NURBS3DVolume morpher. Includes support functions for gradient computations Base class providing supp...
tmp< tensorField > patchDxDbFace(const label patchI, const label cpI)
Get patch dx/db.
void computeParametricCoordinates(const vectorField &points)
label nUSymmetry() const
Get number of variables if CPs are moved symmetrically in U.
autoPtr< labelList > mapPtr_
Map of points-in-box to mesh points.
vector volumeDerivativeU(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt u at point u,v,w.
void confineBoundaryControlPoints()
Confine movement in boundary control points if necessary.
const labelList & getMap()
Get map of points in box to mesh points.
virtual bool writeData(Ostream &os) const
Write the control points to support restart.
scalar computeMaxBoundaryDisplacement(const vectorField &controlPointsMovement, const labelList &patchesToBeMoved)
Compute max. displacement at the boundary.
void findPointsInBox(const vectorField &meshPoints)
Find points within control points box.
tmp< vectorField > computeNewBoundaryPoints(const vectorField &controlPointsMovement, const labelList &patchesToBeMoved, const bool moveCPs=true)
Boundary mesh movement based on given control point movement.
vectorField computeControlPointSensitivities(const pointVectorField &pointSens, const labelList &sensitivityPatchIDs)
const pointVectorField & getParametricCoordinates()
Get parametric coordinates.
tmp< volTensorField > getDxCellsDb(const label cpI)
Get dxCartesiandb for a certain control point on cells.
autoPtr< labelList > reverseMapPtr_
Map of mesh points to points-in-box.
label nVSymmetry() const
Get number of variables if CPs are moved symmetrically in V.
tmp< pointTensorField > getDxDb(const label cpI)
Get dxCartesiandb for a certain control point.
void makeFolders()
Create folders to store cps and derivatives.
void writeCps(const fileName &="cpsFile", const bool transform=true) const
tmp< tensorField > patchDxDb(const label patchI, const label cpI)
Get patch dx/db.
tmp< vectorField > getPointsInBox()
Get mesh points that reside within the control points box.
vector volumeDerivativeV(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt v at point u,v,w.
tmp< tensorField > dndbBasedSensitivities(const label patchI, const label cpI, bool DimensionedNormalSens=true)
void confineControlPoint(const label cpI)
Confine all three movements for a prescribed control point.
Vector< label > nSymmetry() const
vectorField cps_
The volumetric B-Splines control points.
void boundControlPointMovement(vectorField &controlPointsMovement) const
void continuityRealatedConfinement()
Confine control point movement to maintain user-defined continuity.
label getCPID(const label i, const label j, const label k) const
Get control point ID from its I-J-K coordinates.
scalar volumeDerivativeCP(const vector &u, const label cpI) const
Volume point derivative wrt to control point cpI at point u,v,w.
label nWSymmetry() const
Get number of variables if CPs are moved symmetrically in W.
void writeParamCoordinates() const
Write parametric coordinates.
void determineActiveDesignVariablesAndPoints()
Create lists with active design variables and control points.
tmp< vectorField > computeNewPoints(const vectorField &controlPointsMovement)
Mesh movement based on given control point movement.
tensor JacobianUVW(const vector &u) const
Jacobian matrix wrt to the volume parametric coordinates.
void setControlPoints(const vectorField &newCps)
Set new control points.
const labelList & getReverseMap()
vector volumeDerivativeW(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt w at point u,v,w.
void confineControlPointsDirections()
Confine movement in all control points for user-defined directions.
Output to file stream, using an OSstream.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
const Cmpt & component(const direction) const
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.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Specialized bundling of boolean values as a vector of 3 components, element access using x(),...
bool x() const
The x component.
bool y() const
The y component.
bool z() const
The z component.
Differentiation of the mesh data structure.
tmp< tensorField > cellCenters_d(const label pointI)
vectorField makeFaceCentresAndAreas_d(const pointField &p, const pointField &p_d)
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
T getOrDefaultCompat(const word &keyword, std::initializer_list< std::pair< const char *, int > > compat, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
static const char *const coordinates
The keyword "coordinates".
A face is a list of labels corresponding to mesh vertices.
pointField points(const UList< point > &pts) const
Return the points corresponding to this face.
A class for handling file names.
Mesh data needed to do the Finite Volume discretisation.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
localIOdictionary is derived from IOdictionary but excludes parallel master reading.
Smooth ATC in cells having a point to a set of patches supplied by type.
A patch is a list of labels that address the faces in the global face list.
label start() const
Return start label of this patch in the polyMesh face list.
splitCell * master() const
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
PtrList< coordinateSystem > coordinates(solidRegions.size())
#define FatalIOErrorInLookup(ios, lookupTag, lookupName, lookupTable)
Report an error message using Foam::FatalIOError.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
const word dictName("faMeshDefinition")
#define DebugInfo
Report an information message using Foam::Info.
#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.
const dimensionSet dimless
Dimensionless.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
List< label > labelList
A List of labels.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
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)
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Type gMax(const FieldField< Field, Type > &f)
#define defineRunTimeSelectionTable(baseType, argNames)
Define run-time selection table.
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.