src/meshTools/momentOfInertia/momentOfInertia.C

   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"
  27 #include "polyMeshTetDecomposition.H"
  28
  29 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
  30
  31 void Foam::momentOfInertia::massPropertiesSolid
  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
 188 void Foam::momentOfInertia::massPropertiesShell
 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
 244 void Foam::momentOfInertia::massPropertiesSolid
 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
 264 void Foam::momentOfInertia::massPropertiesShell
 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
 284 Foam::tensor Foam::momentOfInertia::applyParallelAxisTheorem
 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
 302 Foam::tmp<Foam::tensorField> Foam::momentOfInertia::meshInertia
 303 (
 304     const polyMesh& mesh
 305 )
 306 {
 307     tmp<tensorField> tTf = tmp<tensorField>(new tensorField(mesh.nCells()));
 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
 320 Foam::tensor Foam::momentOfInertia::meshInertia
 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 // ************************************************************************* //