Go to the documentation of this file.
49 void Foam::motionSmootherAlgo::testSyncPositions
62 point(GREAT,GREAT,GREAT)
67 if (
mag(syncedFld[i] -
fld[i]) > maxMag)
70 <<
"On point " << i <<
" point:" <<
fld[i]
71 <<
" synchronised point:" << syncedFld[i]
82 const scalar val =
fld[pointi];
84 if ((val > -GREAT) && (val < GREAT))
89 <<
"Problem : point:" << pointi <<
" value:" << val
103 for (
const label
faceId : faceLabels)
117 const edgeList& edges = mesh_.edges();
120 auto& wght = twght.ref();
124 wght[edgei] = 1.0/(edges[edgei].mag(
points)+SMALL);
131 void Foam::motionSmootherAlgo::minSmooth
134 const bitSet& isAffectedPoint,
140 tmp<pointScalarField> tavgFld = avg
147 for (
const label pointi : meshPoints)
149 if (isAffectedPoint.test(pointi))
154 0.5*
fld[pointi] + 0.5*avgFld[pointi]
165 void Foam::motionSmootherAlgo::minSmooth
168 const bitSet& isAffectedPoint,
173 tmp<pointScalarField> tavgFld = avg
182 if (isAffectedPoint.test(pointi) && isInternalPoint_.test(pointi))
187 0.5*
fld[pointi] + 0.5*avgFld[pointi]
199 void Foam::motionSmootherAlgo::scaleField
208 if (isInternalPoint_.test(pointi))
210 fld[pointi] *= scale;
220 void Foam::motionSmootherAlgo::scaleField
228 for (
const label pointi : meshPoints)
232 fld[pointi] *= scale;
239 void Foam::motionSmootherAlgo::subtractField
248 if (isInternalPoint_.test(pointi))
260 void Foam::motionSmootherAlgo::subtractField
268 for (
const label pointi : meshPoints)
278 void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
280 const label nPointIter,
281 const faceSet& wrongFaces,
284 bitSet& isAffectedPoint
287 isAffectedPoint.reset();
288 isAffectedPoint.resize(mesh_.nPoints());
290 faceSet nbrFaces(mesh_,
"checkFaces", wrongFaces);
297 for (label i = 0; i < nPointIter; i++)
299 pointSet nbrPoints(mesh_,
"grownPoints", getPoints(nbrFaces));
301 for (
const label pointi : nbrPoints)
303 for (
const label celli : mesh_.pointCells(pointi))
305 nbrFaces.set(mesh_.cells()[celli]);
308 nbrFaces.sync(mesh_);
310 if (i == nPointIter - 2)
312 for (
const label facei : nbrFaces)
314 isAffectedPoint.set(mesh_.faces()[facei]);
319 affectedFaces = nbrFaces.toc();
325 Foam::motionSmootherAlgo::motionSmootherAlgo
341 displacement_(displacement),
343 oldPoints_(oldPoints),
344 adaptPatchIDs_(adaptPatchIDs),
345 paramDict_(paramDict),
347 isInternalPoint_(mesh_.nPoints(),
true)
381 return adaptPatchIDs_;
393 oldPoints_ = mesh_.points();
409 pointVectorField::Boundary& displacementBf =
414 for (
const label patchi : patchIDs)
416 displacementBf[patchi] ==
417 displacementBf[patchi].patchInternalField();
426 const lduSchedule& patchSchedule = displacement.mesh().globalData().
429 forAll(patchSchedule, patchEvalI)
431 const label patchi = patchSchedule[patchEvalI].patch;
433 if (!adaptPatchSet.
found(patchi))
435 if (patchSchedule[patchEvalI].init)
437 displacementBf[patchi]
442 displacementBf[patchi]
454 for (
const label patchi : patchIDs)
456 displacementBf[patchi] == displacementBf[patchi].patchInternalField();
463 setDisplacementPatchFields(adaptPatchIDs_, displacement_);
499 for (
const label pointi : cppMeshPoints)
501 if (isPatchPoint.test(pointi))
503 displacement[pointi] =
Zero;
510 forAll(ppMeshPoints, patchPointi)
512 displacement[ppMeshPoints[patchPointi]] = patchDisp[patchPointi];
528 setDisplacementPatchFields(patchIDs, displacement);
535 forAll(ppMeshPoints, patchPointi)
537 const vector& newDisp = displacement[ppMeshPoints[patchPointi]];
539 if (
mag(newDisp-patchDisp[patchPointi]) > SMALL)
550 Pout<<
"Written " << nVerts <<
" points that are changed to file "
555 forAll(ppMeshPoints, patchPointi)
557 patchDisp[patchPointi] = displacement[ppMeshPoints[patchPointi]];
564 setDisplacement(adaptPatchIDs_, pp_, patchDisp, displacement_);
576 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
578 pointVectorField::Boundary& displacementBf =
582 forAll(patchSchedule, patchEvalI)
584 const label patchi = patchSchedule[patchEvalI].patch;
586 if (adaptPatchSet.
found(patchi))
588 if (patchSchedule[patchEvalI].init)
590 displacementBf[patchi]
595 displacementBf[patchi]
603 forAll(patchSchedule, patchEvalI)
605 const label patchi = patchSchedule[patchEvalI].patch;
607 if (!adaptPatchSet.
found(patchi))
609 if (patchSchedule[patchEvalI].init)
611 displacementBf[patchi]
616 displacementBf[patchi]
644 Info<<
"Correcting 2-D mesh motion";
646 if (mesh_.globalData().parallel())
649 <<
"2D mesh-motion probably not correct in parallel" <<
endl;
654 const edgeList& edges = mesh_.edges();
659 for (
const label edgei : neIndices)
662 const edge&
e = edges[edgei];
664 point& pStart = newPoints[
e.start()];
666 pStart += pn*(pn & (oldPoints[
e.start()] - pStart));
668 point& pEnd = newPoints[
e.end()];
670 pEnd += pn*(pn & (oldPoints[
e.end()] - pEnd));
680 Pout<<
"motionSmootherAlgo::modifyMotionPoints :"
681 <<
" testing sync of newPoints."
683 testSyncPositions(newPoints, 1
e-6*mesh_.bounds().mag());
692 mesh_.clearTetBasePtIs();
693 pp_.movePoints(mesh_.points());
699 const scalar errorReduction
702 scalar old = paramDict_.get<scalar>(
"errorReduction");
704 paramDict_.remove(
"errorReduction");
705 paramDict_.add(
"errorReduction", errorReduction);
714 const bool smoothMesh,
715 const label nAllowableErrors
733 const bool smoothMesh,
734 const label nAllowableErrors
757 actualPatchTypes.
setSize(pbm.size());
760 actualPatchTypes[patchi] = pbm[patchi].type();
766 const pointVectorField::Boundary& pfld =
767 displacement_.boundaryField();
768 actualPatchFieldTypes.
setSize(pfld.size());
774 actualPatchFieldTypes[patchi] =
779 actualPatchFieldTypes[patchi] = pfld[patchi].type();
789 mesh_.time().timeName(),
795 scale_*displacement_,
796 actualPatchFieldTypes,
803 Pout<<
"scaleMesh : testing sync of totalDisplacement" <<
endl;
809 1
e-6*mesh_.bounds().mag()
816 modifyMotionPoints(tnewPoints.ref());
828 const bool smoothMesh,
829 const label nAllowableErrors
832 if (!smoothMesh && adaptPatchIDs_.empty())
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."
847 Pout<<
"Entering scaleMesh : coupled patches:" <<
endl;
853 refCast<const coupledPolyPatch>(
patches[patchi]);
855 Pout<<
'\t' << patchi <<
'\t' << pp.
name()
858 <<
" forwardT:" << pp.
forwardT().size()
864 const scalar errorReduction = get<scalar>
868 const label nSmoothScale = get<label>
883 Info<<
"Moving mesh using displacement scaling :"
884 <<
" min:" <<
gMin(scale_.primitiveField())
885 <<
" max:" <<
gMax(scale_.primitiveField())
892 mesh_.movePoints(newPoints);
897 faceSet wrongFaces(mesh_,
"wrongFaces", mesh_.nFaces()/100+100);
916 wrongFaces.sync(mesh_);
922 if (
mag(errorReduction) < SMALL)
925 for (
const label facei : wrongFaces)
927 const label own = mesh_.faceOwner()[facei];
928 const cell& ownFaces = mesh_.cells()[own];
930 newWrongFaces.
insert(ownFaces);
932 if (facei < mesh_.nInternalFaces())
934 const label nei = mesh_.faceNeighbour()[facei];
935 const cell& neiFaces = mesh_.cells()[nei];
937 newWrongFaces.
insert(neiFaces);
940 wrongFaces.transfer(newWrongFaces);
941 wrongFaces.sync(mesh_);
948 pointSet usedPoints(mesh_,
"usedPoints", getPoints(wrongFaces));
949 usedPoints.
sync(mesh_);
955 bitSet isAffectedPoint(mesh_.nPoints());
956 getAffectedFacesAndPoints
966 Pout<<
"Faces in error:" << wrongFaces.size()
967 <<
" with points:" << usedPoints.
size()
971 if (adaptPatchIDs_.size())
974 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
980 scaleField(usedPoints, errorReduction, scale_);
986 for (label i = 0; i < nSmoothScale; ++i)
988 if (adaptPatchIDs_.size())
1006 minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
1022 Pout<<
"scale_ after smoothing :"
1038 for (
const label patchi : adaptPatchIDs_)
1042 !isA<fixedValuePointPatchVectorField>
1044 displacement_.boundaryField()[patchi]
1050 <<
" has wrong boundary condition "
1051 << displacement_.boundaryField()[patchi].type()
1052 <<
" on field " << displacement_.name() <<
nl
1053 <<
"Only type allowed is "
1054 << fixedValuePointPatchVectorField::typeName
1063 isInternalPoint_.unset(pp_.meshPoints());
int debug
Static debugging option.
List< label > labelList
A List of labels.
vectorField pointField
pointField is a vectorField.
label size() const noexcept
The number of elements in table.
virtual const pointField & points() const
Return raw points.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
const pointBoundaryMesh & boundary() const
Return reference to boundary mesh.
const word & name() const
Return name.
List< edge > edgeList
A List of edges.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
void correctBoundaryConditions(pointVectorField &) const
Special correctBoundaryConditions which evaluates fixedValue.
const vector & planeNormal() const
Return plane normal.
A class for managing temporary objects.
static constexpr const zero Zero
Global zero (0)
void correct()
Take over existing mesh position.
void updateMesh()
Update for new mesh topology.
const indirectPrimitivePatch & patch() const
Reference to patch.
const labelList & adaptPatchIDs() const
Patch labels that are being adapted.
virtual const tensorField & forwardT() const
Return face transformation tensor.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
static const pointConstraints & New(const pointMesh &mesh, Args &&... args)
Get existing or create a new MeshObject.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
Ostream & endl(Ostream &os)
Add newline and flush stream.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
prefixOSstream Pout
An Ostream wrapper for parallel output to std::cout.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Mesh consisting of general polyhedral cells.
#define forAll(list, i)
Loop across all elements in list.
const objectRegistry & db() const
Return the local objectRegistry.
virtual bool parallel() const
Are the cyclic planes parallel.
cellMask correctBoundaryConditions()
Class applies a two-dimensional correction to mesh motion point field.
virtual bool separated() const
Are the planes separated.
~motionSmootherAlgo()
Destructor.
messageStream Info
Information stream (uses stdout - output is on the master only)
virtual const fileName & name() const
Return the name of the stream.
virtual void sync(const polyMesh &mesh)
Sync set across coupled patches. Adds coupled points to set.
const pointMesh & mesh() const
Return the mesh reference.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Mesh representing a set of points created from polyMesh.
errorManip< error > abort(error &err)
const dictionary & paramDict() const
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
errorManipArg< error, int > exit(error &err, const int errNo=1)
Output to file stream, using an OSstream.
const TargetType * isA(const Type &t)
Check if dynamic_cast to TargetType is possible.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
const polyMesh & mesh() const
Reference to mesh.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
const dimensionedScalar e
Elementary charge.
label nPoints() const
Number of mesh points.
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
const polyBoundaryMesh & patches
bool required() const
Is 2D correction required, i.e. is the mesh a wedge or slab.
void movePoints()
Update for new mesh geometry.
const pointMesh & pMesh() const
Reference to pointMesh.
static const Vector< scalar > zero
void constrain(GeometricField< Type, pointPatchField, pointMesh > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
const labelList & normalEdgeIndices() const
Return indices of normal edges.
void correctPoints(pointField &p) const
Correct motion points.
Type gMin(const FieldField< Field, Type > &f)
const globalMeshData & globalData() const
Return parallel info.
vector point
Point is a vector.
bool found(const Key &key) const
Return true if hashed entry is found in table.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys and label hasher.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
void setSize(const label newSize)
Alias for resize(const label)
fileName path() const
The complete path.
defineTypeNameAndDebug(combustionModel, 0)
const word & name() const
The patch name.
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
labelList pointLabels(nPoints, -1)
Type gMax(const FieldField< Field, Type > &f)
A list of faces which address into the list of points.