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];
167 if (paramCoors[globalPointIndex] != vector::zero)
169 actualMap[curIndex] = map[pI];
170 reverseMap[globalPointIndex] = curIndex;
175 reverseMap[globalPointIndex] = -1;
180 actualMap.setSize(curIndex);
182 reduce(curIndex, sumOp<label>());
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;
214 boolList dropOffPoints(map.size(),
false);
215 label nDropedPoints(0);
218 tmp<vectorField> tsplinesBasedCoors(
coordinates(paramCoors));
219 vectorField& splinesBasedCoors = tsplinesBasedCoors.ref();
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;
260 res.component(0) > tolerance_
261 || res.component(1) > tolerance_
262 || res.component(2) > tolerance_
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);
280 reduce(maxIterNeeded, maxOp<label>());
282 label nParameterizedPoints = map.size() - nDropedPoints;
288 map.setSize(nParameterizedPoints);
293 if (!dropOffPoints[pI])
295 map[curIndex] = mapOld[pI];
296 reverseMap[mapOld[pI]] = curIndex;
301 paramCoors[mapOld[pI]] = vector::zero;
302 reverseMap[mapOld[pI]] = -1;
306 reduce(nDropedPoints, sumOp<label>());
307 reduce(nParameterizedPoints, sumOp<label>());
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]]);
324 reduce(maxDiff, maxOp<scalar>());
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 " 340 tmp<vectorField> tPoints
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;
432 confineInertControlPoints();
447 const label pU = basisU_.degree();
448 const label pV = basisV_.degree();
449 const label pW = basisW_.degree();
450 forAll(activeControlPoints_, cpI)
452 if (activeControlPoints_[cpI])
455 label i(-1), j(-1),
k(-1);
456 getIJK(i, j,
k, cpI);
458 const scalar uMin(knotsU[i]);
459 const scalar uMax(knotsU[i + pU + 1]);
460 const scalar vMin(knotsV[j]);
461 const scalar vMax(knotsV[j + pV + 1]);
462 const scalar wMin(knotsW[
k]);
463 const scalar wMax(knotsW[
k + pW + 1]);
464 bool foundParamPt(
false);
465 for (
const label paramI : map)
467 const vector& paramCoors = parametricCoors()[paramI];
470 paramCoors.x() >= uMin && paramCoors.x() < uMax
471 && paramCoors.y() >= vMin && paramCoors.y() < vMax
472 && paramCoors.z() >= wMin && paramCoors.z() < wMax
479 reduce(foundParamPt, orOp<bool>());
482 activeControlPoints_[cpI] =
false;
483 activeDesignVariables_[3*cpI] =
false;
484 activeDesignVariables_[3*cpI + 1] =
false;
485 activeDesignVariables_[3*cpI + 2] =
false;
486 Info<<
"Disabling control " << cpI
488 <<
"(" << 3*cpI <<
", " << 3*cpI+1 <<
", " << 3*cpI+2 <<
")" 489 <<
" since they does not parameterize any mesh point" 499 const label nCPsU = basisU_.nCPs();
500 const label nCPsV = basisV_.nCPs();
501 const label nCPsW = basisW_.nCPs();
504 if (confineBoundaryControlPoints_)
507 for (label iCPw = 0; iCPw < nCPsW; iCPw += nCPsW - 1)
509 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
511 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
513 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
518 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
520 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
522 for (label iCPu = 0; iCPu < nCPsU; iCPu += nCPsU - 1)
524 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
529 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
531 for (label iCPv = 0; iCPv < nCPsV; iCPv += nCPsV - 1)
533 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
535 confineControlPoint(getCPID(iCPu, iCPv, iCPw));
545 const label nCPsU = basisU_.nCPs();
546 const label nCPsV = basisV_.nCPs();
547 const label nCPsW = basisW_.nCPs();
551 forAll(confineUMinCPs_, iCPu)
553 const boolVector& confineSlice = confineUMinCPs_[iCPu];
555 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
557 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
559 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
564 forAll(confineUMaxCPs_, sliceI)
566 const boolVector& confineSlice = confineUMaxCPs_[sliceI];
567 label iCPu = nCPsU - 1 - sliceI;
569 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
571 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
573 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
580 forAll(confineVMinCPs_, iCPv)
582 const boolVector& confineSlice = confineVMinCPs_[iCPv];
584 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
586 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
588 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
593 forAll(confineVMaxCPs_, sliceI)
595 const boolVector& confineSlice = confineVMaxCPs_[sliceI];
596 label iCPv = nCPsV - 1 - sliceI;
598 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
600 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
602 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
609 forAll(confineWMinCPs_, iCPw)
611 const boolVector& confineSlice = confineWMinCPs_[iCPw];
613 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
615 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
617 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
622 forAll(confineWMaxCPs_, sliceI)
624 const boolVector& confineSlice = confineWMaxCPs_[sliceI];
625 label iCPw = nCPsW - 1 - sliceI;
627 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
629 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
631 confineControlPoint(getCPID(iCPu, iCPv, iCPw), confineSlice);
640 for (label cpI = 0; cpI < cps_.size(); ++cpI)
642 if (confineUMovement_) activeDesignVariables_[3*cpI] =
false;
643 if (confineVMovement_) activeDesignVariables_[3*cpI + 1] =
false;
644 if (confineWMovement_) activeDesignVariables_[3*cpI + 2] =
false;
651 if (cpI < 0 || cpI > cps_.size() -1)
654 <<
"Attempted to confine control point movement for a control point " 655 <<
" ID which is out of bounds" 660 activeDesignVariables_[3*cpI] =
false;
661 activeDesignVariables_[3*cpI + 1] =
false;
662 activeDesignVariables_[3*cpI + 2] =
false;
670 const boolVector& confineDirections
673 if (cpI < 0 || cpI > cps_.size() -1)
676 <<
"Attempted to confine control point movement for a control point " 677 <<
" ID which is out of bounds" 682 if (confineDirections.x()) activeDesignVariables_[3*cpI] =
false;
683 if (confineDirections.y()) activeDesignVariables_[3*cpI + 1] =
false;
684 if (confineDirections.z()) activeDesignVariables_[3*cpI + 2] =
false;
693 const dictionary&
dict,
695 bool computeParamCoors
704 fileName(
"uniform")/fileName(
"volumetricBSplines"),
706 IOobject::READ_IF_PRESENT,
717 maxIter_(
dict.getOrDefault<label>(
"maxIterations", 10)),
718 tolerance_(
dict.getOrDefault<scalar>(
"tolerance", 1.
e-10)),
719 nMaxBound_(
dict.getOrDefault<scalar>(
"nMaxBoundIterations", 4)),
722 reverseMapPtr_(nullptr),
723 parametricCoordinatesPtr_(nullptr),
727 dict.getOrDefaultCompat<bool>
729 "confineUMovement", {{
"confineX1movement", 1912}}, false
736 "confineVMovement", {{
"confineX2movement", 1912}}, false
743 "confineWMovement", {{
"confineX3movement", 1912}}, false
746 confineBoundaryControlPoints_
754 "confineUMinCPs", {{
"boundUMinCPs", 1912}}, boolVectorList()
761 "confineUMaxCPs", {{
"boundUMaxCPs", 1912}}, boolVectorList()
768 "confineVMinCPs", {{
"boundVMinCPs", 1912}}, boolVectorList()
775 "confineVMaxCPs", {{
"boundVMaxCPs", 1912}}, boolVectorList()
782 "confineWMinCPs", {{
"boundWMinCPs", 1912}}, boolVectorList()
789 "confineWMaxCPs", {{
"boundWMaxCPs", 1912}}, boolVectorList()
792 activeControlPoints_(0),
793 activeDesignVariables_(0),
794 cpsFolder_(
"controlPoints"),
803 (confineUMinCPs_.size() + confineUMaxCPs_.size() >= basisU_.nCPs())
804 || (confineVMinCPs_.size() + confineVMaxCPs_.size() >= basisV_.nCPs())
805 || (confineWMinCPs_.size() + confineWMaxCPs_.size() >= basisW_.nCPs())
809 <<
"Number of control point slices to be kept frozen at " 810 <<
"the boundaries is invalid \n" 811 <<
"Number of control points in u " << basisU_.nCPs() <<
"\n" 812 <<
"Number of control points in v " << basisV_.nCPs() <<
"\n" 813 <<
"Number of control points in w " << basisW_.nCPs() <<
"\n" 818 if (
found(
"controlPoints"))
825 basisU_.nCPs()*basisV_.nCPs()*basisW_.nCPs()
839 const dictionary&
dict,
841 bool computeParamCoors
844 const word modelType(
dict.
get<word>(
"type"));
846 Info<<
"NURBS3DVolume type : " << modelType <<
endl;
848 auto* ctorPtr = dictionaryConstructorTable(modelType);
857 *dictionaryConstructorTablePtr_
861 return autoPtr<NURBS3DVolume>(ctorPtr(
dict,
mesh, computeParamCoors));
876 scalar minX1 =
min(cps_.component(0));
877 scalar maxX1 =
max(cps_.component(0));
878 scalar minX2 =
min(cps_.component(1));
879 scalar maxX2 =
max(cps_.component(1));
880 scalar minX3 =
min(cps_.component(2));
881 scalar maxX3 =
max(cps_.component(2));
883 scalar oneOverDenomX(1./(maxX1 - minX1));
884 scalar oneOverDenomY(1./(maxX2 - minX2));
885 scalar oneOverDenomZ(1./(maxX3 - minX3));
889 paramCoors[pI].x() = (
points[pI].x() - minX1)*oneOverDenomX;
890 paramCoors[pI].y() = (
points[pI].y() - minX2)*oneOverDenomY;
891 paramCoors[pI].z() = (
points[pI].z() - minX3)*oneOverDenomZ;
897 label nDropedPoints(0);
900 tmp<vectorField> tsplinesBasedCoors(
coordinates(paramCoors));
901 vectorField& splinesBasedCoors = tsplinesBasedCoors.ref();
905 Info<<
"Mapping of mesh points to parametric space for box " << name_
908 label maxIterNeeded(0);
912 label nBoundIters(0);
913 vector res(GREAT, GREAT, GREAT);
916 vector& uVec = paramCoors[pI];
917 vector& coorPointI = splinesBasedCoors[pI];
918 uVec += ((
inv(JacobianUVW(uVec))) & (
points[pI] - coorPointI));
926 if (nBoundIters > nMaxBound_)
928 dropOffPoints[pI] =
true;
941 res.component(0) > tolerance_
942 || res.component(1) > tolerance_
943 || res.component(2) > tolerance_
949 <<
"Mapping to parametric space for point " << pI
950 <<
" failed." <<
endl 951 <<
"Residual after " << maxIter_ + 1 <<
" iterations : " 953 <<
"parametric coordinates " << paramCoors[pI]
955 <<
"Local system coordinates " <<
points[pI] <<
endl 956 <<
"Threshold residual per direction : " << tolerance_
959 maxIterNeeded =
max(maxIterNeeded, iter);
961 reduce(maxIterNeeded, maxOp<label>());
963 label nParameterizedPoints =
points.
size() - nDropedPoints;
965 reduce(nDropedPoints, sumOp<label>());
966 reduce(nParameterizedPoints, sumOp<label>());
967 Info<<
"Found " << nDropedPoints
968 <<
" to discard from morphing boxes" <<
endl;
969 Info<<
"Keeping " << nParameterizedPoints
970 <<
" parameterized points in boxes" <<
endl;
982 const label degreeU = basisU_.degree();
983 const label degreeV = basisV_.degree();
984 const label degreeW = basisW_.degree();
986 const label nCPsU = basisU_.nCPs();
987 const label nCPsV = basisV_.nCPs();
988 const label nCPsW = basisW_.nCPs();
992 for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
994 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
995 for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
997 const scalar basisVW = basisW*basisV_.basisValue(iCPv, degreeV, v);
998 for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
1001 cps_[getCPID(iCPu, iCPv, iCPw)]
1002 *basisU_.basisDerivativeU(iCPu, degreeU, u)
1019 const label degreeU = basisU_.degree();
1020 const label degreeV = basisV_.degree();
1021 const label degreeW = basisW_.degree();
1023 const label nCPsU = basisU_.nCPs();
1024 const label nCPsV = basisV_.nCPs();
1025 const label nCPsW = basisW_.nCPs();
1029 for (label iCPw = 0; iCPw < nCPsW; ++iCPw)
1031 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
1032 for (label iCPv = 0; iCPv < nCPsV; ++iCPv)
1034 const scalar basisWDeriV =
1035 basisW*basisV_.basisDerivativeU(iCPv, degreeV, v);
1036 for (label iCPu = 0; iCPu < nCPsU; ++iCPu)
1039 cps_[getCPID(iCPu, iCPv, iCPw)]
1040 *basisU_.basisValue(iCPu, degreeU, u)
1057 const label degreeU = basisU_.degree();
1058 const label degreeV = basisV_.degree();
1059 const label degreeW = basisW_.degree();
1061 const label nCPsU = basisU_.nCPs();
1062 const label nCPsV = basisV_.nCPs();
1063 const label nCPsW = basisW_.nCPs();
1067 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
1069 const scalar derivW(basisW_.basisDerivativeU(iCPw, degreeW, w));
1070 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
1072 const scalar derivWBasisV =
1073 derivW*basisV_.basisValue(iCPv, degreeV, v);
1074 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
1077 cps_[getCPID(iCPu, iCPv, iCPw)]
1078 *basisU_.basisValue(iCPu, degreeU, u)
1093 const scalar u = uVector.x();
1094 const scalar v = uVector.y();
1095 const scalar w = uVector.z();
1097 vector uDeriv = volumeDerivativeU(u, v, w);
1098 vector vDeriv = volumeDerivativeV(u, v, w);
1099 vector wDeriv = volumeDerivativeW(u, v, w);
1101 tensor Jacobian(uDeriv, vDeriv, wDeriv,
true);
1113 const scalar u = uVector.x();
1114 const scalar v = uVector.y();
1115 const scalar w = uVector.z();
1117 const label nCPsU = basisU_.nCPs();
1118 const label nCPsV = basisV_.nCPs();
1120 const label degreeU = basisU_.degree();
1121 const label degreeV = basisV_.degree();
1122 const label degreeW = basisW_.degree();
1124 label iCPw = cpI/label(nCPsU*nCPsV);
1125 label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
1126 label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
1132 basisU_.basisValue(iCPu, degreeU, u)
1133 *basisV_.basisValue(iCPv, degreeV, v)
1134 *basisW_.basisValue(iCPw, degreeW, w);
1149 const vectorField& parametricCoordinates = getParametricCoordinates();
1151 forAll(controlPointDerivs, cpI)
1153 forAll(sensitivityPatchIDs, pI)
1155 const label patchI = sensitivityPatchIDs[pI];
1162 const label whichPointInBox = reverseMapPtr_()[
globalIndex];
1166 if (whichPointInBox != -1)
1168 controlPointDerivs[cpI] +=
1186 return controlPointDerivs;
1197 computeControlPointSensitivities
1215 const vectorField& parametricCoordinates = getParametricCoordinates();
1219 const labelList& reverseMap = reverseMapPtr_();
1221 forAll(controlPointDerivs, cpI)
1223 forAll(sensitivityPatchIDs, pI)
1225 const label patchI = sensitivityPatchIDs[pI];
1227 const label patchStart =
patch.start();
1233 const face& fGlobal = mesh_.faces()[fI + patchStart];
1240 const label whichPointInBox = reverseMap[
globalIndex];
1243 if (whichPointInBox != -1)
1247 facePointDerivs[pI] =
1249 * volumeDerivativeCP
1264 controlPointDerivs[cpI] += patchSens[fI] & fCtrs_d;
1271 return controlPointDerivs;
1285 const vectorField& parametricCoordinates = getParametricCoordinates();
1289 const labelList& reverseMap = reverseMapPtr_();
1292 const label patchStart =
patch.start();
1296 const face& fGlobal = mesh_.faces()[fI + patchStart];
1303 const label whichPointInBox = reverseMap[
globalIndex];
1306 if (whichPointInBox != -1)
1310 facePointDerivs[pI] =
1326 cpSens += faceSens[fI] & fCtrs_d;
1339 bool DimensionedNormalSens
1349 const label patchStart = ppatch.
start();
1350 const labelList& reverseMap = reverseMapPtr_();
1353 const vectorField& parametricCoordinates = getParametricCoordinates();
1358 const face& fGlobal = mesh_.faces()[fI + patchStart];
1365 const label whichPointInBox = reverseMap[
globalIndex];
1368 if (whichPointInBox != -1)
1372 facePointDerivs[pI] =
1390 if (DimensionedNormalSens)
1392 dndbSens[fI] = dNdbSens[1];
1396 dndbSens[fI] = dNdbSens[2];
1411 const vectorField& parametricCoordinates = getParametricCoordinates();
1414 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1419 auto& dxdb = tdxdb.ref();
1423 const label globalIndex = meshPoints[pI];
1424 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1427 if (whichPointInBox != -1)
1430 transformationTensorDxDb(globalIndex)
1433 parametricCoordinates[globalIndex],
1450 const vectorField& parametricCoordinates = getParametricCoordinates();
1453 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1454 const label patchStart =
patch.start();
1458 auto& dxdb = tdxdb.ref();
1461 deltaBoundary deltaBound(mesh_);
1465 const face& fGlobal = mesh_.faces()[fI + patchStart];
1466 const pointField facePoints = fGlobal.points(mesh_.points());
1471 const label globalIndex = fGlobal[pI];
1472 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1475 if (whichPointInBox != -1)
1479 facePointDerivs[pI] =
1480 transformationTensorDxDb(globalIndex)
1483 parametricCoordinates[globalIndex],
1490 deltaBound.makeFaceCentresAndAreas_d
1506 const label
nPoints = mapPtr_().size();
1508 auto&
points = tpoints.ref();
1512 const label globalPI = mapPtr_()[pI];
1525 const label degreeU = basisU_.degree();
1526 const label degreeV = basisV_.degree();
1527 const label degreeW = basisW_.degree();
1529 const label nCPsU = basisU_.nCPs();
1530 const label nCPsV = basisV_.nCPs();
1531 const label nCPsW = basisW_.nCPs();
1533 const scalar u = uVector.x();
1534 const scalar v = uVector.y();
1535 const scalar w = uVector.z();
1538 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
1540 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
1541 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
1543 const scalar basisVW =
1544 basisW*basisV_.basisValue(iCPv, degreeV, v);
1545 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
1548 cps_[getCPID(iCPu, iCPv, iCPw)]
1549 *basisU_.basisValue(iCPu, degreeU, u)
1565 const vectorField& paramCoors = getParametricCoordinates();
1569 cps_ += controlPointsMovement;
1570 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1573 tmp<vectorField> tparameterizedPoints =
coordinates(paramCoors);
1574 const vectorField& parameterizedPoints = tparameterizedPoints();
1577 tmp<vectorField> tnewPoints(
new vectorField(mesh_.points()));
1581 forAll(parameterizedPoints, pI)
1583 newPoints[map[pI]] = transformPointToCartesian(parameterizedPoints[pI]);
1587 updateLocalCoordinateSystem(newPoints);
1589 <<
"Max mesh movement equal to " 1600 const bool updateCPs
1604 const vectorField& paramCoors = getParametricCoordinates();
1607 cps_ += controlPointsMovement;
1611 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1615 tmp<vectorField> tnewPoints(
new vectorField(mesh_.points()));
1619 for (
const label patchI : patchesToBeMoved)
1621 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1624 for (
const label globalIndex : meshPoints)
1626 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1629 if (whichPointInBox != -1)
1631 newPoints[globalIndex] =
1632 transformPointToCartesian
1636 paramCoors[globalIndex]
1646 updateLocalCoordinateSystem(newPoints);
1651 cps_ -= controlPointsMovement;
1655 <<
"Max mesh movement equal to " 1669 const label nCPsU = basisU_.nCPs();
1670 const label nCPsV = basisV_.nCPs();
1672 return k*nCPsU*nCPsV + j*nCPsU + i;
1684 const label nCPsU = basisU_.nCPs();
1685 const label nCPsV = basisV_.nCPs();
1686 k = cpID/(nCPsU*nCPsV);
1687 const label inKplaneID = (cpID -
k*(nCPsU*nCPsV));
1688 j = inKplaneID/nCPsU;
1689 i = inKplaneID - j*nCPsU;
1695 if (cps_.size() != newCps.
size())
1698 <<
"Attempting to replace control points with a set of " 1711 forAll(controlPointsMovement, cpI)
1713 if (!activeDesignVariables_[3*cpI])
1715 controlPointsMovement[cpI].x() =
Zero;
1717 if (!activeDesignVariables_[3*cpI + 1])
1719 controlPointsMovement[cpI].y() =
Zero;
1721 if (!activeDesignVariables_[3*cpI + 2])
1723 controlPointsMovement[cpI].z() =
Zero;
1738 const vectorField& paramCoors = getParametricCoordinates();
1740 cps_ += controlPointsMovement;
1742 scalar maxDisplacement(
Zero);
1743 for (
const label patchI : patchesToBeMoved)
1745 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1748 for (
const label globalIndex : meshPoints)
1750 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1753 if (whichPointInBox != -1)
1756 transformPointToCartesian
1760 paramCoors[globalIndex]
1767 mag(newPoint - mesh_.points()[globalIndex])
1775 return maxDisplacement;
1783 findPointsInBox(localSystemCoordinates_);
1785 tmp<vectorField> pointsInBox
1798 findPointsInBox(localSystemCoordinates_);
1807 if (!reverseMapPtr_)
1809 findPointsInBox(localSystemCoordinates_);
1812 return reverseMapPtr_();
1819 if (!parametricCoordinatesPtr_)
1825 findPointsInBox(localSystemCoordinates_);
1827 computeParametricCoordinates(getPointsInBox()());
1830 return parametricCoordinatesPtr_();
1837 const vectorField& parametricCoordinates = getParametricCoordinates();
1840 tmp<pointTensorField> tDxDb
1847 mesh_.time().timeName(),
1862 for (
const label globalIndex : map)
1865 transformationTensorDxDb(globalIndex)
1868 parametricCoordinates[globalIndex],
1883 const vectorField& parametricCoordinates = getParametricCoordinates();
1886 tmp<volTensorField> tDxDb
1893 mesh_.time().timeName(),
1904 deltaBoundary deltaBound(mesh_);
1910 for (
const label globalIndex : map)
1913 transformationTensorDxDb(globalIndex)
1916 parametricCoordinates[globalIndex],
1919 const labelList& pointCellsI = pointCells[globalIndex];
1920 tmp<tensorField> tC_d = deltaBound.cellCenters_d(globalIndex);
1924 const label cellI = pointCellsI[cI];
1925 DxDb[cellI] += C_d[cI] & pointDxDb;
1930 forAll(mesh_.boundary(), pI)
1932 const fvPatch&
patch = mesh_.boundary()[pI];
1933 if (!isA<coupledFvPatch>(
patch))
1935 DxDb.boundaryFieldRef()[pI] = patchDxDbFace(pI, cpI);
1940 DxDb.correctBoundaryConditions();
1948 label nU(basisU_.nCPs());
1949 return label(nU % 2 == 0 ? 0.5*nU : 0.5*(nU - 1) + 1);
1955 label nV(basisV_.nCPs());
1956 return label(nV % 2 == 0 ? 0.5*nV : 0.5*(nV - 1) + 1);
1962 label nW(basisW_.nCPs());
1963 return label(nW % 2 == 0 ? 0.5*nW : 0.5*(nW - 1) + 1);
1969 return Vector<label>(nUSymmetry(), nVSymmetry(), nWSymmetry());
1979 const label nCPsU = basisU_.nCPs();
1980 const label nCPsV = basisV_.nCPs();
1985 forAll(cpsInCartesian, cpI)
1987 cpsInCartesian[cpI] = transformPointToCartesian(cps_[cpI]);
1991 Info<<
"Writing control point positions to file" <<
endl;
1995 OFstream cpsFile(
"optimisation"/cpsFolder_/name_ + baseName +
".csv");
1998 <<
"\"Points : 0\", \"Points : 1\", \"Points : 2\"," 1999 <<
"\"i\", \"j\", \"k\"," 2000 <<
"\"active : 0\", \"active : 1\", \"active : 2\"" <<
endl;
2002 forAll(cpsInCartesian, cpI)
2004 const label iCPw = cpI/label(nCPsU*nCPsV);
2005 const label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
2006 const label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
2009 << cpsInCartesian[cpI].x() <<
", " 2010 << cpsInCartesian[cpI].y() <<
", " 2011 << cpsInCartesian[cpI].z() <<
", " 2015 << activeDesignVariables_[3*cpI] <<
", " 2016 << activeDesignVariables_[3*cpI + 1] <<
", " 2017 << activeDesignVariables_[3*cpI + 2] <<
endl;
2025 parametricCoordinatesPtr_().write();
2031 cps_.writeEntry(
"controlPoints",
os);
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
void continuityRealatedConfinement()
Confine control point movement to maintain user-defined continuity.
vector volumeDerivativeU(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt u at point u,v,w.
label nWSymmetry() const
Get number of variables if CPs are moved symmetrically in W.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
void size(const label n)
Older name for setAddressableSize.
tmp< vectorField > getPointsInBox()
Get mesh points that reside within the control points box.
A class for handling file names.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
autoPtr< labelList > reverseMapPtr_
Map of mesh points to points-in-box.
tmp< tensorField > patchDxDb(const label patchI, const label cpI)
Get patch dx/db.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
void confineControlPoint(const label cpI)
Confine all three movements for a prescribed control point.
label start() const noexcept
Return start label of this patch in the polyMesh face list.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static const pointMesh & New(const polyMesh &mesh, Args &&... args)
Get existing or create a new MeshObject. Registered with typeName.
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
GeometricField< tensor, fvPatchField, volMesh > volTensorField
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
vectorField cps_
The volumetric B-Splines control points.
bool bound(vector &vec, scalar minValue=1e-7, scalar maxValue=0.999999) const
Bound components to certain limits.
const word dictName("faMeshDefinition")
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
void getIJK(label &i, label &j, label &k, const label cpID) const
Get I-J-K coordinates from control point ID.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
tmp< tensorField > dndbBasedSensitivities(const label patchI, const label cpI, bool DimensionedNormalSens=true)
Part of control point sensitivities related to the face normal variations.
autoPtr< labelList > mapPtr_
Map of points-in-box to mesh points.
tmp< tensorField > patchDxDbFace(const label patchI, const label cpI)
Get patch dx/db.
Differentiation of the mesh data structure.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
vectorField makeFaceCentresAndAreas_d(const pointField &p, const pointField &p_d)
Given a face and the points to be moved in the normal direction, find faceArea, faceCentre and unitVe...
label k
Boltzmann constant.
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
label nUSymmetry() const
Get number of variables if CPs are moved symmetrically in U.
List< labelList > labelListList
List of labelList.
const pointVectorField & getParametricCoordinates()
Get parametric coordinates.
NURBS3DVolume morpher. Includes support functions for gradient computations Base class providing supp...
Macros for easy insertion into run-time selection tables.
GeometricField< tensor, pointPatchField, pointMesh > pointTensorField
static autoPtr< NURBS3DVolume > New(const dictionary &dict, const fvMesh &mesh, bool computeParamCoors=true)
Return a reference to the selected NURBS model.
void confineInertControlPoints()
Deactivate control points if they affect no mesh point.
#define forAll(list, i)
Loop across all elements in list.
vectorField computeControlPointSensitivities(const pointVectorField &pointSens, const labelList &sensitivityPatchIDs)
Control point sensitivities computed using point-based surface sensitivities.
void determineActiveDesignVariablesAndPoints()
Create lists with active design variables and control points.
void makeFolders()
Create folders to store cps and derivatives.
pointField points(const UList< point > &pts) const
Return the points corresponding to this face.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
tensor JacobianUVW(const vector &u) const
Jacobian matrix wrt to the volume parametric coordinates.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
void writeParamCoordinates() const
Write parametric coordinates.
virtual bool writeData(Ostream &os) const
Write the control points to support restart.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void setControlPoints(const vectorField &newCps)
Set new control points.
tmp< vectorField > computeNewBoundaryPoints(const vectorField &controlPointsMovement, const labelList &patchesToBeMoved, const bool moveCPs=true)
Boundary mesh movement based on given control point movement.
void confineControlPointsDirections()
Confine movement in all control points for user-defined directions.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
vector volumeDerivativeV(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt v at point u,v,w.
void findPointsInBox(const vectorField &meshPoints)
Find points within control points box.
label getCPID(const label i, const label j, const label k) const
Get control point ID from its I-J-K coordinates.
dimensioned< tensor > dimensionedTensor
Dimensioned tensor obtained from generic dimensioned type.
void boundControlPointMovement(vectorField &controlPointsMovement) const
Bound control points movement in the boundary control points and in certain directions if needed...
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
#define DebugInfo
Report an information message using Foam::Info.
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross-product operators.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
label nVSymmetry() const
Get number of variables if CPs are moved symmetrically in V.
void writeCps(const fileName &="cpsFile", const bool transform=true) const
Write control points on a cartesian coordinates system for visualization.
Vector< label > nSymmetry() const
Get number of variables per direction, if CPs are moved symmetrically.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
Type gMax(const FieldField< Field, Type > &f)
OBJstream os(runTime.globalPath()/outputName)
tmp< vectorField > coordinates(const vectorField &uVector) const
Compute cartesian coordinates based on control points and parametric coordinates. ...
defineTypeNameAndDebug(combustionModel, 0)
void computeParametricCoordinates(const vectorField &points)
Compute parametric coordinates in order to match a given set of coordinates, based on the cps of the ...
tmp< vectorField > computeNewPoints(const vectorField &controlPointsMovement)
Mesh movement based on given control point movement.
const labelList & getMap()
Get map of points in box to mesh points.
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
void confineBoundaryControlPoints()
Confine movement in boundary control points if necessary.
PtrList< coordinateSystem > coordinates(solidRegions.size())
vector point
Point is a vector.
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Bound the given scalar field if it has gone unbounded.
#define WarningInFunction
Report a warning using Foam::Warning.
vector volumeDerivativeW(const scalar u, const scalar v, const scalar w) const
Volume point derivative wrt w at point u,v,w.
NURBS3DVolume(const dictionary &dict, const fvMesh &mesh, bool computeParamCoors=true)
Construct from dictionary.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
Automatically write from objectRegistry::writeObject()
const std::string patch
OpenFOAM patch number as a std::string.
scalar computeMaxBoundaryDisplacement(const vectorField &controlPointsMovement, const labelList &patchesToBeMoved)
Compute max. displacement at the boundary.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
const labelList & getReverseMap()
Get map of mesh points to points in box. Return -1 if point is outside the box.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
tmp< pointTensorField > getDxDb(const label cpI)
Get dxCartesiandb for a certain control point.
List< label > labelList
A List of labels.
tmp< volTensorField > getDxCellsDb(const label cpI)
Get dxCartesiandb for a certain control point on cells.
A class for managing temporary objects.
A patch is a list of labels that address the faces in the global face list.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Tensor of scalars, i.e. Tensor<scalar>.
#define FatalIOErrorInLookup(ios, lookupTag, lookupName, lookupTable)
Report an error message using Foam::FatalIOError.
List< bool > boolList
A List of bools.
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
static autoPtr< controlPointsDefinition > New(NURBS3DVolume &box)
Return a reference to the selected controlPointsDefinition model.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
T getOrDefaultCompat(const word &keyword, std::initializer_list< std::pair< const char *, int >> compat, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value using any compatibility names if needed...
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)
scalar volumeDerivativeCP(const vector &u, const label cpI) const
Volume point derivative wrt to control point cpI at point u,v,w.