60 { outputFormatType::PLAIN,
"plain" },
61 { outputFormatType::DICTIONARY,
"dictionary" },
71 { scalingType::LENGTH,
"length" },
72 { scalingType::FORCE,
"force" },
73 { scalingType::MOMENT,
"moment" },
85static inline Ostream& putPlain(Ostream&
os,
const vector& v)
87 return os << v.x() <<
' ' << v.y() <<
' ' << v.z();
95static void writeList(Ostream&
os,
const string& header,
const UList<T>& list)
97 const label len = list.size();
100 os << header.c_str() <<
nl;
106 for (label i=0; i < len; ++i)
133 inputName_(
"positions.in"),
134 outputName_(
"forces.out"),
135 logName_(
"movement.log"),
136 inputFormat_(lumpedPointState::inputFormatType::DICTIONARY),
137 outputFormat_(outputFormatType::DICTIONARY),
147 const dictionary&
dict,
161 inputName_(
"positions.in"),
162 outputName_(
"forces.out"),
163 logName_(
"movement.log"),
177 return (
timeIndex >= lastTrigger_ + calcFrequency_);
193 auto&
points0 = tpoints0.ref();
198 originalIds_.
clear();
199 controllers_.clear();
200 patchControls_.clear();
204 Map<label> pointIdMap;
211 <<
"Incorrect number of pointLabels. Had "
212 << originalIds_.size() <<
" for " <<
points0.
size() <<
" points"
221 for (
const label
id : originalIds_)
223 if (!pointIdMap.insert(
id, pointi))
231 if (!duplicates.empty())
234 <<
"Found duplicate point ids "
240 const dictionary* dictptr =
dict.
findDict(
"controllers");
243 controllers_ = HashPtrTable<lumpedPointController>(*dictptr);
248 (*iter)->remapPointLabels(
points0.
size(), pointIdMap);
256 controllers_.clear();
264 const dictionary& commDict =
dict.
subDict(
"communication");
265 coupler_.readDict(commDict);
267 calcFrequency_ = commDict.
getOrDefault<label>(
"calcFrequency", 1);
271 commDict.readEntry(
"inputName", inputName_);
272 commDict.readEntry(
"outputName", outputName_);
273 commDict.readIfPresent(
"logName", logName_);
276 outputFormat_ = formatNames.get(
"outputFormat", commDict);
281 const dictionary* scaleDict =
nullptr;
283 if ((scaleDict = commDict.findDict(
"scaleInput")) !=
nullptr)
285 for (
int i=0; i < scaleInput_.size(); ++i)
287 const word&
key = scalingNames[scalingType(i)];
291 scaleDict->readIfPresent(key, scaleInput_[i])
292 && scaleInput_[i] > 0
295 Info<<
"Using input " <<
key <<
" multiplier: "
296 << scaleInput_[i] <<
nl;
301 if ((scaleDict = commDict.findDict(
"scaleOutput")) !=
nullptr)
303 for (
int i=0; i < scaleOutput_.size(); ++i)
305 const word&
key = scalingNames[scalingType(i)];
309 scaleDict->readIfPresent(key, scaleOutput_[i])
310 && scaleOutput_[i] > 0
313 Info<<
"Using output " <<
key <<
" multiplier: "
314 << scaleOutput_[i] <<
nl;
319 state0_ = lumpedPointState
331 state0_.scalePoints(scaleInput_[scalingType::LENGTH]);
342 return hasPatchControl(pp.index());
348 const pointPatch& fpatch
351 return hasInterpolator(fpatch.index());
360 const auto ctrlIter = patchControls_.cfind(pp.index());
362 if (!ctrlIter.good())
365 <<
"No controllers for patch " << pp.name()
369 const patchControl& ctrl = *ctrlIter;
371 for (
const word& ctrlName : ctrl.names_)
373 const auto iter = controllers_.cfind(ctrlName);
378 <<
"No controller: " << ctrlName <<
nl
379 <<
" For patch " << pp.name()
394 const pointField& lumpedCentres0 = state0().points();
396 const label patchIndex = pp.index();
401 patchControl& ctrl = patchControls_(patchIndex);
402 ctrl.names_ = ctrlNames;
404 labelList& faceToPoint = ctrl.faceToPoint_;
405 faceToPoint.
resize(pp.size(), -1);
407 checkPatchControl(pp);
409 const faceList& faces = pp.boundaryMesh().mesh().faces();
414 for (
const word& ctrlName : ctrl.names_)
416 const auto iter = controllers_.
cfind(ctrlName);
421 <<
"No controller: " << ctrlName <<
nl
430 if (ctrl.names_.size() && subsetPointIds.empty())
433 <<
"Controllers specified, but without any points" <<
nl
438 treeDataPoint treePoints
441 subsetPointIds.sortedToc(),
442 subsetPointIds.size()
446 treeBoundBox bb(lumpedCentres0);
449 indexedOctree<treeDataPoint> ppTree
458 const scalar searchDistSqr(
sqr(GREAT));
460 const label patchStart = pp.start();
463 const point fc(faces[patchStart + patchFacei].centre(
points0));
466 faceToPoint[patchFacei] =
467 treePoints.pointLabel
469 ppTree.findNearest(fc, searchDistSqr).index()
475 Pout<<
"Added face mapping for patch: " << patchIndex <<
endl;
482 const pointPatch& fpatch,
487 const pointField& lumpedCentres0 = state0().points();
489 const label patchIndex = fpatch.index();
491 patchControl& ctrl = patchControls_(patchIndex);
493 List<lumpedPointInterpolator>& interpList = ctrl.interp_;
497 Map<labelHashSet> adjacency;
502 for (
const word& ctrlName : ctrl.names_)
504 const auto iter = controllers_.
cfind(ctrlName);
509 <<
"No controller: " << ctrlName <<
nl
525 adjacency(prev).insert(curr);
526 adjacency(curr).insert(prev);
528 else if (!adjacency.found(curr))
530 adjacency(curr).
clear();
535 if (ctrl.names_.empty())
538 const label len = state0().size();
540 for (label i=0; i < len; ++i)
542 const label curr = i;
546 const label prev = i-1;
547 adjacency(prev).insert(curr);
548 adjacency(curr).insert(prev);
550 else if (!adjacency.found(curr))
552 adjacency(curr).clear();
557 if (ctrl.names_.size() && subsetPointIds.empty())
560 <<
"Controllers specified, but without any points" <<
nl
566 Map<labelPairList> adjacencyPairs;
572 DynamicList<labelPair> pairs;
576 DynamicList<scalar> acuteAngles;
580 const label nearest = iter.key();
582 labelList neighbours(iter.val().sortedToc());
584 const label len = neighbours.size();
589 const point& nearPt = lumpedCentres0[nearest];
591 for (label j=1; j < len; ++j)
593 for (label i=0; i < j; ++i)
595 labelPair neiPair(neighbours[i], neighbours[j]);
600 lumpedCentres0[neiPair.first()],
601 lumpedCentres0[neiPair.second()]
605 if (tri.pointToBarycentric(tri.a(), bary) > SMALL)
608 pairs.append(neiPair);
616 acuteAngles.append(-(ab & ac));
621 if (pairs.size() > 1)
628 adjacencyPairs.insert(nearest, pairs);
634 Info<<
"Adjacency table for patch: " << fpatch.name() <<
nl;
636 for (
const label own : adjacency.sortedToc())
639 for (
const label nei : adjacency[own].
sortedToc())
644 if (originalIds_.size())
646 Info<<
" # " << originalIds_[own] <<
" =>";
647 for (
const label nei : adjacency[own].
sortedToc())
649 Info<<
' ' << originalIds_[nei];
658 treeDataPoint treePoints
661 subsetPointIds.sortedToc(),
662 subsetPointIds.size()
665 treeBoundBox bb(lumpedCentres0);
668 indexedOctree<treeDataPoint> ppTree
680 const scalar searchDistSqr(
sqr(GREAT));
682 const labelList& meshPoints = fpatch.meshPoints();
684 interpList.
resize(meshPoints.size());
686 DynamicList<scalar> unsortedNeiWeightDist;
687 DynamicList<label> unsortedNeighbours;
689 forAll(meshPoints, pointi)
694 const label nearest =
695 treePoints.pointLabel
697 ppTree.findNearest(ptOnMesh, searchDistSqr).index()
700 interpList[pointi].nearest(nearest);
717 const point& nearPt = lumpedCentres0[nearest];
721 const labelPairList& adjacentPairs = adjacencyPairs[nearest];
723 unsortedNeighbours = adjacency[nearest].toc();
724 unsortedNeiWeightDist.resize(unsortedNeighbours.size());
726 forAll(unsortedNeighbours, nbri)
728 unsortedNeiWeightDist[nbri] =
729 magSqr(ptOnMesh - lumpedCentres0[unsortedNeighbours[nbri]]);
738 for (
const label nbri : distOrder)
740 const label nextPointi = unsortedNeighbours[nbri];
742 const point& nextPt = lumpedCentres0[nextPointi];
746 const scalar weight =
747 (toMeshPt.vec() & toNextPt.unitVec()) / toNextPt.mag();
757 unsortedNeiWeightDist[nbri] = weight;
760 distOrder[ngood] = nbri;
765 distOrder.resize(ngood);
767 if (distOrder.size() < 1)
772 UIndirectList<label> neighbours(unsortedNeighbours, distOrder);
773 UIndirectList<scalar> neiWeight(unsortedNeiWeightDist, distOrder);
775 bool useFirst =
true;
777 if (neighbours.size() > 1 && adjacentPairs.size())
782 neiPointid.set(neighbours);
784 for (
const labelPair& ends : adjacentPairs)
786 label nei1 = ends.first();
787 label nei2 = ends.second();
789 if (!neiPointid.test(nei1) || !neiPointid.test(nei2))
794 else if (neighbours.find(nei2) < neighbours.find(nei1))
798 std::swap(nei1, nei2);
802 triFace triF(nearest, nei1, nei2);
806 triF.tri(lumpedCentres0).pointToBarycentric(ptOnMesh, bary)
812 interpList[pointi].set(triF, bary);
823 interpList[pointi].next(neighbours.first(), neiWeight.first());
831 const polyMesh& pmesh
834 const label nLumpedPoints = state0().
size();
836 List<scalar> zoneAreas(nLumpedPoints,
Zero);
838 if (patchControls_.empty())
841 <<
"Attempted to calculate areas without setMapping()"
846 const polyBoundaryMesh&
patches = pmesh.boundaryMesh();
849 if (isA<fvMesh>(pmesh))
851 const fvMesh&
mesh = dynamicCast<const fvMesh>(pmesh);
859 const label patchIndex = iter.key();
860 const patchControl& ctrl = iter.val();
862 const labelList& faceToPoint = ctrl.faceToPoint_;
864 const polyPatch& pp =
patches[patchIndex];
869 const label pointIndex = faceToPoint[patchFacei];
878 zoneAreas[pointIndex] +=
mag(patchSf[patchIndex][patchFacei]);
891 const polyMesh& pmesh,
892 List<vector>& forces,
893 List<vector>& moments
896 const label nLumpedPoints = state0().size();
898 forces.resize(nLumpedPoints);
899 moments.resize(nLumpedPoints);
901 if (patchControls_.empty())
904 <<
"Attempted to calculate forces without setMapping()"
907 forces.resize(nLumpedPoints,
Zero);
908 moments.resize(nLumpedPoints,
Zero);
917 const pointField& lumpedCentres = state().points();
919 const polyBoundaryMesh&
patches = pmesh.boundaryMesh();
921 const word pName(forcesDict_.getOrDefault<word>(
"p",
"p"));
922 scalar pRef = forcesDict_.getOrDefault<scalar>(
"pRef", 0);
923 scalar rhoRef = forcesDict_.getOrDefault<scalar>(
"rhoRef", 1);
927 PtrMap<vectorField> forceOnPatches;
932 if (isA<fvMesh>(pmesh) && pPtr)
934 const fvMesh&
mesh = dynamicCast<const fvMesh>(pmesh);
954 const label patchIndex = iter.key();
955 const patchControl& ctrl = iter.val();
957 const labelList& faceToPoint = ctrl.faceToPoint_;
959 if (!forceOnPatches.found(patchIndex))
969 * patchSf[patchIndex] * (patchPress[patchIndex] - pRef)
974 const vectorField& forceOnPatch = *forceOnPatches[patchIndex];
976 const polyPatch& pp =
patches[patchIndex];
981 const label pointIndex = faceToPoint[patchFacei];
990 forces[pointIndex] += forceOnPatch[patchFacei];
997 patchCf[patchIndex][patchFacei]
998 - lumpedCentres[pointIndex]
1002 moments[pointIndex] += lever ^ forceOnPatch[patchFacei];
1021 const pointPatch& fpatch,
1025 return pointsDisplacement(state(), fpatch,
points0);
1032 const lumpedPointState& state,
1033 const pointPatch& fpatch,
1037 const label patchIndex = fpatch.index();
1040 const pointField& lumpedCentres0 = state0().points();
1043 const pointField& lumpedCentres = state.points();
1046 const tensorField& localToGlobal = state.rotations();
1048 const labelList& meshPoints = fpatch.meshPoints();
1053 auto& disp = tdisp.ref();
1056 const List<lumpedPointInterpolator>& interpList
1057 = patchControls_[patchIndex].interp_;
1059 forAll(meshPoints, pointi)
1061 const lumpedPointInterpolator& interp = interpList[pointi];
1065 const vector origin0 = interp.interpolate(lumpedCentres0);
1066 const vector origin = interp.interpolate(lumpedCentres);
1067 const tensor rotTensor = interp.interpolate(localToGlobal);
1069 disp[pointi] = (rotTensor & (
p0 - origin0)) + origin -
p0;
1079 const lumpedPointState& state,
1080 const pointPatch& fpatch,
1084 const label patchIndex = fpatch.index();
1087 const pointField& lumpedCentres0 = state0().points();
1090 const pointField& lumpedCentres = state.points();
1093 const tensorField& localToGlobal = state.rotations();
1095 const labelList& meshPoints = fpatch.meshPoints();
1100 auto& disp = tdisp.ref();
1103 const List<lumpedPointInterpolator>& interpList =
1104 patchControls_[patchIndex].interp_;
1106 forAll(meshPoints, pointi)
1108 const lumpedPointInterpolator& interp = interpList[pointi];
1112 const vector origin0 = interp.interpolate(lumpedCentres0);
1113 const vector origin = interp.interpolate(lumpedCentres);
1114 const tensor rotTensor = interp.interpolate(localToGlobal);
1116 disp[pointi] = (rotTensor & (
p0 - origin0)) + origin;
1132 lumpedPointState prev = state_;
1134 const bool status = state_.readData
1137 coupler().resolveFile(inputName_),
1138 state0().rotationOrder(),
1142 scalePoints(state_);
1144 state_.relax(relax_, prev);
1153 const UList<vector>& forces,
1154 const UList<vector>& moments,
1155 const outputFormatType fmt,
1156 const Tuple2<scalar, scalar>* timesWritten
1159 const bool writeMoments = (moments.size() == forces.size());
1161 if (fmt == outputFormatType::PLAIN)
1163 os <<
"########" <<
nl;
1166 os <<
"# Time value=" << timesWritten->first() <<
nl
1167 <<
"# Time prev=" << timesWritten->second() <<
nl;
1169 os <<
"# size=" << this->size() <<
nl
1170 <<
"# columns (points) (forces)";
1179 bool report =
false;
1180 scalar scaleLength = scaleOutput_[scalingType::LENGTH];
1181 scalar scaleForce = scaleOutput_[scalingType::FORCE];
1182 scalar scaleMoment = scaleOutput_[scalingType::MOMENT];
1184 if (scaleLength > 0)
1204 if (scaleMoment > 0)
1216 os <<
"# scaling points=" << scaleLength
1217 <<
" forces=" << scaleForce;
1221 os <<
" moments=" << scaleMoment;
1227 os <<
"########" <<
nl;
1231 const point& pt = state0().points()[i];
1233 putPlain(
os, scaleLength * pt) <<
' ';
1235 if (i < forces.size())
1237 const vector val(scaleForce * forces[i]);
1242 putPlain(
os, vector::zero);
1248 if (i < moments.size())
1250 const vector val(scaleMoment * moments[i]);
1255 putPlain(
os, vector::zero);
1268 os <<
"////////" <<
nl;
1291 const UList<vector>& forces,
1292 const UList<vector>& moments,
1293 const Tuple2<scalar, scalar>* timesWritten
1305 coupler().resolveFile(outputName_)
1308 writeData(
os, forces, moments, outputFormat_, timesWritten);
1315 coupler().resolveFile(logName_),
1320 writeData(
os, forces, moments, outputFormatType::PLAIN, timesWritten);
Inter-processor communication reduction functions.
const dimensionSet & dimensions() const
Return dimensions.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
EnumType get(const word &enumName) const
The enumeration corresponding to the given name.
GeometricBoundaryField< vector, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
bool set(const Key &key)
Same as insert (no value to overwrite)
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
void resize(const label len)
Adjust allocated size of list.
void clear()
Clear the list, i.e. set size to zero.
bool inplaceClip(T &val) const
Inplace clip value by the min/max limits.
static MinMax< scalar > zero_one()
A 0-1 range corresponding to the pTraits zero, one.
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
dictionary subOrEmptyDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX, const bool mandatory=false) const
dictionary * findDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX)
Find and return a sub-dictionary pointer if present.
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Foam::dictionary writeDict() const
Write to dictionary.
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
const surfaceVectorField & Sf() const
Return cell face area vectors.
bool hasPatchControl(const label patchIndex) const
Check if patch control exists for specified patch.
static const word canonicalName
The canonical name ("lumpedPointMovement") for the dictionary.
tmp< pointField > pointsDisplacement(const pointPatch &fpatch, const pointField &points0) const
Displace points according to the current state.
static const Enum< outputFormatType > formatNames
Names for the output format types.
tmp< pointField > pointsPosition(const lumpedPointState &state, const pointPatch &fpatch, const pointField &points0) const
The points absolute position according to specified state.
void readDict(const dictionary &dict)
Update settings from dictionary.
void setInterpolator(const pointPatch &fpatch, const pointField &points0)
Check if patch control exists for specified patch.
void setPatchControl(const polyPatch &pp, const wordList &ctrlNames, const pointField &points0)
Define pressure-zones mapping for faces in the specified patches.
bool readState()
Read state from file, applying relaxation as requested.
void couplingCompleted(const label timeIndex) const
Register that coupling is completed at this calcFrequency.
bool couplingPending(const label timeIndex) const
Check if coupling is pending (according to the calcFrequency)
static const Enum< scalingType > scalingNames
Names for the scaling types.
void checkPatchControl(const polyPatch &pp) const
Check if patch control exists for specified patch.
outputFormatType
Output format types.
bool writeData(Ostream &os, const UList< vector > &forces, const UList< vector > &moments, const outputFormatType fmt=outputFormatType::PLAIN, const Tuple2< scalar, scalar > *timesWritten=nullptr) const
Write points, forces, moments. Only call from the master process.
static int debug
Debug switch.
bool hasInterpolator(const pointPatch &fpatch) const
Check if patch control exists for specified patch.
lumpedPointMovement()
Default construct.
bool forcesAndMoments(const polyMesh &pmesh, List< vector > &forces, List< vector > &moments) const
The forces and moments acting on each pressure-zone.
scalingType
Output format types.
List< scalar > areas(const polyMesh &pmesh) const
The areas for each pressure-zone.
static const Enum< inputFormatType > formatNames
Names for the input format types.
static const Enum< eulerOrder > eulerOrderNames
splitCell * master() const
Tensor of scalars, i.e. Tensor<scalar>.
A class for managing temporary objects.
@ END_STATEMENT
End entry [isseparator].
@ BEGIN_LIST
Begin list [isseparator].
@ END_LIST
End list [isseparator].
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
const volScalarField & p0
const polyBoundaryMesh & patches
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
List< label > sortedToc(const UList< bool > &bools)
Return the (sorted) values corresponding to 'true' entries.
int debugSwitch(const char *name, const int deflt=0)
Lookup debug switch or add default value.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
void writeList(vtk::formatter &fmt, const UList< uint8_t > &values)
Write a list of uint8_t values.
const dimensionSet dimPressure
Pair< label > labelPair
A pair of labels.
List< word > wordList
A List of words.
List< labelPair > labelPairList
List of labelPairs.
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
line< point, const point & > linePointRef
A line using referred points.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Barycentric2D< scalar > barycentric2D
A scalar version of the templated Barycentric2D.
vector point
Point is a vector.
Ostream & endl(Ostream &os)
Add newline and flush stream.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
triangle< point, const point & > triPointRef
A triangle using referred points.
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static void writeData(Ostream &os, const Type &val)
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
labelList pointLabels(nPoints, -1)
#define forAll(list, i)
Loop across all elements in list.
#define forAllIters(container, iter)
Iterate across all elements in the container object.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
pointField points0(pointIOField(IOobject("points", mesh.time().constant(), polyMesh::meshSubDir, mesh, IOobject::MUST_READ, IOobject::NO_WRITE, false)))