surfaceInertia.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  surfaceInertia
26 
27 Description
28  Calculates the inertia tensor, principal axes and moments of a
29  command line specified triSurface. Inertia can either be of the
30  solid body or of a thin shell.
31 
32 \*---------------------------------------------------------------------------*/
33 
34 #include "argList.H"
35 #include "ListOps.H"
36 #include "triSurface.H"
37 #include "OFstream.H"
38 #include "meshTools.H"
39 #include "Random.H"
40 #include "transform.H"
41 #include "IOmanip.H"
42 #include "Pair.H"
43 #include "momentOfInertia.H"
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
47 using namespace Foam;
48 
49 int main(int argc, char *argv[])
50 {
52  (
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."
56  );
57 
59  argList::validArgs.append("surfaceFile");
61  (
62  "shellProperties",
63  "inertia of a thin shell"
64  );
65 
67  (
68  "density",
69  "scalar",
70  "Specify density, "
71  "kg/m3 for solid properties, kg/m2 for shell properties"
72  );
73 
75  (
76  "referencePoint",
77  "vector",
78  "Inertia relative to this point, not the centre of mass"
79  );
80 
81  argList args(argc, argv);
82 
83  const fileName surfFileName = args[1];
84  const scalar density = args.optionLookupOrDefault("density", 1.0);
85 
86  vector refPt = Zero;
87  bool calcAroundRefPt = args.optionReadIfPresent("referencePoint", refPt);
88 
89  triSurface surf(surfFileName);
90 
91  scalar m = 0.0;
92  vector cM = Zero;
93  tensor J = Zero;
94 
95  if (args.optionFound("shellProperties"))
96  {
97  momentOfInertia::massPropertiesShell(surf, density, m, cM, J);
98  }
99  else
100  {
101  momentOfInertia::massPropertiesSolid(surf, density, m, cM, J);
102  }
103 
104  if (m < 0)
105  {
107  << "Negative mass detected, the surface may be inside-out." << endl;
108  }
109 
110  vector eVal = eigenValues(J);
111 
112  tensor eVec = eigenVectors(J);
113 
114  label pertI = 0;
115 
116  Random rand(57373);
117 
118  while ((magSqr(eVal) < VSMALL) && pertI < 10)
119  {
121  << "No eigenValues found, shape may have symmetry, "
122  << "perturbing inertia tensor diagonal" << endl;
123 
124  J.xx() *= 1.0 + SMALL*rand.scalar01();
125  J.yy() *= 1.0 + SMALL*rand.scalar01();
126  J.zz() *= 1.0 + SMALL*rand.scalar01();
127 
128  eVal = eigenValues(J);
129 
130  eVec = eigenVectors(J);
131 
132  pertI++;
133  }
134 
135  bool showTransform = true;
136 
137  if
138  (
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))
142  )
143  {
144  // Make the eigenvectors a right handed orthogonal triplet
145  eVec = tensor
146  (
147  eVec.x(),
148  eVec.y(),
149  eVec.z() * sign((eVec.x() ^ eVec.y()) & eVec.z())
150  );
151 
152  // Finding the most natural transformation. Using Lists
153  // rather than tensors to allow indexed permutation.
154 
155  // Cartesian basis vectors - right handed orthogonal triplet
156  List<vector> cartesian(3);
157 
158  cartesian[0] = vector(1, 0, 0);
159  cartesian[1] = vector(0, 1, 0);
160  cartesian[2] = vector(0, 0, 1);
161 
162  // Principal axis basis vectors - right handed orthogonal
163  // triplet
164  List<vector> principal(3);
165 
166  principal[0] = eVec.x();
167  principal[1] = eVec.y();
168  principal[2] = eVec.z();
169 
170  scalar maxMagDotProduct = -GREAT;
171 
172  // Matching axis indices, first: cartesian, second:principal
173 
174  Pair<label> match(-1, -1);
175 
176  forAll(cartesian, cI)
177  {
178  forAll(principal, pI)
179  {
180  scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
181 
182  if (magDotProduct > maxMagDotProduct)
183  {
184  maxMagDotProduct = magDotProduct;
185 
186  match.first() = cI;
187 
188  match.second() = pI;
189  }
190  }
191  }
192 
193  scalar sense = sign
194  (
195  cartesian[match.first()] & principal[match.second()]
196  );
197 
198  if (sense < 0)
199  {
200  // Invert the best match direction and swap the order of
201  // the other two vectors
202 
203  List<vector> tPrincipal = principal;
204 
205  tPrincipal[match.second()] *= -1;
206 
207  tPrincipal[(match.second() + 1) % 3] =
208  principal[(match.second() + 2) % 3];
209 
210  tPrincipal[(match.second() + 2) % 3] =
211  principal[(match.second() + 1) % 3];
212 
213  principal = tPrincipal;
214 
215  vector tEVal = eVal;
216 
217  tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
218 
219  tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
220 
221  eVal = tEVal;
222  }
223 
224  label permutationDelta = match.second() - match.first();
225 
226  if (permutationDelta != 0)
227  {
228  // Add 3 to the permutationDelta to avoid negative indices
229 
230  permutationDelta += 3;
231 
232  List<vector> tPrincipal = principal;
233 
234  vector tEVal = eVal;
235 
236  for (label i = 0; i < 3; i++)
237  {
238  tPrincipal[i] = principal[(i + permutationDelta) % 3];
239 
240  tEVal[i] = eVal[(i + permutationDelta) % 3];
241  }
242 
243  principal = tPrincipal;
244 
245  eVal = tEVal;
246  }
247 
248  label matchedAlready = match.first();
249 
250  match =Pair<label>(-1, -1);
251 
252  maxMagDotProduct = -GREAT;
253 
254  forAll(cartesian, cI)
255  {
256  if (cI == matchedAlready)
257  {
258  continue;
259  }
260 
261  forAll(principal, pI)
262  {
263  if (pI == matchedAlready)
264  {
265  continue;
266  }
267 
268  scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
269 
270  if (magDotProduct > maxMagDotProduct)
271  {
272  maxMagDotProduct = magDotProduct;
273 
274  match.first() = cI;
275 
276  match.second() = pI;
277  }
278  }
279  }
280 
281  sense = sign
282  (
283  cartesian[match.first()] & principal[match.second()]
284  );
285 
286  if (sense < 0 || (match.second() - match.first()) != 0)
287  {
288  principal[match.second()] *= -1;
289 
290  List<vector> tPrincipal = principal;
291 
292  tPrincipal[(matchedAlready + 1) % 3] =
293  principal[(matchedAlready + 2) % 3]*-sense;
294 
295  tPrincipal[(matchedAlready + 2) % 3] =
296  principal[(matchedAlready + 1) % 3]*-sense;
297 
298  principal = tPrincipal;
299 
300  vector tEVal = eVal;
301 
302  tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
303 
304  tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
305 
306  eVal = tEVal;
307  }
308 
309  eVec = tensor(principal[0], principal[1], principal[2]);
310 
311  // {
312  // tensor R = rotationTensor(vector(1, 0, 0), eVec.x());
313 
314  // R = rotationTensor(R & vector(0, 1, 0), eVec.y()) & R;
315 
316  // Info<< "R = " << nl << R << endl;
317 
318  // Info<< "R - eVec.T() " << R - eVec.T() << endl;
319  // }
320  }
321  else
322  {
324  << "Non-unique eigenvectors, cannot compute transformation "
325  << "from Cartesian axes" << endl;
326 
327  showTransform = false;
328  }
329 
330  // calculate the total surface area
331 
332  scalar surfaceArea = 0;
333 
334  forAll(surf, facei)
335  {
336  const labelledTri& f = surf[facei];
337 
338  if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
339  {
341  << "Illegal triangle " << facei << " vertices " << f
342  << " coords " << f.points(surf.points()) << endl;
343  }
344  else
345  {
346  surfaceArea += triPointRef
347  (
348  surf.points()[f[0]],
349  surf.points()[f[1]],
350  surf.points()[f[2]]
351  ).mag();
352  }
353  }
354 
355  Info<< nl << setprecision(12)
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;
364 
365  if (showTransform)
366  {
367  Info<< "Transform tensor from reference state (orientation):" << nl
368  << eVec.T() << 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"
373  << endl;
374 
375  Info<< nl
376  << "Entries for sixDoFRigidBodyDisplacement boundary condition:"
377  << nl
378  << " mass " << m << token::END_STATEMENT << nl
379  << " centreOfMass " << cM << token::END_STATEMENT << nl
380  << " momentOfInertia " << eVal << token::END_STATEMENT << nl
381  << " orientation " << eVec.T() << token::END_STATEMENT
382  << endl;
383  }
384 
385  if (calcAroundRefPt)
386  {
387  Info<< nl << "Inertia tensor relative to " << refPt << ": " << nl
389  << endl;
390  }
391 
392  OFstream str("axes.obj");
393 
394  Info<< nl << "Writing scaled principal axes at centre of mass of "
395  << surfFileName << " to " << str.name() << endl;
396 
397  scalar scale = mag(cM - surf.points()[0])/eVal.component(findMin(eVal));
398 
399  meshTools::writeOBJ(str, cM);
400  meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
401  meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
402  meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
403 
404  for (label i = 1; i < 4; i++)
405  {
406  str << "l " << 1 << ' ' << i + 1 << endl;
407  }
408 
409  Info<< nl << "End" << nl << endl;
410 
411  return 0;
412 }
413 
414 
415 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
const Cmpt & z() const
Definition: VectorI.H:87
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const Cmpt & yy() const
Definition: TensorI.H:181
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A class for handling file names.
Definition: fileName.H:69
const Cmpt & x() const
Definition: VectorI.H:75
Output to file stream.
Definition: OFstream.H:81
bool optionReadIfPresent(const word &opt, T &) const
Read a value from the named option if present.
Definition: argListI.H:198
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.
Definition: Ostream.H:253
static void noParallel()
Remove the parallel options.
Definition: argList.C:146
const Type & first() const
Return first.
Definition: Pair.H:87
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:154
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Vector< Cmpt > x() const
Definition: TensorI.H:279
Vector< Cmpt > y() const
Definition: TensorI.H:286
T optionLookupOrDefault(const word &opt, const T &deflt) const
Read a value from the named option if present.
Definition: argListI.H:237
Various functions to operate on Lists.
Omanip< int > setprecision(const int i)
Definition: IOmanip.H:205
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 /)
Definition: fileName.C:179
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
3D tensor transformation operations.
Extract command arguments and options from the supplied argc and argv parameters. ...
Definition: argList.H:102
static void addOption(const word &opt, const string &param="", const string &usage="")
Add to an option to validOptions with usage information.
Definition: argList.C:93
const Cmpt & xx() const
Definition: TensorI.H:153
const Cmpt & y() const
Definition: VectorI.H:81
static const zero Zero
Definition: zero.H:91
Triangle with additional region number.
Definition: labelledTri.H:57
Vector< Cmpt > z() const
Definition: TensorI.H:293
Simple random number generator.
Definition: Random.H:49
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:262
Tensor< Cmpt > T() const
Return transpose.
Definition: TensorI.H:321
const Cmpt & zz() const
Definition: TensorI.H:209
dimensionedTensor eigenVectors(const dimensionedTensor &dt)
pointField points(const pointField &) const
Return the points corresponding to this face.
Definition: triFaceI.H:128
const Cmpt & component(const direction) const
Definition: VectorSpaceI.H:101
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
Definition: triPointRef.H:44
messageStream Info
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.
Definition: argList.C:83
static void addNote(const string &)
Add extra notes for the usage information.
Definition: argList.C:124
Triangulated surface description with patch information.
Definition: triSurface.H:65
Foam::argList args(argc, argv)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Namespace for OpenFOAM.