src/tetFiniteElement/tetFemCalculus/tetFec.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Class of static functions to calculate explicit finite element derivatives.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "tetCell.H"
  30 #include "tetCellList.H"
  31 #include "tetPointRef.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  39
  40
  41 template<class Type>
  42 tmp
  43 <
  44     GeometricField
  45     <
  46         typename outerProduct<vector, Type>::type,
  47         tetPolyPatchField,
  48         tetPointMesh
  49     >
  50 >
  51 tetFec::grad
  52 (
  53     const GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi
  54 )
  55 {
  56     typedef typename outerProduct<vector, Type>::type GradType;
  57
  58     const tetPolyMesh& tetMesh = psi.mesh();
  59
  60     tmp<GeometricField<GradType, tetPolyPatchField, tetPointMesh> > tFemGrad
  61     (
  62         new GeometricField<GradType, tetPolyPatchField, tetPointMesh>
  63         (
  64             IOobject
  65             (
  66                 "grad("+psi.name()+')',
  67                 psi.instance(),
  68                 tetMesh(),
  69                 IOobject::NO_READ,
  70                 IOobject::NO_WRITE
  71             ),
  72             tetMesh,
  73             dimensioned<GradType>
  74             (
  75                 "zero",
  76                 psi.dimensions()/dimLength,
  77                 pTraits<GradType>::zero
  78             )
  79         )
  80     );
  81
  82     GeometricField<GradType, tetPolyPatchField, tetPointMesh>& femGrad =
  83         tFemGrad();
  84
  85     pointField points = tetMesh.points();
  86     cellShapeList tetCellShapes = tetMesh.tetCells();
  87
  88     scalarField weights (points.size(), 0.0);
  89
  90     forAll (tetCellShapes, tetI)
  91     {
  92         cellShape& curShape = tetCellShapes[tetI];
  93
  94         tetPointRef curTetrahedron
  95         (
  96             points[curShape[0]],
  97             points[curShape[1]],
  98             points[curShape[2]],
  99             points[curShape[3]]
 100         );
 101
 102         GradType tetGrad =
 103           - (1.0/3.0)*
 104             (
 105                 curTetrahedron.Sa()*psi.internalField()[curShape[0]]
 106               + curTetrahedron.Sb()*psi.internalField()[curShape[1]]
 107               + curTetrahedron.Sc()*psi.internalField()[curShape[2]]
 108               + curTetrahedron.Sd()*psi.internalField()[curShape[3]]
 109             )/curTetrahedron.mag();
 110
 111         forAll (curShape, pointI)
 112         {
 113             scalar weight =
 114                 curTetrahedron.mag()/
 115                 mag
 116                 (
 117                     points[curShape[pointI]]
 118                   - curShape.centre(points)
 119                 );
 120
 121             femGrad.internalField()[curShape[pointI]] += weight*tetGrad;
 122
 123             weights[curShape[pointI]] += weight;
 124         }
 125     }
 126
 127     femGrad.internalField() /= weights;
 128
 129     return tFemGrad;
 130 }
 131
 132
 133 template<class Type>
 134 tmp
 135 <
 136     GeometricField
 137     <
 138         typename outerProduct<vector, Type>::type,
 139         elementPatchField,
 140         elementMesh
 141     >
 142 >
 143 tetFec::elementGrad
 144 (
 145     const GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi
 146 )
 147 {
 148     typedef typename outerProduct<vector, Type>::type GradType;
 149
 150     const tetPolyMesh& tetMesh = psi.mesh();
 151
 152     const polyMesh& mesh = tetMesh();
 153
 154
 155     tmp<GeometricField<GradType, elementPatchField, elementMesh> > tElemGrad
 156     (
 157         new GeometricField<GradType, elementPatchField, elementMesh>
 158         (
 159             IOobject
 160             (
 161                 "grad("+psi.name()+')',
 162                 psi.instance(),
 163                 mesh,
 164                 IOobject::NO_READ,
 165                 IOobject::NO_WRITE
 166             ),
 167             tetMesh,
 168             dimensioned<GradType>
 169             (
 170                 "zero",
 171                 psi.dimensions()/dimLength,
 172                 pTraits<GradType>::zero
 173             )
 174         )
 175     );
 176
 177     GeometricField<GradType, elementPatchField, elementMesh>& elemGrad =
 178         tElemGrad();
 179
 180     pointField points = tetMesh.points();
 181
 182     scalarField weights (tetMesh.nCells(), 0.0);
 183
 184     const vectorField& C = mesh.cellCentres();
 185
 186
 187     for (label cellI = 0; cellI < tetMesh.nCells(); cellI++)
 188     {
 189         tetCellList tets = tetMesh.tets(cellI);
 190
 191         forAll (tets, tetI)
 192         {
 193             tetPointRef curTetrahedron = tets[tetI].tet(points);
 194
 195             cellShape curShape = tets[tetI].tetCellShape();
 196
 197             GradType tetGrad =
 198               - (1.0/3.0)*
 199                 (
 200                     curTetrahedron.Sa()*psi.internalField()[curShape[0]]
 201                   + curTetrahedron.Sb()*psi.internalField()[curShape[1]]
 202                   + curTetrahedron.Sc()*psi.internalField()[curShape[2]]
 203                   + curTetrahedron.Sd()*psi.internalField()[curShape[3]]
 204                 )/curTetrahedron.mag();
 205
 206             scalar weight = mag(C[cellI] - curShape.centre(points));
 207
 208             elemGrad.internalField()[cellI] += weight*tetGrad;
 209             weights[cellI] += weight;
 210         }
 211     }
 212
 213     elemGrad.internalField() /= weights;
 214
 215     return tElemGrad;
 216 }
 217
 218
 219
 220 template<class Type>
 221 tmp<GeometricField<Type, elementPatchField, elementMesh> > tetFec::ddt
 222 (
 223     const GeometricField<Type, elementPatchField, elementMesh>& ef
 224 )
 225 {
 226     const polyMesh& mesh = ef.mesh()();
 227
 228     dimensionedScalar rDeltaT = 1.0/mesh.time().deltaT();
 229
 230     IOobject ddtIOobject
 231     (
 232         "ddt("+ef.name()+')',
 233         mesh.time().timeName(),
 234         mesh,
 235         IOobject::NO_READ,
 236         IOobject::NO_WRITE
 237     );
 238
 239     return tmp<GeometricField<Type, elementPatchField, elementMesh> >
 240     (
 241         new GeometricField<Type, elementPatchField, elementMesh>
 242         (
 243             ddtIOobject,
 244             rDeltaT*(ef - ef.oldTime())
 245         )
 246     );
 247 }
 248
 249 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 250
 251 } // End namespace Foam
 252
 253 // ************************************************************************* //