49 int main(
int argc,
char *argv[])
53 "Calculates the inertia tensor and principal axes and moments " 54 "of the specified surface.\n" 55 "Inertia can either be of the solid body or of a thin shell." 63 "inertia of a thin shell" 71 "kg/m3 for solid properties, kg/m2 for shell properties" 78 "Inertia relative to this point, not the centre of mass" 107 <<
"Negative mass detected, the surface may be inside-out." <<
endl;
118 while ((
magSqr(eVal) < VSMALL) && pertI < 10)
121 <<
"No eigenValues found, shape may have symmetry, " 122 <<
"perturbing inertia tensor diagonal" <<
endl;
124 J.
xx() *= 1.0 + SMALL*rand.scalar01();
125 J.
yy() *= 1.0 + SMALL*rand.scalar01();
126 J.
zz() *= 1.0 + SMALL*rand.scalar01();
135 bool showTransform =
true;
139 (
mag(eVec.
x() ^ eVec.
y()) > (1.0 - SMALL))
140 && (
mag(eVec.
y() ^ eVec.
z()) > (1.0 - SMALL))
141 && (
mag(eVec.
z() ^ eVec.
x()) > (1.0 - SMALL))
149 eVec.
z() *
sign((eVec.
x() ^ eVec.
y()) & eVec.
z())
158 cartesian[0] =
vector(1, 0, 0);
159 cartesian[1] =
vector(0, 1, 0);
160 cartesian[2] =
vector(0, 0, 1);
166 principal[0] = eVec.
x();
167 principal[1] = eVec.
y();
168 principal[2] = eVec.
z();
170 scalar maxMagDotProduct = -GREAT;
180 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
182 if (magDotProduct > maxMagDotProduct)
184 maxMagDotProduct = magDotProduct;
195 cartesian[match.first()] & principal[match.second()]
205 tPrincipal[match.second()] *= -1;
207 tPrincipal[(match.second() + 1) % 3] =
208 principal[(match.second() + 2) % 3];
210 tPrincipal[(match.second() + 2) % 3] =
211 principal[(match.second() + 1) % 3];
213 principal = tPrincipal;
217 tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
219 tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
224 label permutationDelta = match.second() - match.first();
226 if (permutationDelta != 0)
230 permutationDelta += 3;
236 for (
label i = 0; i < 3; i++)
238 tPrincipal[i] = principal[(i + permutationDelta) % 3];
240 tEVal[i] = eVal[(i + permutationDelta) % 3];
243 principal = tPrincipal;
248 label matchedAlready = match.first();
252 maxMagDotProduct = -GREAT;
256 if (cI == matchedAlready)
263 if (pI == matchedAlready)
268 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
270 if (magDotProduct > maxMagDotProduct)
272 maxMagDotProduct = magDotProduct;
283 cartesian[match.first()] & principal[match.second()]
286 if (sense < 0 || (match.second() - match.first()) != 0)
288 principal[match.second()] *= -1;
292 tPrincipal[(matchedAlready + 1) % 3] =
293 principal[(matchedAlready + 2) % 3]*-sense;
295 tPrincipal[(matchedAlready + 2) % 3] =
296 principal[(matchedAlready + 1) % 3]*-sense;
298 principal = tPrincipal;
302 tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
304 tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
309 eVec =
tensor(principal[0], principal[1], principal[2]);
324 <<
"Non-unique eigenvectors, cannot compute transformation " 325 <<
"from Cartesian axes" <<
endl;
327 showTransform =
false;
332 scalar surfaceArea = 0;
338 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
341 <<
"Illegal triangle " << facei <<
" vertices " << f
342 <<
" coords " << f.
points(surf.points()) << endl;
356 <<
"Density: " << density <<
nl 357 <<
"Mass: " << m <<
nl 358 <<
"Centre of mass: " << cM <<
nl 359 <<
"Surface area: " << surfaceArea <<
nl 360 <<
"Inertia tensor around centre of mass: " <<
nl << J <<
nl 361 <<
"eigenValues (principal moments): " << eVal <<
nl 362 <<
"eigenVectors (principal axes): " <<
nl 363 << eVec.
x() <<
nl << eVec.
y() <<
nl << eVec.
z() <<
endl;
367 Info<<
"Transform tensor from reference state (orientation):" <<
nl 369 <<
"Rotation tensor required to transform " 370 "from the body reference frame to the global " 371 "reference frame, i.e.:" <<
nl 372 <<
"globalVector = orientation & bodyLocalVector" 376 <<
"Entries for sixDoFRigidBodyDisplacement boundary condition:" 387 Info<<
nl <<
"Inertia tensor relative to " << refPt <<
": " <<
nl 394 Info<<
nl <<
"Writing scaled principal axes at centre of mass of " 395 << surfFileName <<
" to " << str.
name() <<
endl;
404 for (
label i = 1; i < 4; i++)
406 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.
bool optionReadIfPresent(const word &opt, T &) const
Read a value from the named option if present.
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.
static void noParallel()
Remove the parallel options.
const Type & first() const
Return first.
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.
T optionLookupOrDefault(const word &opt, const T &deflt) const
Read a value from the named option if present.
Various functions to operate on Lists.
Omanip< int > setprecision(const int i)
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
word name() const
Return file name (part beyond last /)
An ordered pair of two objects of type <T> with first() and second() elements.
bool optionFound(const word &opt) const
Return true if the named option is found.
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.
Triangle with additional region number.
Simple random number generator.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
Tensor< Cmpt > T() const
Return transpose.
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
pointField points(const pointField &) const
Return the points corresponding to this face.
const Cmpt & component(const direction) const
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
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)
Tensor< scalar > tensor
Tensor of scalars.