56 int 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." 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" 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 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));
413 for (label i = 1; i < 4; i++)
415 str <<
"l " << 1 <<
' ' << i + 1 <<
endl;
dimensionedScalar sign(const dimensionedScalar &ds)
bool match(const UList< wordRe > &patterns, const std::string &text)
Return true if text matches one of the regular expressions.
static void addNote(const string ¬e)
Add extra notes for the usage information.
A class for handling file names.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Output to file stream, using an OSstream.
constexpr char nl
The newline '\n' character (0x0a)
static void massPropertiesSolid(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static void addBoolOption(const word &optName, const string &usage="", bool advanced=false)
Add a bool option to validOptions with usage information.
static void noParallel()
Remove the parallel options.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Various functions to operate on Lists.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
#define forAll(list, i)
Loop across all elements in list.
label findMin(const ListType &input, label start=0)
Linear search for the index of the min element, similar to std::min_element but for lists and returns...
Omanip< int > setprecision(const int i)
Extract command arguments and options from the supplied argc and argv parameters. ...
static void addOption(const word &optName, const string ¶m="", const string &usage="", bool advanced=false)
Add an option to validOptions with usage information.
A triFace with additional (region) index.
Istream and Ostream manipulators taking arguments.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
static tensor applyParallelAxisTheorem(scalar mass, const vector &cM, const tensor &J, const vector &refPt)
T get(const label index) const
Get a value from the argument at index.
#define WarningInFunction
Report a warning using Foam::Warning.
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
static void addArgument(const string &argName, const string &usage="")
Append a (mandatory) argument to validArgs.
triangle< point, const point & > triPointRef
A triangle using referred points.
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Triangulated surface description with patch information.
Foam::argList args(argc, argv)
Tensor of scalars, i.e. Tensor<scalar>.
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J, bool doReduce=false)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
bool found(const word &optName) const
Return true if the named option is found.
static constexpr const zero Zero
Global zero (0)