46 void Foam::motionSmootherAlgo::testSyncPositions
59 point(great,great,great)
64 if (
mag(syncedFld[i] -
fld[i]) > maxMag)
67 <<
"On point " << i <<
" point:" <<
fld[i]
68 <<
" synchronised point:" << syncedFld[i]
79 const scalar val =
fld[pointi];
81 if ((val > -great) && (val < great))
86 <<
"Problem : point:" << pointi <<
" value:" << val
102 const face&
f = mesh_.faces()[iter.key()];
106 usedPoints.insert(
f[fp]);
119 const edgeList& edges = mesh_.edges();
121 tmp<scalarField> twght(
new scalarField(edges.size()));
126 wght[edgeI] = 1.0/(edges[edgeI].mag(
points)+small);
133 void Foam::motionSmootherAlgo::minSmooth
136 const PackedBoolList& isAffectedPoint,
142 tmp<pointScalarField> tavgFld = avg
151 label pointi = meshPoints[i];
152 if (isAffectedPoint.get(pointi) == 1)
157 0.5*
fld[pointi] + 0.5*avgFld[pointi]
168 void Foam::motionSmootherAlgo::minSmooth
171 const PackedBoolList& isAffectedPoint,
176 tmp<pointScalarField> tavgFld = avg
185 if (isAffectedPoint.get(pointi) == 1 && isInternalPoint(pointi))
190 0.5*
fld[pointi] + 0.5*avgFld[pointi]
202 void Foam::motionSmootherAlgo::scaleField
211 if (isInternalPoint(iter.key()))
213 fld[iter.key()] *= scale;
223 void Foam::motionSmootherAlgo::scaleField
233 label pointi = meshPoints[i];
237 fld[pointi] *= scale;
244 void Foam::motionSmootherAlgo::subtractField
253 if (isInternalPoint(iter.key()))
255 fld[iter.key()] =
max(0.0,
fld[iter.key()]-
f);
265 void Foam::motionSmootherAlgo::subtractField
275 label pointi = meshPoints[i];
285 bool Foam::motionSmootherAlgo::isInternalPoint(
const label pointi)
const
287 return isInternalPoint_.get(pointi) == 1;
291 void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
293 const label nPointIter,
294 const faceSet& wrongFaces,
297 PackedBoolList& isAffectedPoint
300 isAffectedPoint.setSize(mesh_.nPoints());
303 faceSet nbrFaces(mesh_,
"checkFaces", wrongFaces);
310 for (
label i = 0; i < nPointIter; i++)
312 pointSet nbrPoints(mesh_,
"grownPoints", getPoints(nbrFaces.toc()));
316 const labelList& pCells = mesh_.pointCells(iter.key());
320 const cell& cFaces = mesh_.cells()[pCells[pCelli]];
324 nbrFaces.insert(cFaces[cFacei]);
328 nbrFaces.sync(mesh_);
330 if (i == nPointIter - 2)
334 const face&
f = mesh_.faces()[iter.key()];
337 isAffectedPoint.set(
f[fp], 1);
343 affectedFaces = nbrFaces.toc();
364 displacement_(displacement),
366 oldPoints_(oldPoints),
367 adaptPatchIDs_(adaptPatchIDs),
368 paramDict_(paramDict),
369 isInternalPoint_(mesh_.
nPoints(), 1)
403 return adaptPatchIDs_;
415 oldPoints_ = mesh_.points();
441 displacementBf[
patchi].patchInternalField();
453 forAll(patchSchedule, patchEvalI)
459 if (patchSchedule[patchEvalI].init)
483 displacementBf[
patchi].patchInternalField();
490 setDisplacementPatchFields(adaptPatchIDs_, displacement_);
519 forAll(ppMeshPoints, patchPointi)
521 displacement[ppMeshPoints[patchPointi]] = patchDisp[patchPointi];
537 setDisplacementPatchFields(patchIDs, displacement);
544 forAll(ppMeshPoints, patchPointi)
546 const vector& newDisp = displacement[ppMeshPoints[patchPointi]];
548 if (
mag(newDisp-patchDisp[patchPointi]) > small)
550 const point& pt = mesh.
points()[ppMeshPoints[patchPointi]];
559 Pout<<
"Written " << nVerts <<
" points that are changed to file "
564 forAll(ppMeshPoints, patchPointi)
566 patchDisp[patchPointi] = displacement[ppMeshPoints[patchPointi]];
573 setDisplacement(adaptPatchIDs_, pp_, patchDisp, displacement_);
585 const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
591 forAll(patchSchedule, patchEvalI)
597 if (patchSchedule[patchEvalI].init)
612 forAll(patchSchedule, patchEvalI)
618 if (patchSchedule[patchEvalI].init)
653 Info<<
"Correcting 2-D mesh motion";
655 if (mesh_.globalData().parallel())
658 <<
"2D mesh-motion probably not correct in parallel" <<
endl;
663 const edgeList& edges = mesh_.edges();
670 const edge&
e = edges[neIndices[i]];
672 point& pStart = newPoints[
e.start()];
674 pStart += pn*(pn & (oldPoints[
e.start()] - pStart));
676 point& pEnd = newPoints[
e.end()];
678 pEnd += pn*(pn & (oldPoints[
e.end()] - pEnd));
688 Pout<<
"motionSmootherAlgo::modifyMotionPoints :"
689 <<
" testing sync of newPoints."
691 testSyncPositions(newPoints, 1
e-6*mesh_.bounds().mag());
700 mesh_.clearTetBasePtIs();
707 const scalar errorReduction
710 scalar oldErrorReduction = paramDict_.lookup<scalar>(
"errorReduction");
712 paramDict_.remove(
"errorReduction");
713 paramDict_.add(
"errorReduction", errorReduction);
715 return oldErrorReduction;
722 const bool smoothMesh,
723 const label nAllowableErrors
741 const bool smoothMesh,
742 const label nAllowableErrors
775 displacement_.boundaryField();
782 actualPatchFieldTypes[
patchi] =
803 scale_*displacement_,
804 actualPatchFieldTypes,
811 Pout<<
"scaleMesh : testing sync of totalDisplacement" <<
endl;
817 1
e-6*mesh_.bounds().mag()
824 modifyMotionPoints(tnewPoints.
ref());
836 const bool smoothMesh,
837 const label nAllowableErrors
840 if (!smoothMesh && adaptPatchIDs_.empty())
843 <<
"You specified both no movement on the internal mesh points"
844 <<
" (smoothMesh = false)" <<
nl
845 <<
"and no movement on the patch (adaptPatchIDs is empty)" <<
nl
846 <<
"Hence nothing to adapt."
855 Pout<<
"Entering scaleMesh : coupled patches:" <<
endl;
870 const scalar errorReduction =
871 paramDict.
lookup<scalar>(
"errorReduction");
873 const label nSmoothScale =
887 Info<<
"Moving mesh using displacement scaling :"
888 <<
" min:" <<
gMin(scale_.primitiveField())
889 <<
" max:" <<
gMax(scale_.primitiveField())
896 mesh_.movePoints(newPoints);
901 faceSet wrongFaces(mesh_,
"wrongFaces", mesh_.nFaces()/100+100);
902 checkMesh(
false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
911 wrongFaces.
sync(mesh_);
917 if (
mag(errorReduction) < small)
922 label own = mesh_.faceOwner()[iter.key()];
923 const cell& ownFaces = mesh_.cells()[own];
927 newWrongFaces.
insert(ownFaces[cfI]);
930 if (iter.key() < mesh_.nInternalFaces())
932 label nei = mesh_.faceNeighbour()[iter.key()];
933 const cell& neiFaces = mesh_.cells()[nei];
937 newWrongFaces.
insert(neiFaces[cfI]);
942 wrongFaces.
sync(mesh_);
949 pointSet usedPoints(mesh_,
"usedPoints", getPoints(wrongFaces));
950 usedPoints.
sync(mesh_);
959 getAffectedFacesAndPoints
969 Pout<<
"Faces in error:" << wrongFaces.
size()
970 <<
" with points:" << usedPoints.
size()
974 if (adaptPatchIDs_.size())
977 scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
983 scaleField(usedPoints, errorReduction, scale_);
989 for (
label i = 0; i < nSmoothScale; i++)
991 if (adaptPatchIDs_.size())
1009 minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
1025 Pout<<
"scale_ after smoothing :"
1041 forAll(adaptPatchIDs_, i)
1047 !isA<fixedValuePointPatchVectorField>
1049 displacement_.boundaryField()[
patchi]
1055 <<
" has wrong boundary condition "
1056 << displacement_.boundaryField()[
patchi].type()
1057 <<
" on field " << displacement_.name() <<
nl
1058 <<
"Only type allowed is "
1059 << fixedValuePointPatchVectorField::typeName
1068 const labelList& meshPoints = pp_.meshPoints();
1072 isInternalPoint_.unset(meshPoints[i]);
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
static pointConstraints & New(const pointMesh &mesh)
const Mesh & mesh() const
Return mesh.
Generic GeometricBoundaryField class.
void evaluate()
Evaluate boundary conditions.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
bool insert(const Key &key)
Insert a new entry.
void transfer(HashTable< T, Key, Hash > &)
Transfer the contents of the argument table into this table.
label size() const
Return number of elements in table.
bool found(const Key &) const
Return true if hashedEntry is found in table.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
const objectRegistry & db() const
Return the local objectRegistry.
void setSize(const label)
Reset size of List.
const fileName & name() const
Return the name of the stream.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
label size() const
Return the number of elements in the UPtrList.
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 const transformer & transform() const =0
Return transformation between the coupled patches.
A list of keyword definitions, which are a keyword followed by any number of values (e....
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
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.
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.
motionSmootherAlgo(polyMesh &, pointMesh &, indirectPrimitivePatch &pp, pointVectorField &displacement, pointScalarField &scale, pointField &oldPoints, const labelList &adaptPatchIDs, const dictionary ¶mDict)
Construct from mesh, patches to work on and smoothing parameters.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
void topoChange()
Update for new mesh topology.
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.
bool scaleMesh(const labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
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
void correctBoundaryConditions(pointVectorField &) const
Special correctBoundaryConditions which evaluates fixedValue.
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
fileName path(const word &instance, const fileName &local="") const
Return complete path with alternative instance and local.
const word & name() const
Return name.
const pointMesh & mesh() const
Return the mesh reference.
void constrain(PointField< Type > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
void constrainCorners(PointField< Type > &pf) const
Apply patch-patch constraints only.
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.
Mesh consisting of general polyhedral cells.
const globalMeshData & globalData() const
Return parallel info.
virtual const pointField & points() const
Return raw points.
A class for managing temporary objects.
T & ref() const
Return non-const reference or generate a fatal error.
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 FatalErrorInFunction
Report an error message using Foam::FatalError.
gmvFile<< "tracers "<< particles.size()<< nl;{ pointField positions(particles.size());label particlei=0;forAllConstIter(Cloud< passiveParticle >, particles, iter) { positions[particlei++]=iter().position(mesh);} for(i=0;i< pTraits< point >::nComponents;i++) { forAll(positions, particlei) { gmvFile<< component(positions[particlei], i)<< ' ';} gmvFile<< nl;}}forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const fvPatchList & patches
U correctBoundaryConditions()
#define WarningInFunction
Report a warning using Foam::Warning.
errorManipArg< error, int > exit(error &err, const int errNo=1)
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
errorManip< error > abort(error &err)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
PointField< scalar > pointScalarField
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
prefixOSstream Pout(cout, "Pout")
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Type gMin(const FieldField< Field, Type > &f)
Type gMax(const FieldField< Field, Type > &f)
labelList pointLabels(nPoints, -1)