49void Foam::motionSmootherAlgo::testSyncPositions
51 const pointField&
fld,
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)
105 usedPoints.insert(mesh_.faces()[
faceId]);
117 const edgeList& edges = mesh_.edges();
120 auto& wght = twght.ref();
124 wght[edgei] = 1.0/(edges[edgei].mag(
points)+SMALL);
131void 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]
165void 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]
199void Foam::motionSmootherAlgo::scaleField
208 if (isInternalPoint_.test(pointi))
210 fld[pointi] *= scale;
220void Foam::motionSmootherAlgo::scaleField
228 for (
const label pointi : meshPoints)
232 fld[pointi] *= scale;
239void Foam::motionSmootherAlgo::subtractField
248 if (isInternalPoint_.test(pointi))
260void Foam::motionSmootherAlgo::subtractField
268 for (
const label pointi : meshPoints)
278void 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();
341 displacement_(displacement),
343 oldPoints_(oldPoints),
344 adaptPatchIDs_(adaptPatchIDs),
345 paramDict_(paramDict),
347 isInternalPoint_(mesh_.
nPoints(), true)
381 return adaptPatchIDs_;
393 oldPoints_ = mesh_.points();
414 for (
const label patchi : patchIDs)
416 displacementBf[patchi] ==
417 displacementBf[patchi].patchInternalField();
427 displacement.
mesh().globalData().patchSchedule();
429 for (
const auto& schedEval : patchSchedule)
431 const label patchi = schedEval.patch;
433 if (!adaptPatchSet.
found(patchi))
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();
581 for (
const auto& schedEval : patchSchedule)
583 const label patchi = schedEval.patch;
585 if (adaptPatchSet.
found(patchi))
589 displacementBf[patchi]
594 displacementBf[patchi]
602 for (
const auto& schedEval : patchSchedule)
604 const label patchi = schedEval.patch;
606 if (!adaptPatchSet.
found(patchi))
610 displacementBf[patchi]
615 displacementBf[patchi]
643 Info<<
"Correcting 2-D mesh motion";
645 if (mesh_.globalData().parallel())
648 <<
"2D mesh-motion probably not correct in parallel" <<
endl;
653 const edgeList& edges = mesh_.edges();
658 for (
const label edgei : neIndices)
661 const edge&
e = edges[edgei];
663 point& pStart = newPoints[
e.start()];
665 pStart += pn*(pn & (oldPoints[
e.start()] - pStart));
669 pEnd += pn*(pn & (oldPoints[
e.
end()] - pEnd));
679 Pout<<
"motionSmootherAlgo::modifyMotionPoints :"
680 <<
" testing sync of newPoints."
682 testSyncPositions(newPoints, 1
e-6*mesh_.bounds().mag());
691 mesh_.clearTetBasePtIs();
692 pp_.movePoints(mesh_.points());
698 const scalar errorReduction
701 scalar old = paramDict_.get<scalar>(
"errorReduction");
703 paramDict_.remove(
"errorReduction");
704 paramDict_.add(
"errorReduction", errorReduction);
713 const bool smoothMesh,
714 const label nAllowableErrors
732 const bool smoothMesh,
733 const label nAllowableErrors
759 actualPatchTypes[patchi] = pbm[patchi].type();
766 displacement_.boundaryField();
773 actualPatchFieldTypes[patchi] =
778 actualPatchFieldTypes[patchi] = pfld[patchi].type();
788 mesh_.time().timeName(),
794 scale_*displacement_,
795 actualPatchFieldTypes,
802 Pout<<
"scaleMesh : testing sync of totalDisplacement" <<
endl;
808 1
e-6*mesh_.bounds().mag()
815 modifyMotionPoints(tnewPoints.
ref());
827 const bool smoothMesh,
828 const label nAllowableErrors
831 if (!smoothMesh && adaptPatchIDs_.empty())
834 <<
"You specified both no movement on the internal mesh points"
835 <<
" (smoothMesh = false)" <<
nl
836 <<
"and no movement on the patch (adaptPatchIDs is empty)" <<
nl
837 <<
"Hence nothing to adapt."
846 Pout<<
"Entering scaleMesh : coupled patches:" <<
endl;
852 refCast<const coupledPolyPatch>(
patches[patchi]);
854 Pout<<
'\t' << patchi <<
'\t' << pp.
name()
863 const scalar errorReduction = get<scalar>
867 const label nSmoothScale = get<label>
882 Info<<
"Moving mesh using displacement scaling :"
883 <<
" min:" <<
gMin(scale_.primitiveField())
884 <<
" max:" <<
gMax(scale_.primitiveField())
891 mesh_.movePoints(newPoints);
896 faceSet wrongFaces(mesh_,
"wrongFaces", mesh_.nFaces()/100+100);
915 wrongFaces.
sync(mesh_);
921 if (
mag(errorReduction) < SMALL)
924 for (
const label facei : wrongFaces)
926 const label own = mesh_.faceOwner()[facei];
927 const cell& ownFaces = mesh_.cells()[own];
929 newWrongFaces.
insert(ownFaces);
931 if (facei < mesh_.nInternalFaces())
933 const label nei = mesh_.faceNeighbour()[facei];
934 const cell& neiFaces = mesh_.cells()[nei];
936 newWrongFaces.
insert(neiFaces);
939 wrongFaces.transfer(newWrongFaces);
940 wrongFaces.sync(mesh_);
947 pointSet usedPoints(mesh_,
"usedPoints", getPoints(wrongFaces));
948 usedPoints.
sync(mesh_);
954 bitSet isAffectedPoint(mesh_.nPoints());
955 getAffectedFacesAndPoints
965 Pout<<
"Faces in error:" << wrongFaces.
size()
966 <<
" with points:" << usedPoints.
size()
970 if (adaptPatchIDs_.size())
973 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
979 scaleField(usedPoints, errorReduction, scale_);
985 for (label i = 0; i < nSmoothScale; ++i)
987 if (adaptPatchIDs_.size())
1005 minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
1021 Pout<<
"scale_ after smoothing :"
1037 for (
const label patchi : adaptPatchIDs_)
1041 !isA<fixedValuePointPatchVectorField>
1043 displacement_.boundaryField()[patchi]
1049 <<
" has wrong boundary condition "
1050 << displacement_.boundaryField()[patchi].type()
1051 <<
" on field " << displacement_.name() <<
nl
1052 <<
"Only type allowed is "
1053 << fixedValuePointPatchVectorField::typeName
1062 isInternalPoint_.unset(pp_.meshPoints());
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
const Mesh & mesh() const
Return mesh.
void evaluate()
Evaluate boundary conditions.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
bool found(const Key &key) const
Return true if hashed entry is found in table.
label size() const noexcept
The number of elements in table.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
const word & name() const noexcept
Return the object name.
const objectRegistry & db() const noexcept
Return the local objectRegistry.
fileName path() const
The complete path.
void setSize(const label n)
Alias for resize()
Output to file stream, using an OSstream.
virtual const fileName & name() const
Read/write access to the name of the stream.
A list of faces which address into the list of points.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
bool found(const T &val, label pos=0) const
True if the value if found in the list.
iterator end() noexcept
Return an iterator to end traversing the UList.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
A cell is defined as a list of faces with extra functionality.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
virtual bool separated() const
Are the planes separated.
virtual bool parallel() const
Are the cyclic planes parallel.
virtual const tensorField & forwardT() const
Return face transformation tensor.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
A FixedValue boundary condition for pointField.
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
~motionSmootherAlgo()
Destructor.
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
const labelList & adaptPatchIDs() const
Patch labels that are being adapted.
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
void correct()
Take over existing mesh position.
const polyMesh & mesh() const
Reference to mesh.
void movePoints()
Update for new mesh geometry.
const indirectPrimitivePatch & patch() const
Reference to patch.
const pointMesh & pMesh() const
Reference to pointMesh.
const dictionary & paramDict() const
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
void updateMesh()
Update for new mesh topology.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
const word & name() const noexcept
The patch name.
const pointMesh & mesh() const noexcept
Return the mesh reference.
Mesh representing a set of points created from polyMesh.
const pointBoundaryMesh & boundary() const
Return reference to boundary mesh.
virtual void sync(const polyMesh &mesh)
Sync set across coupled patches. Adds coupled points to set.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
const globalMeshData & globalData() const
Return parallel info.
virtual const pointField & points() const
Return raw points.
label nPoints() const noexcept
Number of mesh points.
A class for managing temporary objects.
Class applies a two-dimensional correction to mesh motion point field.
bool required() const
Is 2D correction required, i.e. is the mesh a wedge or slab.
void correctPoints(pointField &p) const
Correct motion points.
const labelList & normalEdgeIndices() const
Return indices of normal edges.
const vector & planeNormal() const
Return plane normal.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
List< label > labelList
A List of labels.
const TargetType * isA(const Type &t)
Check if dynamic_cast to TargetType is possible.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
List< edge > edgeList
A List of edges.
Type gMin(const FieldField< Field, Type > &f)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
labelList pointLabels(nPoints, -1)
cellMask correctBoundaryConditions()
#define forAll(list, i)
Loop across all elements in list.
A non-counting (dummy) refCount.