1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "motionSmoother.H"
28 #include "twoDPointCorrector.H"
31 #include "fixedValuePointPatchFields.H"
32 #include "pointConstraint.H"
33 #include "syncTools.H"
34 #include "meshTools.H"
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
42 defineTypeNameAndDebug(motionSmoother, 0);
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 // From pointPatchInterpolation
50 void Foam::motionSmoother::makePatchPatchAddressing()
54 Pout<< "motionSmoother::makePatchPatchAddressing() : "
55 << "constructing boundary addressing"
59 const polyBoundaryMesh& bm = mesh_.boundaryMesh();
60 const pointBoundaryMesh& pbm = pMesh_.boundary();
62 // first count the total number of patch-patch points
64 label nPatchPatchPoints = 0;
68 if(!isA<emptyPolyPatch>(bm[patchi]))
70 nPatchPatchPoints += bm[patchi].boundaryPoints().size();
75 // Go through all patches and mark up the external edge points
76 Map<label> patchPatchPointSet(2*nPatchPatchPoints);
78 labelList patchPatchPoints(nPatchPatchPoints);
80 List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
86 if(!isA<emptyPolyPatch>(bm[patchi]))
88 const labelList& bp = bm[patchi].boundaryPoints();
89 const labelList& meshPoints = bm[patchi].meshPoints();
93 label ppp = meshPoints[bp[pointi]];
95 Map<label>::iterator iter = patchPatchPointSet.find(ppp);
97 if (iter == patchPatchPointSet.end())
99 patchPatchPointSet.insert(ppp, pppi);
100 patchPatchPoints[pppi] = ppp;
101 pbm[patchi].applyConstraint
104 patchPatchPointConstraints[pppi]
110 pbm[patchi].applyConstraint
113 patchPatchPointConstraints[iter()]
120 nPatchPatchPoints = pppi;
121 patchPatchPoints.setSize(nPatchPatchPoints);
122 patchPatchPointConstraints.setSize(nPatchPatchPoints);
124 patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
125 patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
127 label nConstraints = 0;
129 forAll(patchPatchPointConstraints, i)
131 if (patchPatchPointConstraints[i].first() != 0)
133 patchPatchPointConstraintPoints_[nConstraints] =
136 patchPatchPointConstraintTensors_[nConstraints] =
137 patchPatchPointConstraints[i].constraintTransformation();
143 patchPatchPointConstraintPoints_.setSize(nConstraints);
144 patchPatchPointConstraintTensors_.setSize(nConstraints);
149 OFstream str(mesh_.db().path()/"constraintPoints.obj");
151 Pout<< "Dumping " << patchPatchPointConstraintPoints_.size()
152 << " constraintPoints to " << str.name() << endl;
153 forAll(patchPatchPointConstraintPoints_, i)
155 label pointI = patchPatchPointConstraintPoints_[i];
157 meshTools::writeOBJ(str, mesh_.points()[pointI]);
160 Pout<< "motionSmoother::makePatchPatchAddressing() : "
161 << "finished constructing boundary addressing"
167 void Foam::motionSmoother::checkFld(const pointScalarField& fld)
171 const scalar val = fld[pointI];
173 if ((val > -GREAT) && (val < GREAT))
177 FatalErrorIn("motionSmoother::checkFld")
178 << "Problem : point:" << pointI << " value:" << val
179 << abort(FatalError);
185 Foam::labelHashSet Foam::motionSmoother::getPoints
187 const labelHashSet& faceLabels
190 labelHashSet usedPoints(mesh_.nPoints()/100);
192 forAllConstIter(labelHashSet, faceLabels, iter)
194 const face& f = mesh_.faces()[iter.key()];
198 usedPoints.insert(f[fp]);
206 // Smooth on selected points (usually patch points)
207 void Foam::motionSmoother::minSmooth
209 const PackedList<1>& isAffectedPoint,
210 const labelList& meshPoints,
211 const pointScalarField& fld,
212 pointScalarField& newFld
215 tmp<pointScalarField> tavgFld = avg
218 scalarField(mesh_.nEdges(), 1.0), // uniform weighting
219 false // fld is not position.
221 const pointScalarField& avgFld = tavgFld();
223 forAll(meshPoints, i)
225 label pointI = meshPoints[i];
226 if (isAffectedPoint.get(pointI) == 1)
231 0.5*fld[pointI] + 0.5*avgFld[pointI]
236 newFld.correctBoundaryConditions();
237 applyCornerConstraints(newFld);
241 // Smooth on all internal points
242 void Foam::motionSmoother::minSmooth
244 const PackedList<1>& isAffectedPoint,
245 const pointScalarField& fld,
246 pointScalarField& newFld
249 tmp<pointScalarField> tavgFld = avg
252 scalarField(mesh_.nEdges(), 1.0), // uniform weighting
253 false // fld is not position.
255 const pointScalarField& avgFld = tavgFld();
259 if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
264 0.5*fld[pointI] + 0.5*avgFld[pointI]
269 newFld.correctBoundaryConditions();
270 applyCornerConstraints(newFld);
274 // Scale on selected points
275 void Foam::motionSmoother::scaleField
277 const labelHashSet& pointLabels,
279 pointScalarField& fld
282 forAllConstIter(labelHashSet, pointLabels, iter)
284 if (isInternalPoint(iter.key()))
286 fld[iter.key()] *= scale;
289 fld.correctBoundaryConditions();
290 applyCornerConstraints(fld);
294 // Scale on selected points (usually patch points)
295 void Foam::motionSmoother::scaleField
297 const labelList& meshPoints,
298 const labelHashSet& pointLabels,
300 pointScalarField& fld
303 forAll(meshPoints, i)
305 label pointI = meshPoints[i];
307 if (pointLabels.found(pointI))
309 fld[pointI] *= scale;
315 bool Foam::motionSmoother::isInternalPoint(const label pointI) const
317 return isInternalPoint_.get(pointI) == 1;
321 void Foam::motionSmoother::getAffectedFacesAndPoints
323 const label nPointIter,
324 const faceSet& wrongFaces,
326 labelList& affectedFaces,
327 PackedList<1>& isAffectedPoint
330 isAffectedPoint.setSize(mesh_.nPoints());
333 faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
335 // Find possible points influenced by nPointIter iterations of
336 // scaling and smoothing by doing pointCellpoint walk.
337 // Also update faces-to-be-checked to extend one layer beyond the points
338 // that will get updated.
340 for (label i = 0; i < nPointIter; i++)
342 pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
344 forAllConstIter(pointSet, nbrPoints, iter)
346 const labelList& pCells = mesh_.pointCells()[iter.key()];
348 forAll(pCells, pCellI)
350 const cell& cFaces = mesh_.cells()[pCells[pCellI]];
352 forAll(cFaces, cFaceI)
354 nbrFaces.insert(cFaces[cFaceI]);
358 nbrFaces.sync(mesh_);
360 if (i == nPointIter - 2)
362 forAllConstIter(faceSet, nbrFaces, iter)
364 const face& f = mesh_.faces()[iter.key()];
367 isAffectedPoint.set(f[fp], 1);
373 affectedFaces = nbrFaces.toc();
377 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
379 Foam::motionSmoother::motionSmoother
383 indirectPrimitivePatch& pp,
384 const labelList& adaptPatchIDs,
385 const dictionary& paramDict
391 adaptPatchIDs_(adaptPatchIDs),
392 paramDict_(paramDict),
398 mesh_.time().timeName(),
410 mesh_.time().timeName(),
416 dimensionedScalar("scale", dimless, 1.0)
418 oldPoints_(mesh_.points()),
419 isInternalPoint_(mesh_.nPoints(), 1),
420 twoDCorrector_(mesh_)
426 Foam::motionSmoother::motionSmoother
429 indirectPrimitivePatch& pp,
430 const labelList& adaptPatchIDs,
431 const pointVectorField& displacement,
432 const dictionary& paramDict
436 pMesh_(const_cast<pointMesh&>(displacement.mesh())),
438 adaptPatchIDs_(adaptPatchIDs),
439 paramDict_(paramDict),
445 mesh_.time().timeName(),
457 mesh_.time().timeName(),
463 dimensionedScalar("scale", dimless, 1.0)
465 oldPoints_(mesh_.points()),
466 isInternalPoint_(mesh_.nPoints(), 1),
467 twoDCorrector_(mesh_)
473 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
475 Foam::motionSmoother::~motionSmoother()
479 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
481 const Foam::polyMesh& Foam::motionSmoother::mesh() const
487 const Foam::pointMesh& Foam::motionSmoother::pMesh() const
493 const Foam::indirectPrimitivePatch& Foam::motionSmoother::patch() const
499 const Foam::labelList& Foam::motionSmoother::adaptPatchIDs() const
501 return adaptPatchIDs_;
505 Foam::pointVectorField& Foam::motionSmoother::displacement()
507 return displacement_;
511 const Foam::pointVectorField& Foam::motionSmoother::displacement() const
513 return displacement_;
517 const Foam::pointScalarField& Foam::motionSmoother::scale() const
523 const Foam::pointField& Foam::motionSmoother::oldPoints() const
529 void Foam::motionSmoother::correct()
531 oldPoints_ = mesh_.points();
535 // No need to update twoDmotion corrector since only holds edge labels
536 // which will remain the same as before. So unless the mesh was distorted
537 // severely outside of motionSmoother there will be no need.
541 void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
543 // See comment in .H file about shared points.
544 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
546 forAll(patches, patchI)
548 const polyPatch& pp = patches[patchI];
552 const labelList& meshPoints = pp.meshPoints();
554 forAll(meshPoints, i)
556 displacement_[meshPoints[i]] = vector::zero;
561 const labelList& ppMeshPoints = pp_.meshPoints();
563 // Set internal point data from displacement on combined patch points.
564 forAll(ppMeshPoints, patchPointI)
566 displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
569 // Adapt the fixedValue bc's (i.e. copy internal point data to
570 // boundaryField for all affected patches)
571 forAll(adaptPatchIDs_, i)
573 label patchI = adaptPatchIDs_[i];
575 displacement_.boundaryField()[patchI] ==
576 displacement_.boundaryField()[patchI].patchInternalField();
579 // Make consistent with non-adapted bc's by evaluating those now and
580 // resetting the displacement from the values.
581 // Note that we're just doing a correctBoundaryConditions with
582 // fixedValue bc's first.
583 labelHashSet adaptPatchSet(adaptPatchIDs_);
585 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
587 forAll(patchSchedule, patchEvalI)
589 label patchI = patchSchedule[patchEvalI].patch;
591 if (!adaptPatchSet.found(patchI))
593 if (patchSchedule[patchEvalI].init)
595 displacement_.boundaryField()[patchI]
596 .initEvaluate(Pstream::scheduled);
600 displacement_.boundaryField()[patchI]
601 .evaluate(Pstream::scheduled);
606 // Multi-patch constraints
607 applyCornerConstraints(displacement_);
609 // Correct for problems introduced by corner constraints
610 syncTools::syncPointList
614 maxMagEqOp(), // combine op
615 vector::zero, // null value
616 false // no separation
619 // Adapt the fixedValue bc's (i.e. copy internal point data to
620 // boundaryField for all affected patches) to take the changes caused
621 // by multi-corner constraints into account.
622 forAll(adaptPatchIDs_, i)
624 label patchI = adaptPatchIDs_[i];
626 displacement_.boundaryField()[patchI] ==
627 displacement_.boundaryField()[patchI].patchInternalField();
632 OFstream str(mesh_.db().path()/"changedPoints.obj");
634 forAll(ppMeshPoints, patchPointI)
636 const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
638 if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
640 const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
642 meshTools::writeOBJ(str, pt);
644 //Pout<< "Point:" << pt
645 // << " oldDisp:" << patchDisp[patchPointI]
646 // << " newDisp:" << newDisp << endl;
649 Pout<< "Written " << nVerts << " points that are changed to file "
650 << str.name() << endl;
653 // Now reset input displacement
654 forAll(ppMeshPoints, patchPointI)
656 patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
661 // correctBoundaryConditions with fixedValue bc's first.
662 void Foam::motionSmoother::correctBoundaryConditions
664 pointVectorField& displacement
667 labelHashSet adaptPatchSet(adaptPatchIDs_);
669 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
671 // 1. evaluate on adaptPatches
672 forAll(patchSchedule, patchEvalI)
674 label patchI = patchSchedule[patchEvalI].patch;
676 if (adaptPatchSet.found(patchI))
678 if (patchSchedule[patchEvalI].init)
680 displacement.boundaryField()[patchI]
681 .initEvaluate(Pstream::blocking);
685 displacement.boundaryField()[patchI]
686 .evaluate(Pstream::blocking);
692 // 2. evaluate on non-AdaptPatches
693 forAll(patchSchedule, patchEvalI)
695 label patchI = patchSchedule[patchEvalI].patch;
697 if (!adaptPatchSet.found(patchI))
699 if (patchSchedule[patchEvalI].init)
701 displacement.boundaryField()[patchI]
702 .initEvaluate(Pstream::blocking);
706 displacement.boundaryField()[patchI]
707 .evaluate(Pstream::blocking);
712 // Multi-patch constraints
713 applyCornerConstraints(displacement);
715 // Correct for problems introduced by corner constraints
716 syncTools::syncPointList
720 maxMagEqOp(), // combine op
721 vector::zero, // null value
722 false // no separation
727 Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
729 pointField& newPoints
732 // Correct for 2-D motion
733 if (twoDCorrector_.required())
735 Info<< "Correct-ing 2-D mesh motion";
737 if (mesh_.globalData().parallel())
739 WarningIn("motionSmoother::movePoints(pointField& newPoints)")
740 << "2D mesh-motion probably not correct in parallel" << endl;
743 // We do not want to move 3D planes so project all points onto those
744 const pointField& oldPoints = mesh_.points();
745 const edgeList& edges = mesh_.edges();
747 const labelList& neIndices = twoDCorrector().normalEdgeIndices();
748 const vector& pn = twoDCorrector().planeNormal();
752 const edge& e = edges[neIndices[i]];
754 point& pStart = newPoints[e.start()];
756 pStart += pn*(pn & (oldPoints[e.start()] - pStart));
758 point& pEnd = newPoints[e.end()];
760 pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
763 // Correct tangentially
764 twoDCorrector_.correctPoints(newPoints);
765 Info<< " ...done" << endl;
770 Pout<< "motionSmoother::movePoints : testing sync of newPoints."
775 minEqOp<point>(), // combine op
776 vector(GREAT,GREAT,GREAT), // null
781 tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
783 //!!! Workaround for movePoints bug
784 const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh()).movePoints(newPoints);
786 pp_.movePoints(mesh_.points());
792 Foam::scalar Foam::motionSmoother::setErrorReduction
794 const scalar errorReduction
797 scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
799 paramDict_.remove("errorReduction");
800 paramDict_.add("errorReduction", errorReduction);
802 return oldErrorReduction;
806 bool Foam::motionSmoother::scaleMesh
808 labelList& checkFaces,
809 const bool smoothMesh,
810 const label nAllowableErrors
813 List<labelPair> emptyBaffles;
824 bool Foam::motionSmoother::scaleMesh
826 labelList& checkFaces,
827 const List<labelPair>& baffles,
828 const bool smoothMesh,
829 const label nAllowableErrors
832 if (!smoothMesh && adaptPatchIDs_.size() == 0)
834 FatalErrorIn("motionSmoother::scaleMesh(const bool")
835 << "You specified both no movement on the internal mesh points"
836 << " (smoothMesh = false)" << nl
837 << "and no movement on the patch (adaptPatchIDs is empty)" << nl
838 << "Hence nothing to adapt."
844 // Had a problem with patches moved non-synced. Check transformations.
845 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
847 Pout<< "Entering scaleMesh : coupled patches:" << endl;
848 forAll(patches, patchI)
850 if (patches[patchI].coupled())
852 const coupledPolyPatch& pp =
853 refCast<const coupledPolyPatch>(patches[patchI]);
855 Pout<< '\t' << patchI << '\t' << pp.name()
856 << " parallel:" << pp.parallel()
857 << " separated:" << pp.separated()
858 << " forwardT:" << pp.forwardT().size()
864 const scalar errorReduction =
865 readScalar(paramDict_.lookup("errorReduction"));
866 const label nSmoothScale =
867 readLabel(paramDict_.lookup("nSmoothScale"));
870 // Note: displacement_ should already be synced already from setDisplacement
871 // but just to make sure.
872 syncTools::syncPointList
877 vector::zero, // null value
878 false // no separation
881 // Set newPoints as old + scale*displacement
882 pointField newPoints;
884 // Create overall displacement with same b.c.s as displacement_
885 pointVectorField totalDisplacement
890 mesh_.time().timeName(),
896 scale_*displacement_,
897 displacement_.boundaryField().types()
899 correctBoundaryConditions(totalDisplacement);
903 Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
908 vector::zero, // null value
913 newPoints = oldPoints_ + totalDisplacement.internalField();
916 Info<< "Moving mesh using diplacement scaling :"
917 << " min:" << gMin(scale_.internalField())
918 << " max:" << gMax(scale_.internalField())
923 movePoints(newPoints);
925 // Check. Returns parallel number of incorrect faces.
926 faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
927 checkMesh(false, mesh_, paramDict_, checkFaces, baffles, wrongFaces);
929 if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
935 // Sync across coupled faces by extending the set.
936 wrongFaces.sync(mesh_);
939 // if errorReduction is set to zero, extend wrongFaces
940 // to face-Cell-faces to ensure quick return to previously valid mesh
942 if (mag(errorReduction) < SMALL)
944 labelHashSet newWrongFaces(wrongFaces);
945 forAllConstIter(labelHashSet, wrongFaces, iter)
947 label own = mesh_.faceOwner()[iter.key()];
948 const cell& ownFaces = mesh_.cells()[own];
950 forAll(ownFaces, cfI)
952 newWrongFaces.insert(ownFaces[cfI]);
955 if (iter.key() < mesh_.nInternalFaces())
957 label nei = mesh_.faceNeighbour()[iter.key()];
958 const cell& neiFaces = mesh_.cells()[nei];
960 forAll(neiFaces, cfI)
962 newWrongFaces.insert(neiFaces[cfI]);
966 wrongFaces.transfer(newWrongFaces);
967 wrongFaces.sync(mesh_);
971 // Find out points used by wrong faces and scale displacement.
972 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
974 pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
975 usedPoints.sync(mesh_);
979 // Grow a few layers to determine
980 // - points to be smoothed
981 // - faces to be checked in next iteration
982 PackedList<1> isAffectedPoint(mesh_.nPoints(), 0);
983 getAffectedFacesAndPoints
985 nSmoothScale, // smoothing iterations
986 wrongFaces, // error faces
993 Pout<< "Faces in error:" << wrongFaces.size()
994 << " with points:" << usedPoints.size()
998 if (adaptPatchIDs_.size() != 0)
1000 // Scale conflicting patch points
1001 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
1005 // Scale conflicting internal points
1006 scaleField(usedPoints, errorReduction, scale_);
1009 for (label i = 0; i < nSmoothScale; i++)
1011 if (adaptPatchIDs_.size() != 0)
1013 // Smooth patch values
1014 pointScalarField oldScale(scale_);
1026 // Smooth internal values
1027 pointScalarField oldScale(scale_);
1028 minSmooth(isAffectedPoint, oldScale, scale_);
1033 syncTools::syncPointList
1038 -GREAT, // null value
1039 false // no separation
1045 Pout<< "scale_ after smoothing :"
1046 << " min:" << Foam::gMin(scale_)
1047 << " max:" << Foam::gMax(scale_)
1056 void Foam::motionSmoother::updateMesh()
1058 const pointBoundaryMesh& patches = pMesh_.boundary();
1060 // Check whether displacement has fixed value b.c. on adaptPatchID
1061 forAll(adaptPatchIDs_, i)
1063 label patchI = adaptPatchIDs_[i];
1067 !isA<fixedValuePointPatchVectorField>
1069 displacement_.boundaryField()[patchI]
1075 "motionSmoother::motionSmoother"
1076 ) << "Patch " << patches[patchI].name()
1077 << " has wrong boundary condition "
1078 << displacement_.boundaryField()[patchI].type()
1079 << " on field " << displacement_.name() << nl
1080 << "Only type allowed is "
1081 << fixedValuePointPatchVectorField::typeName
1082 << exit(FatalError);
1087 // Determine internal points. Note that for twoD there are no internal
1088 // points so we use the points of adaptPatchIDs instead
1090 twoDCorrector_.updateMesh();
1092 const labelList& meshPoints = pp_.meshPoints();
1094 forAll(meshPoints, i)
1096 isInternalPoint_.set(meshPoints[i], 0);
1099 // Calculate master edge addressing
1100 isMasterEdge_ = syncTools::getMasterEdges(mesh_);
1102 makePatchPatchAddressing();
1106 // Specialisation of applyCornerConstraints for scalars because
1107 // no constraint need be applied
1109 void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
1111 GeometricField<scalar, pointPatchField, pointMesh>& pf
1116 // ************************************************************************* //