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);
173 principal[0] = eVec.x();
174 principal[1] = eVec.y();
175 principal[2] = eVec.z();
179 scalar maxMagDotProduct = -GREAT;
185 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
187 if (maxMagDotProduct < magDotProduct)
189 maxMagDotProduct = magDotProduct;
199 cartesian[
match.first()] & principal[
match.second()]
209 tPrincipal[
match.second()] *= -1;
211 tPrincipal[(
match.second() + 1) % 3] =
212 principal[(
match.second() + 2) % 3];
214 tPrincipal[(
match.second() + 2) % 3] =
215 principal[(
match.second() + 1) % 3];
217 principal = tPrincipal;
221 tEVal[(
match.second() + 1) % 3] = eVal[(
match.second() + 2) % 3];
223 tEVal[(
match.second() + 2) % 3] = eVal[(
match.second() + 1) % 3];
228 label permutationDelta =
match.second() -
match.first();
230 if (permutationDelta != 0)
234 permutationDelta += 3;
240 for (label i = 0; i < 3; i++)
242 tPrincipal[i] = principal[(i + permutationDelta) % 3];
244 tEVal[i] = eVal[(i + permutationDelta) % 3];
247 principal = tPrincipal;
252 const label matchedAlready =
match.first();
255 maxMagDotProduct = -GREAT;
259 if (cI == matchedAlready)
266 if (pI == matchedAlready)
271 scalar magDotProduct =
mag(cartesian[cI] & principal[pI]);
273 if (maxMagDotProduct < magDotProduct)
275 maxMagDotProduct = magDotProduct;
285 cartesian[
match.first()] & principal[
match.second()]
288 if (sense < 0 || (
match.second() -
match.first()) != 0)
290 principal[
match.second()] *= -1;
294 tPrincipal[(matchedAlready + 1) % 3] =
295 principal[(matchedAlready + 2) % 3]*-sense;
297 tPrincipal[(matchedAlready + 2) % 3] =
298 principal[(matchedAlready + 1) % 3]*-sense;
300 principal = tPrincipal;
304 tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
306 tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
311 eVec =
tensor(principal[0], principal[1], principal[2]);
326 <<
"Non-unique eigenvectors, cannot compute transformation " 327 <<
"from Cartesian axes" <<
endl;
329 showTransform =
false;
334 scalar surfaceArea = 0;
344 <<
"Illegal triangle " << facei <<
" vertices " <<
f 345 <<
" coords " <<
f.points(surf.points()) <<
endl;
351 surfaceArea += tri.
mag();
357 averageNormal.normalise();
361 <<
"Density: " << density <<
nl 362 <<
"Mass: " << m <<
nl 363 <<
"Centre of mass: " << cM <<
nl 364 <<
"Surface area: " << surfaceArea <<
nl 365 <<
"Average normal: " << averageNormal <<
nl 366 <<
"Inertia tensor around centre of mass: " <<
nl << J <<
nl 367 <<
"eigenValues (principal moments): " << eVal <<
nl 368 <<
"eigenVectors (principal axes): " <<
nl 369 << eVec.x() <<
nl << eVec.y() <<
nl << eVec.z() <<
endl;
373 Info<<
"Transform tensor from reference state (orientation):" <<
nl 375 <<
"Rotation tensor required to transform " 376 "from the body reference frame to the global " 377 "reference frame, i.e.:" <<
nl 378 <<
"globalVector = orientation & bodyLocalVector" 382 <<
"Entries for sixDoFRigidBodyDisplacement boundary condition:" 386 .writeEntry(
"mass", m)
387 .writeEntry(
"centreOfMass", cM)
388 .writeEntry(
"momentOfInertia", eVal)
389 .writeEntry(
"orientation", eVec.T());
394 Info<<
nl <<
"Inertia tensor relative to " << refPt <<
": " <<
nl 404 Info<<
nl <<
"Writing scaled principal axes at centre of mass of " 405 << surfFileName <<
" to " << str.name() <<
endl;
408 mag(cM - surf.points()[0])/eVal.component(
findMin(eVal));
415 for (label i = 1; i < 4; i++)
417 str <<
"l " << 1 <<
' ' << i + 1 <<
endl;
dimensionedScalar sign(const dimensionedScalar &ds)
static void addNote(const string ¬e)
Add extra notes for the usage information.
A triangle primitive used to calculate face normals and swept volumes. Uses referred points...
A class for handling file names.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A 1D vector of objects of type <T> with a fixed length <N>.
scalar mag() const
The magnitude of the triangle area.
Output to file stream as an OSstream, normally using std::ofstream for the actual output...
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.
bool match(const UList< wordRe > &selectors, const std::string &text)
True if text matches one of the selector expressions.
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)
static vector areaNormal(const Point &p0, const Point &p1, const Point &p2)
The area normal for a triangle defined by three points (right-hand rule). Magnitude equal to area of ...
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.
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)