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
1346 auto& dndbSens = tdndbSens.ref();
1350 const label patchStart = ppatch.
start();
1351 const labelList& reverseMap = reverseMapPtr_();
1354 const vectorField& parametricCoordinates = getParametricCoordinates();
1359 const face& fGlobal = mesh_.faces()[fI + patchStart];
1366 const label whichPointInBox = reverseMap[
globalIndex];
1369 if (whichPointInBox != -1)
1373 facePointDerivs[pI] =
1391 if (DimensionedNormalSens)
1393 dndbSens[fI] = dNdbSens[1];
1397 dndbSens[fI] = dNdbSens[2];
1412 const vectorField& parametricCoordinates = getParametricCoordinates();
1415 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1420 auto& dxdb = tdxdb.ref();
1424 const label globalIndex = meshPoints[pI];
1425 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1428 if (whichPointInBox != -1)
1431 transformationTensorDxDb(globalIndex)
1434 parametricCoordinates[globalIndex],
1451 const vectorField& parametricCoordinates = getParametricCoordinates();
1454 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1455 const label patchStart =
patch.start();
1459 auto& dxdb = tdxdb.ref();
1462 deltaBoundary deltaBound(mesh_);
1466 const face& fGlobal = mesh_.faces()[fI + patchStart];
1467 const pointField facePoints = fGlobal.points(mesh_.points());
1472 const label globalIndex = fGlobal[pI];
1473 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1476 if (whichPointInBox != -1)
1480 facePointDerivs[pI] =
1481 transformationTensorDxDb(globalIndex)
1484 parametricCoordinates[globalIndex],
1491 deltaBound.makeFaceCentresAndAreas_d
1507 const label
nPoints = mapPtr_().size();
1509 auto&
points = tpoints.ref();
1513 const label globalPI = mapPtr_()[pI];
1526 const label degreeU = basisU_.degree();
1527 const label degreeV = basisV_.degree();
1528 const label degreeW = basisW_.degree();
1530 const label nCPsU = basisU_.nCPs();
1531 const label nCPsV = basisV_.nCPs();
1532 const label nCPsW = basisW_.nCPs();
1534 const scalar u = uVector.x();
1535 const scalar v = uVector.y();
1536 const scalar w = uVector.z();
1539 for (label iCPw = 0; iCPw < nCPsW; iCPw++)
1541 const scalar basisW(basisW_.basisValue(iCPw, degreeW, w));
1542 for (label iCPv = 0; iCPv < nCPsV; iCPv++)
1544 const scalar basisVW =
1545 basisW*basisV_.basisValue(iCPv, degreeV, v);
1546 for (label iCPu = 0; iCPu < nCPsU; iCPu++)
1549 cps_[getCPID(iCPu, iCPv, iCPw)]
1550 *basisU_.basisValue(iCPu, degreeU, u)
1566 const vectorField& paramCoors = getParametricCoordinates();
1570 cps_ += controlPointsMovement;
1571 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1574 tmp<vectorField> tparameterizedPoints =
coordinates(paramCoors);
1575 const vectorField& parameterizedPoints = tparameterizedPoints();
1579 auto& newPoints = tnewPoints.ref();
1582 forAll(parameterizedPoints, pI)
1584 newPoints[map[pI]] = transformPointToCartesian(parameterizedPoints[pI]);
1588 updateLocalCoordinateSystem(newPoints);
1590 <<
"Max mesh movement equal to " 1601 const bool updateCPs
1605 const vectorField& paramCoors = getParametricCoordinates();
1608 cps_ += controlPointsMovement;
1612 writeCps(
"cpsBsplines"+mesh_.time().timeName());
1617 auto& newPoints = tnewPoints.ref();
1620 for (
const label patchI : patchesToBeMoved)
1622 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1625 for (
const label globalIndex : meshPoints)
1627 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1630 if (whichPointInBox != -1)
1632 newPoints[globalIndex] =
1633 transformPointToCartesian
1637 paramCoors[globalIndex]
1647 updateLocalCoordinateSystem(newPoints);
1652 cps_ -= controlPointsMovement;
1656 <<
"Max mesh movement equal to " 1670 const label nCPsU = basisU_.nCPs();
1671 const label nCPsV = basisV_.nCPs();
1673 return k*nCPsU*nCPsV + j*nCPsU + i;
1685 const label nCPsU = basisU_.nCPs();
1686 const label nCPsV = basisV_.nCPs();
1687 k = cpID/(nCPsU*nCPsV);
1688 const label inKplaneID = (cpID -
k*(nCPsU*nCPsV));
1689 j = inKplaneID/nCPsU;
1690 i = inKplaneID - j*nCPsU;
1696 if (cps_.size() != newCps.
size())
1699 <<
"Attempting to replace control points with a set of " 1712 forAll(controlPointsMovement, cpI)
1714 if (!activeDesignVariables_[3*cpI])
1716 controlPointsMovement[cpI].x() =
Zero;
1718 if (!activeDesignVariables_[3*cpI + 1])
1720 controlPointsMovement[cpI].y() =
Zero;
1722 if (!activeDesignVariables_[3*cpI + 2])
1724 controlPointsMovement[cpI].z() =
Zero;
1739 const vectorField& paramCoors = getParametricCoordinates();
1741 cps_ += controlPointsMovement;
1743 scalar maxDisplacement(
Zero);
1744 for (
const label patchI : patchesToBeMoved)
1746 const polyPatch&
patch = mesh_.boundaryMesh()[patchI];
1749 for (
const label globalIndex : meshPoints)
1751 const label whichPointInBox = reverseMapPtr_()[globalIndex];
1754 if (whichPointInBox != -1)
1757 transformPointToCartesian
1761 paramCoors[globalIndex]
1768 mag(newPoint - mesh_.points()[globalIndex])
1776 return maxDisplacement;
1784 findPointsInBox(localSystemCoordinates_);
1786 tmp<vectorField> pointsInBox
1799 findPointsInBox(localSystemCoordinates_);
1808 if (!reverseMapPtr_)
1810 findPointsInBox(localSystemCoordinates_);
1813 return reverseMapPtr_();
1820 if (!parametricCoordinatesPtr_)
1826 findPointsInBox(localSystemCoordinates_);
1828 computeParametricCoordinates(getPointsInBox()());
1831 return parametricCoordinatesPtr_();
1838 const vectorField& parametricCoordinates = getParametricCoordinates();
1841 tmp<pointTensorField> tDxDb
1848 mesh_.time().timeName(),
1863 for (
const label globalIndex : map)
1866 transformationTensorDxDb(globalIndex)
1869 parametricCoordinates[globalIndex],
1884 const vectorField& parametricCoordinates = getParametricCoordinates();
1896 deltaBoundary deltaBound(mesh_);
1902 for (
const label globalIndex : map)
1905 transformationTensorDxDb(globalIndex)
1908 parametricCoordinates[globalIndex],
1911 const labelList& pointCellsI = pointCells[globalIndex];
1912 tmp<tensorField> tC_d = deltaBound.cellCenters_d(globalIndex);
1916 const label cellI = pointCellsI[cI];
1917 DxDb[cellI] += C_d[cI] & pointDxDb;
1922 forAll(mesh_.boundary(), pI)
1924 const fvPatch&
patch = mesh_.boundary()[pI];
1925 if (!isA<coupledFvPatch>(
patch))
1927 DxDb.boundaryFieldRef()[pI] = patchDxDbFace(pI, cpI);
1932 DxDb.correctBoundaryConditions();
1940 label nU(basisU_.nCPs());
1941 return label(nU % 2 == 0 ? 0.5*nU : 0.5*(nU - 1) + 1);
1947 label nV(basisV_.nCPs());
1948 return label(nV % 2 == 0 ? 0.5*nV : 0.5*(nV - 1) + 1);
1954 label nW(basisW_.nCPs());
1955 return label(nW % 2 == 0 ? 0.5*nW : 0.5*(nW - 1) + 1);
1961 return Vector<label>(nUSymmetry(), nVSymmetry(), nWSymmetry());
1971 const label nCPsU = basisU_.nCPs();
1972 const label nCPsV = basisV_.nCPs();
1977 forAll(cpsInCartesian, cpI)
1979 cpsInCartesian[cpI] = transformPointToCartesian(cps_[cpI]);
1983 Info<<
"Writing control point positions to file" <<
endl;
1987 OFstream cpsFile(
"optimisation"/cpsFolder_/name_ + baseName +
".csv");
1990 <<
"\"Points : 0\", \"Points : 1\", \"Points : 2\"," 1991 <<
"\"i\", \"j\", \"k\"," 1992 <<
"\"active : 0\", \"active : 1\", \"active : 2\"" <<
endl;
1994 forAll(cpsInCartesian, cpI)
1996 const label iCPw = cpI/label(nCPsU*nCPsV);
1997 const label iCPv = (cpI - iCPw*nCPsU*nCPsV)/nCPsU;
1998 const label iCPu = (cpI - iCPw*nCPsU*nCPsV - iCPv*nCPsU);
2001 << cpsInCartesian[cpI].x() <<
", " 2002 << cpsInCartesian[cpI].y() <<
", " 2003 << cpsInCartesian[cpI].z() <<
", " 2007 << activeDesignVariables_[3*cpI] <<
", " 2008 << activeDesignVariables_[3*cpI + 1] <<
", " 2009 << activeDesignVariables_[3*cpI + 2] <<
endl;
2017 parametricCoordinatesPtr_().write();
2023 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 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.
static tmp< GeometricField< tensor, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< tensor >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
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.
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.
Do not request registration (bool: false)
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Combines List elements. 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.