momentOfInertia.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 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 \*---------------------------------------------------------------------------*/
25 
26 #include "momentOfInertia.H"
28 
29 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
30 
32 (
33  const pointField& pts,
34  const triFaceList& triFaces,
35  scalar density,
36  scalar& mass,
37  vector& cM,
38  tensor& J
39 )
40 {
41  // Reimplemented from: Wm4PolyhedralMassProperties.cpp
42  // File Version: 4.10.0 (2009/11/18)
43 
44  // Geometric Tools, LC
45  // Copyright (c) 1998-2010
46  // Distributed under the Boost Software License, Version 1.0.
47  // http://www.boost.org/LICENSE_1_0.txt
48  // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
49 
50  // Boost Software License - Version 1.0 - August 17th, 2003
51 
52  // Permission is hereby granted, free of charge, to any person or
53  // organization obtaining a copy of the software and accompanying
54  // documentation covered by this license (the "Software") to use,
55  // reproduce, display, distribute, execute, and transmit the
56  // Software, and to prepare derivative works of the Software, and
57  // to permit third-parties to whom the Software is furnished to do
58  // so, all subject to the following:
59 
60  // The copyright notices in the Software and this entire
61  // statement, including the above license grant, this restriction
62  // and the following disclaimer, must be included in all copies of
63  // the Software, in whole or in part, and all derivative works of
64  // the Software, unless such copies or derivative works are solely
65  // in the form of machine-executable object code generated by a
66  // source language processor.
67 
68  // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
69  // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
70  // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
71  // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
72  // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
73  // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
74  // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
75  // USE OR OTHER DEALINGS IN THE SOFTWARE.
76 
77  const scalar r6 = 1.0/6.0;
78  const scalar r24 = 1.0/24.0;
79  const scalar r60 = 1.0/60.0;
80  const scalar r120 = 1.0/120.0;
81 
82  // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
83  scalarField integrals(10, 0.0);
84 
85  forAll(triFaces, i)
86  {
87  const triFace& tri(triFaces[i]);
88 
89  // vertices of triangle i
90  vector v0 = pts[tri[0]];
91  vector v1 = pts[tri[1]];
92  vector v2 = pts[tri[2]];
93 
94  // cross product of edges
95  vector eA = v1 - v0;
96  vector eB = v2 - v0;
97  vector n = eA ^ eB;
98 
99  // compute integral terms
100  scalar tmp0, tmp1, tmp2;
101 
102  scalar f1x, f2x, f3x, g0x, g1x, g2x;
103 
104  tmp0 = v0.x() + v1.x();
105  f1x = tmp0 + v2.x();
106  tmp1 = v0.x()*v0.x();
107  tmp2 = tmp1 + v1.x()*tmp0;
108  f2x = tmp2 + v2.x()*f1x;
109  f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
110  g0x = f2x + v0.x()*(f1x + v0.x());
111  g1x = f2x + v1.x()*(f1x + v1.x());
112  g2x = f2x + v2.x()*(f1x + v2.x());
113 
114  scalar f1y, f2y, f3y, g0y, g1y, g2y;
115 
116  tmp0 = v0.y() + v1.y();
117  f1y = tmp0 + v2.y();
118  tmp1 = v0.y()*v0.y();
119  tmp2 = tmp1 + v1.y()*tmp0;
120  f2y = tmp2 + v2.y()*f1y;
121  f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
122  g0y = f2y + v0.y()*(f1y + v0.y());
123  g1y = f2y + v1.y()*(f1y + v1.y());
124  g2y = f2y + v2.y()*(f1y + v2.y());
125 
126  scalar f1z, f2z, f3z, g0z, g1z, g2z;
127 
128  tmp0 = v0.z() + v1.z();
129  f1z = tmp0 + v2.z();
130  tmp1 = v0.z()*v0.z();
131  tmp2 = tmp1 + v1.z()*tmp0;
132  f2z = tmp2 + v2.z()*f1z;
133  f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
134  g0z = f2z + v0.z()*(f1z + v0.z());
135  g1z = f2z + v1.z()*(f1z + v1.z());
136  g2z = f2z + v2.z()*(f1z + v2.z());
137 
138  // update integrals
139  integrals[0] += n.x()*f1x;
140  integrals[1] += n.x()*f2x;
141  integrals[2] += n.y()*f2y;
142  integrals[3] += n.z()*f2z;
143  integrals[4] += n.x()*f3x;
144  integrals[5] += n.y()*f3y;
145  integrals[6] += n.z()*f3z;
146  integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
147  integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
148  integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
149  }
150 
151  integrals[0] *= r6;
152  integrals[1] *= r24;
153  integrals[2] *= r24;
154  integrals[3] *= r24;
155  integrals[4] *= r60;
156  integrals[5] *= r60;
157  integrals[6] *= r60;
158  integrals[7] *= r120;
159  integrals[8] *= r120;
160  integrals[9] *= r120;
161 
162  // mass
163  mass = integrals[0];
164 
165  // center of mass
166  cM = vector(integrals[1], integrals[2], integrals[3])/mass;
167 
168  // inertia relative to origin
169  J.xx() = integrals[5] + integrals[6];
170  J.xy() = -integrals[7];
171  J.xz() = -integrals[9];
172  J.yx() = J.xy();
173  J.yy() = integrals[4] + integrals[6];
174  J.yz() = -integrals[8];
175  J.zx() = J.xz();
176  J.zy() = J.yz();
177  J.zz() = integrals[4] + integrals[5];
178 
179  // inertia relative to center of mass
180  J -= mass*((cM & cM)*I - cM*cM);
181 
182  // Apply density
183  mass *= density;
184  J *= density;
185 }
186 
187 
189 (
190  const pointField& pts,
191  const triFaceList& triFaces,
192  scalar density,
193  scalar& mass,
194  vector& cM,
195  tensor& J
196 )
197 {
198  // Reset properties for accumulation
199 
200  mass = 0.0;
201  cM = vector::zero;
202  J = tensor::zero;
203 
204  // Find centre of mass
205 
206  forAll(triFaces, i)
207  {
208  const triFace& tri(triFaces[i]);
209 
210  triPointRef t
211  (
212  pts[tri[0]],
213  pts[tri[1]],
214  pts[tri[2]]
215  );
216 
217  scalar triMag = t.mag();
218 
219  cM += triMag*t.centre();
220 
221  mass += triMag;
222  }
223 
224  cM /= mass;
225 
226  mass *= density;
227 
228  // Find inertia around centre of mass
229 
230  forAll(triFaces, i)
231  {
232  const triFace& tri(triFaces[i]);
233 
234  J += triPointRef
235  (
236  pts[tri[0]],
237  pts[tri[1]],
238  pts[tri[2]]
239  ).inertia(cM, density);
240  }
241 }
242 
243 
245 (
246  const triSurface& surf,
247  scalar density,
248  scalar& mass,
249  vector& cM,
250  tensor& J
251 )
252 {
253  triFaceList faces(surf.size());
254 
255  forAll(surf, i)
256  {
257  faces[i] = triFace(surf[i]);
258  }
259 
260  massPropertiesSolid(surf.points(), faces, density, mass, cM, J);
261 }
262 
263 
265 (
266  const triSurface& surf,
267  scalar density,
268  scalar& mass,
269  vector& cM,
270  tensor& J
271 )
272 {
273  triFaceList faces(surf.size());
274 
275  forAll(surf, i)
276  {
277  faces[i] = triFace(surf[i]);
278  }
279 
280  massPropertiesShell(surf.points(), faces, density, mass, cM, J);
281 }
282 
283 
285 (
286  scalar mass,
287  const vector& cM,
288  const tensor& J,
289  const vector& refPt
290 )
291 {
292  // The displacement vector (refPt = cM) is the displacement of the
293  // new reference point from the centre of mass of the body that
294  // the inertia tensor applies to.
295 
296  vector d = (refPt - cM);
297 
298  return J + mass*((d & d)*I - d*d);
299 }
300 
301 
303 (
304  const polyMesh& mesh
305 )
306 {
308 
309  tensorField& tf = tTf();
310 
311  forAll(tf, cI)
312  {
313  tf[cI] = meshInertia(mesh, cI);
314  }
315 
316  return tTf;
317 }
318 
319 
321 (
322  const polyMesh& mesh,
323  label cellI
324 )
325 {
326  List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
327  (
328  mesh,
329  cellI
330  );
331 
332  triFaceList faces(cellTets.size());
333 
334  forAll(cellTets, cTI)
335  {
336  faces[cTI] = cellTets[cTI].faceTriIs(mesh);
337  }
338 
339  scalar m = 0.0;
340  vector cM = vector::zero;
341  tensor J = tensor::zero;
342 
343  massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J);
344 
345  return J;
346 }
347 
348 
349 // ************************************************************************* //
const Cmpt & xy() const
Definition: TensorI.H:167
static tmp< tensorField > meshInertia(const polyMesh &mesh)
const Cmpt & zz() const
Definition: TensorI.H:216
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:96
static const sphericalTensor I(1)
Triangulated surface description with patch information.
Definition: triSurface.H:57
static void massPropertiesSolid(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const Cmpt & yy() const
Definition: TensorI.H:188
label nCells() const
face triFace(3)
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
label n
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
const Cmpt & y() const
Definition: VectorI.H:71
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
#define forAll(list, i)
Definition: UList.H:421
const Cmpt & x() const
Definition: VectorI.H:65
const Cmpt & yx() const
Definition: TensorI.H:181
const Cmpt & z() const
Definition: VectorI.H:77
Point centre() const
Return centre (centroid)
Definition: triangleI.H:89
const tensorField & tf
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static tensor applyParallelAxisTheorem(scalar mass, const vector &cM, const tensor &J, const vector &refPt)
const Cmpt & zy() const
Definition: TensorI.H:209
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const Cmpt & xx() const
Definition: TensorI.H:160
const Cmpt & yz() const
Definition: TensorI.H:195
const Cmpt & xz() const
Definition: TensorI.H:174
const Cmpt & zx() const
Definition: TensorI.H:202
A class for managing temporary objects.
Definition: PtrList.H:118
const Field< PointType > & points() const
Return reference to global points.