56int main(
int argc,
char *argv[])
60 "Calculates the inertia tensor and principal axes and moments"
61 " of the specified surface.\n"
62 "Inertia can either be of the solid body or of a thin shell."
65 argList::noParallel();
66 argList::addArgument(
"input",
"The input surface file");
68 argList::addBoolOption
71 "Inertia of a thin shell"
79 "kg/m3 for solid properties, kg/m2 for shell properties"
86 "Inertia relative to this point, not the centre of mass"
105 momentOfInertia::massPropertiesShell(surf, density, m, cM, J);
109 momentOfInertia::massPropertiesSolid(surf, density, m, cM, J);
115 <<
"Negative mass detected, the surface may be inside-out." <<
endl;
126 while ((
magSqr(eVal) < VSMALL) && pertI < 10)
129 <<
"No eigenValues found, shape may have symmetry, "
130 <<
"perturbing inertia tensor diagonal" <<
endl;
132 J.xx() *= 1.0 + SMALL*rand.sample01<scalar>();
133 J.yy() *= 1.0 + SMALL*rand.sample01<scalar>();
134 J.zz() *= 1.0 + SMALL*rand.sample01<scalar>();
143 bool showTransform =
true;
147 (
mag(eVec.x() ^ eVec.y()) > (1.0 - SMALL))
148 && (
mag(eVec.y() ^ eVec.z()) > (1.0 - SMALL))
149 && (
mag(eVec.z() ^ eVec.x()) > (1.0 - SMALL))
157 eVec.z() *
sign((eVec.x() ^ eVec.y()) & eVec.z())
166 cartesian[0] =
vector(1, 0, 0);
167 cartesian[1] =
vector(0, 1, 0);
168 cartesian[2] =
vector(0, 0, 1);
174 principal[0] = eVec.x();
175 principal[1] = eVec.y();
176 principal[2] = eVec.z();
178 scalar maxMagDotProduct = -GREAT;
188 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
190 if (magDotProduct > maxMagDotProduct)
192 maxMagDotProduct = magDotProduct;
203 cartesian[
match.first()] & principal[
match.second()]
213 tPrincipal[
match.second()] *= -1;
215 tPrincipal[(
match.second() + 1) % 3] =
216 principal[(
match.second() + 2) % 3];
218 tPrincipal[(
match.second() + 2) % 3] =
219 principal[(
match.second() + 1) % 3];
221 principal = tPrincipal;
225 tEVal[(
match.second() + 1) % 3] = eVal[(
match.second() + 2) % 3];
227 tEVal[(
match.second() + 2) % 3] = eVal[(
match.second() + 1) % 3];
232 label permutationDelta =
match.second() -
match.first();
234 if (permutationDelta != 0)
238 permutationDelta += 3;
244 for (label i = 0; i < 3; i++)
246 tPrincipal[i] = principal[(i + permutationDelta) % 3];
248 tEVal[i] = eVal[(i + permutationDelta) % 3];
251 principal = tPrincipal;
256 label matchedAlready =
match.first();
260 maxMagDotProduct = -GREAT;
264 if (cI == matchedAlready)
271 if (pI == matchedAlready)
276 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
278 if (magDotProduct > maxMagDotProduct)
280 maxMagDotProduct = magDotProduct;
291 cartesian[
match.first()] & principal[
match.second()]
294 if (sense < 0 || (
match.second() -
match.first()) != 0)
296 principal[
match.second()] *= -1;
300 tPrincipal[(matchedAlready + 1) % 3] =
301 principal[(matchedAlready + 2) % 3]*-sense;
303 tPrincipal[(matchedAlready + 2) % 3] =
304 principal[(matchedAlready + 1) % 3]*-sense;
306 principal = tPrincipal;
310 tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
312 tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
317 eVec =
tensor(principal[0], principal[1], principal[2]);
332 <<
"Non-unique eigenvectors, cannot compute transformation "
333 <<
"from Cartesian axes" <<
endl;
335 showTransform =
false;
340 scalar surfaceArea = 0;
346 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
349 <<
"Illegal triangle " << facei <<
" vertices " <<
f
350 <<
" coords " <<
f.points(surf.points()) <<
endl;
364 <<
"Density: " << density <<
nl
365 <<
"Mass: " << m <<
nl
366 <<
"Centre of mass: " << cM <<
nl
367 <<
"Surface area: " << surfaceArea <<
nl
368 <<
"Inertia tensor around centre of mass: " <<
nl << J <<
nl
369 <<
"eigenValues (principal moments): " << eVal <<
nl
370 <<
"eigenVectors (principal axes): " <<
nl
371 << eVec.x() <<
nl << eVec.y() <<
nl << eVec.z() <<
endl;
375 Info<<
"Transform tensor from reference state (orientation):" <<
nl
377 <<
"Rotation tensor required to transform "
378 "from the body reference frame to the global "
379 "reference frame, i.e.:" <<
nl
380 <<
"globalVector = orientation & bodyLocalVector"
384 <<
"Entries for sixDoFRigidBodyDisplacement boundary condition:"
388 .writeEntry(
"mass", m)
389 .writeEntry(
"centreOfMass", cM)
390 .writeEntry(
"momentOfInertia", eVal)
391 .writeEntry(
"orientation", eVec.T());
396 Info<<
nl <<
"Inertia tensor relative to " << refPt <<
": " <<
nl
397 << momentOfInertia::applyParallelAxisTheorem(m, cM, J, refPt)
403 Info<<
nl <<
"Writing scaled principal axes at centre of mass of "
404 << surfFileName <<
" to " << str.name() <<
endl;
406 scalar scale =
mag(cM - surf.points()[0])/eVal.component(
findMin(eVal));
408 meshTools::writeOBJ(str, cM);
409 meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
410 meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
411 meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
413 for (label i = 1; i < 4; i++)
415 str <<
"l " << 1 <<
' ' << i + 1 <<
endl;
Istream and Ostream manipulators taking arguments.
Various functions to operate on Lists.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Output to file stream, using an OSstream.
An ordered pair of two objects of type <T> with first() and second() elements.
Extract command arguments and options from the supplied argc and argv parameters.
T get(const label index) const
Get a value from the argument at index.
bool found(const word &optName) const
Return true if the named option is found.
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
A class for handling file names.
A triFace with additional (region) index.
Tensor of scalars, i.e. Tensor<scalar>.
Triangulated surface description with patch information.
A triangle primitive used to calculate face normals and swept volumes.
scalar mag() const
Return scalar magnitude.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define WarningInFunction
Report a warning using Foam::Warning.
bool match(const UList< wordRe > &patterns, const std::string &text)
Return true if text matches one of the regular expressions.
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
label findMin(const ListType &input, label start=0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setprecision(const int i)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.