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
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();
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]
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);
533 OFstream str(
mesh.
db().
path()/
"changedPoints.obj");
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 " 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));
667 point& pEnd = newPoints[
e.end()];
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();
767 actualPatchFieldTypes.
setSize(pfld.size());
770 if (
isA<fixedValuePointPatchField<vector>>(pfld[patchi]))
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()
812 tmp<pointField> tnewPoints(oldPoints_ + totalDisplacement.primitiveField());
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." 844 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
846 Pout<<
"Entering scaleMesh : coupled patches:" <<
endl;
851 const coupledPolyPatch&
pp =
852 refCast<const coupledPolyPatch>(
patches[patchi]);
854 Pout<<
'\t' << patchi <<
'\t' <<
pp.name()
855 <<
" parallel:" <<
pp.parallel()
856 <<
" separated:" <<
pp.separated()
857 <<
" forwardT:" <<
pp.forwardT().size()
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);
908 if (
returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
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 :" 1034 const pointBoundaryMesh&
patches = pMesh_.boundary();
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());
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
const polyBoundaryMesh & pbm
const labelList patchIDs(pbm.patchSet(polyPatchNames, false, true).sortedToc())
"blocking" : (MPI_Bsend, MPI_Recv)
const polyMesh & mesh() const
Reference to mesh.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
labelList pointLabels(nPoints, -1)
label nPoints() const noexcept
Number of mesh points.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Type gMin(const FieldField< Field, Type > &f)
List< edge > edgeList
List of edge.
void correctBoundaryConditions(pointVectorField &) const
Special correctBoundaryConditions which evaluates fixedValue.
const word & name() const noexcept
Return the object name.
static const pointConstraints & New(const pointMesh &mesh, Args &&... args)
Get existing or create a new MeshObject. Registered with typeName.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
constexpr char nl
The newline '\n' character (0x0a)
const pointMesh & pMesh() const
Reference to pointMesh.
const dictionary & paramDict() const
const labelList & normalEdgeIndices() const
Return indices of normal edges.
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
cellMask correctBoundaryConditions()
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool found(const T &val, label pos=0) const
Same as contains()
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
List< lduScheduleEntry > lduSchedule
A List of lduSchedule entries.
GeometricBoundaryField< vector, pointPatchField, pointMesh > Boundary
Type of boundary fields.
Ignore writing from objectRegistry::writeObject()
void movePoints()
Update for new mesh geometry.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
labelList faceLabels(nFaceLabels)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
fileName path() const
The complete path.
virtual const pointField & points() const
Return raw points.
Mesh representing a set of points created from polyMesh.
#define forAll(list, i)
Loop across all elements in list.
void correctPoints(pointField &p) const
Correct motion points.
void correct()
Take over existing mesh position.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
static const char *const typeName
Typename for Field.
A list of faces which address into the list of points.
void evaluate()
Evaluate boundary conditions.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
"scheduled" : (MPI_Send, MPI_Recv)
bool required() const
Is 2D correction required, i.e. is the mesh a wedge or slab.
const polyMesh & mesh() const noexcept
Return the mesh reference.
void setDisplacementPatchFields()
Set patch fields on displacement to be consistent with.
Class applies a two-dimensional correction to mesh motion point field.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label size() const noexcept
The number of entries in the list.
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
const objectRegistry & db() const noexcept
Return the local objectRegistry.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
const globalMeshData & globalData() const
Return parallel info.
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
void modifyMotionPoints(pointField &newPoints) const
Apply optional point constraint (2d correction)
const Mesh & mesh() const noexcept
Return mesh.
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
const vector & planeNormal() const
Return plane normal.
defineTypeNameAndDebug(combustionModel, 0)
const indirectPrimitivePatch & patch() const
Reference to patch.
void setSize(const label newLen)
Same as resize()
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))
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
List< word > wordList
List of word.
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
tmp< pointField > curPoints() const
Get the current points (oldPoints+scale*displacement)
const polyBoundaryMesh & patches
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
messageStream Info
Information stream (stdout output on master, null elsewhere)
void updateMesh()
Update for new mesh topology.
Mesh consisting of general polyhedral cells.
const Type * isA(const U &obj)
Check if dynamic_cast to Type is possible.
const labelList & adaptPatchIDs() const
Patch labels that are being adapted.
List< label > labelList
A List of labels.
A class for managing temporary objects.
void constrain(GeometricField< Type, pointPatchField, pointMesh > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Do not request registration (bool: false)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
~motionSmootherAlgo()
Destructor.
static constexpr const zero Zero
Global zero (0)