49 void Foam::motionSmoother::makePatchPatchAddressing()
53 Pout<<
"motionSmoother::makePatchPatchAddressing() : "
54 <<
"constructing boundary addressing"
59 const pointBoundaryMesh& pbm = pMesh_.
boundary();
63 label nPatchPatchPoints = 0;
67 if(!isA<emptyPolyPatch>(bm[
patchi]))
69 nPatchPatchPoints += bm[
patchi].boundaryPoints().size();
75 Map<label> patchPatchPointSet(2*nPatchPatchPoints);
77 labelList patchPatchPoints(nPatchPatchPoints);
79 List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
85 if(!isA<emptyPolyPatch>(bm[patchi]))
92 label ppp = meshPoints[bp[pointi]];
96 if (iter == patchPatchPointSet.end())
98 patchPatchPointSet.insert(ppp, pppi);
99 patchPatchPoints[pppi] = ppp;
100 pbm[
patchi].applyConstraint
103 patchPatchPointConstraints[pppi]
109 pbm[
patchi].applyConstraint
112 patchPatchPointConstraints[iter()]
119 nPatchPatchPoints = pppi;
120 patchPatchPoints.setSize(nPatchPatchPoints);
121 patchPatchPointConstraints.setSize(nPatchPatchPoints);
123 patchPatchPointConstraintPoints_.
setSize(nPatchPatchPoints);
124 patchPatchPointConstraintTensors_.
setSize(nPatchPatchPoints);
126 label nConstraints = 0;
128 forAll(patchPatchPointConstraints, i)
130 if (patchPatchPointConstraints[i].first() != 0)
132 patchPatchPointConstraintPoints_[nConstraints] =
135 patchPatchPointConstraintTensors_[nConstraints] =
136 patchPatchPointConstraints[i].constraintTransformation();
142 patchPatchPointConstraintPoints_.
setSize(nConstraints);
143 patchPatchPointConstraintTensors_.
setSize(nConstraints);
148 OFstream str(mesh_.
db().
path()/
"constraintPoints.obj");
150 Pout<<
"Dumping " << patchPatchPointConstraintPoints_.
size()
151 <<
" constraintPoints to " << str.name() <<
endl;
152 forAll(patchPatchPointConstraintPoints_, i)
154 label pointI = patchPatchPointConstraintPoints_[i];
159 Pout<<
"motionSmoother::makePatchPatchAddressing() : "
160 <<
"finished constructing boundary addressing"
170 const scalar val = fld[pointI];
172 if ((val > -GREAT) && (val < GREAT))
177 <<
"Problem : point:" << pointI <<
" value:" << val
193 const face&
f = mesh_.faces()[iter.key()];
197 usedPoints.insert(f[fp]);
206 void Foam::motionSmoother::minSmooth
214 tmp<pointScalarField> tavgFld = avg
224 label pointI = meshPoints[i];
225 if (isAffectedPoint.get(pointI) == 1)
230 0.5*fld[pointI] + 0.5*avgFld[pointI]
235 newFld.correctBoundaryConditions();
236 applyCornerConstraints(newFld);
241 void Foam::motionSmoother::minSmooth
248 tmp<pointScalarField> tavgFld = avg
258 if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
263 0.5*fld[pointI] + 0.5*avgFld[pointI]
268 newFld.correctBoundaryConditions();
269 applyCornerConstraints(newFld);
274 void Foam::motionSmoother::scaleField
283 if (isInternalPoint(iter.key()))
285 fld[iter.key()] *= scale;
288 fld.correctBoundaryConditions();
289 applyCornerConstraints(fld);
294 void Foam::motionSmoother::scaleField
304 label pointI = meshPoints[i];
306 if (pointLabels.found(pointI))
308 fld[pointI] *= scale;
314 bool Foam::motionSmoother::isInternalPoint(
const label pointI)
const
316 return isInternalPoint_.get(pointI) == 1;
320 void Foam::motionSmoother::getAffectedFacesAndPoints
322 const label nPointIter,
323 const faceSet& wrongFaces,
329 isAffectedPoint.setSize(mesh_.nPoints());
332 faceSet nbrFaces(mesh_,
"checkFaces", wrongFaces);
339 for (label i = 0; i < nPointIter; i++)
341 pointSet nbrPoints(mesh_,
"grownPoints", getPoints(nbrFaces.toc()));
345 const labelList& pCells = mesh_.pointCells(iter.key());
349 const cell& cFaces = mesh_.cells()[pCells[pCellI]];
353 nbrFaces.insert(cFaces[cFaceI]);
357 nbrFaces.sync(mesh_);
359 if (i == nPointIter - 2)
363 const face& f = mesh_.faces()[iter.key()];
366 isAffectedPoint.set(f[fp], 1);
372 affectedFaces = nbrFaces.toc();
378 Foam::motionSmoother::motionSmoother
390 adaptPatchIDs_(adaptPatchIDs),
391 paramDict_(paramDict),
397 mesh_.time().timeName(),
409 mesh_.time().timeName(),
417 oldPoints_(mesh_.points()),
418 isInternalPoint_(mesh_.nPoints(), 1),
419 twoDCorrector_(mesh_)
425 Foam::motionSmoother::motionSmoother
435 pMesh_(const_cast<pointMesh&>(displacement.
mesh())),
437 adaptPatchIDs_(adaptPatchIDs),
438 paramDict_(paramDict),
444 mesh_.time().timeName(),
456 mesh_.time().timeName(),
464 oldPoints_(mesh_.points()),
465 isInternalPoint_(mesh_.nPoints(), 1),
466 twoDCorrector_(mesh_)
500 return adaptPatchIDs_;
512 return displacement_;
518 return displacement_;
536 oldPoints_ = mesh_.points();
566 const labelList& ppMeshPoints = pp_.meshPoints();
569 forAll(ppMeshPoints, patchPointI)
571 displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
578 label patchI = adaptPatchIDs_[i];
580 displacement_.boundaryField()[patchI] ==
581 displacement_.boundaryField()[patchI].patchInternalField();
590 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
592 forAll(patchSchedule, patchEvalI)
594 label patchI = patchSchedule[patchEvalI].patch;
596 if (!adaptPatchSet.
found(patchI))
598 if (patchSchedule[patchEvalI].init)
600 displacement_.boundaryField()[patchI]
605 displacement_.boundaryField()[patchI]
612 applyCornerConstraints(displacement_);
629 label patchI = adaptPatchIDs_[i];
631 displacement_.boundaryField()[patchI] ==
632 displacement_.boundaryField()[patchI].patchInternalField();
637 OFstream str(mesh_.db().path()/
"changedPoints.obj");
639 forAll(ppMeshPoints, patchPointI)
641 const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
643 if (
mag(newDisp-patchDisp[patchPointI]) > SMALL)
645 const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
654 Pout<<
"Written " << nVerts <<
" points that are changed to file "
659 forAll(ppMeshPoints, patchPointI)
661 patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
674 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
677 forAll(patchSchedule, patchEvalI)
679 label patchI = patchSchedule[patchEvalI].patch;
681 if (adaptPatchSet.
found(patchI))
683 if (patchSchedule[patchEvalI].init)
698 forAll(patchSchedule, patchEvalI)
700 label patchI = patchSchedule[patchEvalI].patch;
702 if (!adaptPatchSet.
found(patchI))
704 if (patchSchedule[patchEvalI].init)
718 applyCornerConstraints(displacement);
738 if (twoDCorrector_.required())
740 Info<<
"Correct-ing 2-D mesh motion";
742 if (mesh_.globalData().parallel())
744 WarningIn(
"motionSmoother::movePoints(pointField& newPoints)")
745 <<
"2D mesh-motion probably not correct in parallel" <<
endl;
750 const edgeList& edges = mesh_.edges();
752 const labelList& neIndices = twoDCorrector().normalEdgeIndices();
753 const vector& pn = twoDCorrector().planeNormal();
757 const edge&
e = edges[neIndices[i]];
761 pStart += pn*(pn & (oldPoints[e.
start()] - pStart));
765 pEnd += pn*(pn & (oldPoints[e.
end()] - pEnd));
769 twoDCorrector_.correctPoints(newPoints);
775 Pout<<
"motionSmoother::movePoints : testing sync of newPoints."
781 vector(GREAT,GREAT,GREAT),
783 1
E-6*mesh_.bounds().mag()
789 pp_.movePoints(mesh_.points());
797 const scalar errorReduction
800 scalar oldErrorReduction =
readScalar(paramDict_.lookup(
"errorReduction"));
802 paramDict_.remove(
"errorReduction");
803 paramDict_.add(
"errorReduction", errorReduction);
805 return oldErrorReduction;
812 const bool smoothMesh,
813 const label nAllowableErrors
831 const bool smoothMesh,
832 const label nAllowableErrors
853 const bool smoothMesh,
854 const label nAllowableErrors
857 if (!smoothMesh && adaptPatchIDs_.empty())
860 <<
"You specified both no movement on the internal mesh points"
861 <<
" (smoothMesh = false)" <<
nl
862 <<
"and no movement on the patch (adaptPatchIDs is empty)" <<
nl
863 <<
"Hence nothing to adapt."
872 Pout<<
"Entering scaleMesh : coupled patches:" <<
endl;
875 if (patches[patchI].coupled())
878 refCast<const coupledPolyPatch>(patches[patchI]);
880 Pout<<
'\t' << patchI <<
'\t' << pp.
name()
889 const scalar errorReduction =
891 const label nSmoothScale =
892 readLabel(paramDict.
lookup(
"nSmoothScale"));
915 mesh_.time().timeName(),
921 scale_*displacement_,
922 displacement_.boundaryField().types()
928 Pout<<
"scaleMesh : testing sync of totalDisplacement" <<
endl;
935 1
E-6*mesh_.bounds().mag()
939 newPoints = oldPoints_ + totalDisplacement.internalField();
942 Info<<
"Moving mesh using diplacement scaling :"
943 <<
" min:" <<
gMin(scale_.internalField())
944 <<
" max:" <<
gMax(scale_.internalField())
949 movePoints(newPoints);
952 faceSet wrongFaces(mesh_,
"wrongFaces", mesh_.nFaces()/100+100);
953 checkMesh(
false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
962 wrongFaces.sync(mesh_);
968 if (
mag(errorReduction) < SMALL)
973 label own = mesh_.faceOwner()[iter.key()];
974 const cell& ownFaces = mesh_.cells()[own];
978 newWrongFaces.
insert(ownFaces[cfI]);
981 if (iter.key() < mesh_.nInternalFaces())
983 label nei = mesh_.faceNeighbour()[iter.key()];
984 const cell& neiFaces = mesh_.cells()[nei];
988 newWrongFaces.
insert(neiFaces[cfI]);
992 wrongFaces.transfer(newWrongFaces);
993 wrongFaces.sync(mesh_);
1000 pointSet usedPoints(mesh_,
"usedPoints", getPoints(wrongFaces));
1001 usedPoints.
sync(mesh_);
1009 getAffectedFacesAndPoints
1019 Pout<<
"Faces in error:" << wrongFaces.size()
1020 <<
" with points:" << usedPoints.
size()
1024 if (adaptPatchIDs_.size())
1027 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
1032 scaleField(usedPoints, errorReduction, scale_);
1035 for (label i = 0; i < nSmoothScale; i++)
1037 if (adaptPatchIDs_.size())
1054 minSmooth(isAffectedPoint, oldScale, scale_);
1071 Pout<<
"scale_ after smoothing :"
1087 forAll(adaptPatchIDs_, i)
1089 label patchI = adaptPatchIDs_[i];
1093 !isA<fixedValuePointPatchVectorField>
1095 displacement_.boundaryField()[patchI]
1101 "motionSmoother::motionSmoother"
1102 ) <<
"Patch " << patches[patchI].name()
1103 <<
" has wrong boundary condition "
1104 << displacement_.boundaryField()[patchI].type()
1105 <<
" on field " << displacement_.name() <<
nl
1106 <<
"Only type allowed is "
1107 << fixedValuePointPatchVectorField::typeName
1116 twoDCorrector_.updateMesh();
1118 const labelList& meshPoints = pp_.meshPoints();
1122 isInternalPoint_.set(meshPoints[i], 0);
1128 makePatchPatchAddressing();
1135 void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
1137 GeometricField<scalar, pointPatchField, pointMesh>& pf