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;
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;
Istream and Ostream manipulators taking arguments.
Various functions to operate on Lists.
#define forAll(list, i)
Loop across all elements in list.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Tensor< Cmpt > T() const
Return transpose.
const Cmpt & component(const direction) const
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.
static void addNote(const string &)
Add extra notes for the usage information.
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
bool optionFound(const word &opt) const
Return true if the named option is found.
bool optionReadIfPresent(const word &opt, T &) const
Read a value from the named option if present.
static SLList< string > validArgs
A list of valid (mandatory) arguments.
T optionLookupOrDefault(const word &opt, const T &deflt) const
Read a value from the named option if present.
A class for handling file names.
word name() const
Return file name (part beyond last /)
Triangle with additional region number.
static tensor applyParallelAxisTheorem(scalar mass, const vector &cM, const tensor &J, const vector &refPt)
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
static void massPropertiesSolid(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
scalar scalar01()
Return a scalar uniformly distributed between zero and one.
Triangulated surface description with patch information.
int main(int argc, char *argv[])
#define WarningInFunction
Report a warning using Foam::Warning.
Omanip< int > setprecision(const int i)
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Tensor< scalar > tensor
Tensor of scalars.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Vector< scalar > vector
A scalar version of the templated Vector.
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensionedVector eigenValues(const dimensionedTensor &dt)
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
triangle< point, const point & > triPointRef
dimensioned< scalar > magSqr(const dimensioned< Type > &)
randomGenerator rndGen(653213)
Foam::argList args(argc, argv)