49 int main(
int argc,
char *argv[])
55 "Calculates the inertia tensor and principal axes and moments " 56 "of the specified surface.\n" 57 "Inertia can either be of the solid body or of a thin shell." 64 "inertia of a thin shell" 72 "kg/m3 for solid properties, kg/m2 for shell properties" 79 "Inertia relative to this point, not the centre of mass" 108 <<
"Negative mass detected, the surface may be inside-out." <<
endl;
119 while ((
magSqr(eVal) < vSmall) && pertI < 10)
122 <<
"No eigenValues found, shape may have symmetry, " 123 <<
"perturbing inertia tensor diagonal" <<
endl;
125 J.
xx() *= 1.0 + small*rand.scalar01();
126 J.
yy() *= 1.0 + small*rand.scalar01();
127 J.
zz() *= 1.0 + small*rand.scalar01();
136 bool showTransform =
true;
140 (
mag(eVec.
x() ^ eVec.
y()) > (1.0 - small))
141 && (
mag(eVec.
y() ^ eVec.
z()) > (1.0 - small))
142 && (
mag(eVec.
z() ^ eVec.
x()) > (1.0 - small))
150 eVec.
z() *
sign((eVec.
x() ^ eVec.
y()) & eVec.
z())
159 cartesian[0] =
vector(1, 0, 0);
160 cartesian[1] =
vector(0, 1, 0);
161 cartesian[2] =
vector(0, 0, 1);
167 principal[0] = eVec.
x();
168 principal[1] = eVec.
y();
169 principal[2] = eVec.
z();
171 scalar maxMagDotProduct = -great;
181 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
183 if (magDotProduct > maxMagDotProduct)
185 maxMagDotProduct = magDotProduct;
196 cartesian[match.first()] & principal[match.second()]
206 tPrincipal[match.second()] *= -1;
208 tPrincipal[(match.second() + 1) % 3] =
209 principal[(match.second() + 2) % 3];
211 tPrincipal[(match.second() + 2) % 3] =
212 principal[(match.second() + 1) % 3];
214 principal = tPrincipal;
218 tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
220 tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
225 label permutationDelta = match.second() - match.first();
227 if (permutationDelta != 0)
231 permutationDelta += 3;
237 for (
label i = 0; i < 3; i++)
239 tPrincipal[i] = principal[(i + permutationDelta) % 3];
241 tEVal[i] = eVal[(i + permutationDelta) % 3];
244 principal = tPrincipal;
249 label matchedAlready = match.first();
253 maxMagDotProduct = -great;
257 if (cI == matchedAlready)
264 if (pI == matchedAlready)
269 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
271 if (magDotProduct > maxMagDotProduct)
273 maxMagDotProduct = magDotProduct;
284 cartesian[match.first()] & principal[match.second()]
287 if (sense < 0 || (match.second() - match.first()) != 0)
289 principal[match.second()] *= -1;
293 tPrincipal[(matchedAlready + 1) % 3] =
294 principal[(matchedAlready + 2) % 3]*-sense;
296 tPrincipal[(matchedAlready + 2) % 3] =
297 principal[(matchedAlready + 1) % 3]*-sense;
299 principal = tPrincipal;
303 tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
305 tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
310 eVec =
tensor(principal[0], principal[1], principal[2]);
325 <<
"Non-unique eigenvectors, cannot compute transformation " 326 <<
"from Cartesian axes" <<
endl;
328 showTransform =
false;
333 scalar surfaceArea = 0;
339 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
342 <<
"Illegal triangle " << facei <<
" vertices " << f
343 <<
" coords " << f.
points(surf.points()) << endl;
357 <<
"Density: " << density <<
nl 358 <<
"Mass: " << m <<
nl 359 <<
"Centre of mass: " << cM <<
nl 360 <<
"Surface area: " << surfaceArea <<
nl 361 <<
"Inertia tensor around centre of mass: " <<
nl << J <<
nl 362 <<
"eigenValues (principal moments): " << eVal <<
nl 363 <<
"eigenVectors (principal axes): " <<
nl 364 << eVec.
x() <<
nl << eVec.
y() <<
nl << eVec.
z() <<
endl;
368 Info<<
"Transform tensor from reference state (orientation):" <<
nl 370 <<
"Rotation tensor required to transform " 371 "from the body reference frame to the global " 372 "reference frame, i.e.:" <<
nl 373 <<
"globalVector = orientation & bodyLocalVector" 377 <<
"Entries for sixDoFRigidBodyDisplacement boundary condition:" 388 Info<<
nl <<
"Inertia tensor relative to " << refPt <<
": " <<
nl 395 Info<<
nl <<
"Writing scaled principal axes at centre of mass of " 396 << surfFileName <<
" to " << str.
name() <<
endl;
405 for (
label i = 1; i < 4; i++)
407 str <<
"l " << 1 <<
' ' << i + 1 <<
endl;
dimensionedScalar sign(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
Tensor< Cmpt > T() const
Return transpose.
static void massPropertiesSolid(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
dimensionedVector eigenValues(const dimensionedTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool optionFound(const word &opt) const
Return true if the named option is found.
static SLList< string > validArgs
A list of valid (mandatory) arguments.
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
Vector< scalar > vector
A scalar version of the templated Vector.
bool optionReadIfPresent(const word &opt, T &) const
Read a value from the named option if present.
Various functions to operate on Lists.
T optionLookupOrDefault(const word &opt, const T &deflt) const
Read a value from the named option if present.
Omanip< int > setprecision(const int i)
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
An ordered pair of two objects of type <T> with first() and second() elements.
pointField points(const pointField &) const
Return the points corresponding to this face.
Extract command arguments and options from the supplied argc and argv parameters. ...
static void addOption(const word &opt, const string ¶m="", const string &usage="")
Add to an option to validOptions with usage information.
word name() const
Return file name (part beyond last /)
Triangle with additional region number.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
static tensor applyParallelAxisTheorem(scalar mass, const vector &cM, const tensor &J, const vector &refPt)
#define WarningInFunction
Report a warning using Foam::Warning.
triangle< point, const point & > triPointRef
const Cmpt & component(const direction) const
dimensioned< scalar > mag(const dimensioned< Type > &)
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
static void addNote(const string &)
Add extra notes for the usage information.
Triangulated surface description with patch information.
Foam::argList args(argc, argv)
const Type & first() const
Return first.
Tensor< scalar > tensor
Tensor of scalars.