src/tetDecompositionFiniteElement/tetFemCalculus/tetFec.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Class of static functions to calculate explicit finite element derivatives.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "tetCell.H"
  31 #include "tetCellList.H"
  32 #include "tetPointRef.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  40
  41
  42 template<class Type>
  43 tmp
  44 <
  45     GeometricField
  46     <
  47         typename outerProduct<vector, Type>::type,
  48         tetPolyPatchField,
  49         tetPointMesh
  50     >
  51 >
  52 tetFec::grad
  53 (
  54     const GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi
  55 )
  56 {
  57     typedef typename outerProduct<vector, Type>::type GradType;
  58
  59     const tetPolyMesh& tetMesh = psi.mesh();
  60
  61     tmp<GeometricField<GradType, tetPolyPatchField, tetPointMesh> > tFemGrad
  62     (
  63         new GeometricField<GradType, tetPolyPatchField, tetPointMesh>
  64         (
  65             IOobject
  66             (
  67                 "grad("+psi.name()+')',
  68                 psi.instance(),
  69                 tetMesh(),
  70                 IOobject::NO_READ,
  71                 IOobject::NO_WRITE
  72             ),
  73             tetMesh,
  74             dimensioned<GradType>
  75             (
  76                 "zero",
  77                 psi.dimensions()/dimLength,
  78                 pTraits<GradType>::zero
  79             )
  80         )
  81     );
  82
  83     GeometricField<GradType, tetPolyPatchField, tetPointMesh>& femGrad =
  84         tFemGrad();
  85
  86     pointField points = tetMesh.points();
  87     cellShapeList tetCellShapes = tetMesh.tetCells();
  88
  89     scalarField weights (points.size(), 0.0);
  90
  91     forAll (tetCellShapes, tetI)
  92     {
  93         cellShape& curShape = tetCellShapes[tetI];
  94
  95         tetPointRef curTetrahedron
  96         (
  97             points[curShape[0]],
  98             points[curShape[1]],
  99             points[curShape[2]],
 100             points[curShape[3]]
 101         );
 102
 103         GradType tetGrad =
 104           - (1.0/3.0)*
 105             (
 106                 curTetrahedron.Sa()*psi.internalField()[curShape[0]]
 107               + curTetrahedron.Sb()*psi.internalField()[curShape[1]]
 108               + curTetrahedron.Sc()*psi.internalField()[curShape[2]]
 109               + curTetrahedron.Sd()*psi.internalField()[curShape[3]]
 110             )/curTetrahedron.mag();
 111
 112         forAll (curShape, pointI)
 113         {
 114             scalar weight =
 115                 curTetrahedron.mag()/
 116                 mag
 117                 (
 118                     points[curShape[pointI]]
 119                   - curShape.centre(points)
 120                 );
 121
 122             femGrad.internalField()[curShape[pointI]] += weight*tetGrad;
 123
 124             weights[curShape[pointI]] += weight;
 125         }
 126     }
 127
 128     femGrad.internalField() /= weights;
 129
 130     return tFemGrad;
 131 }
 132
 133
 134 template<class Type>
 135 tmp
 136 <
 137     GeometricField
 138     <
 139         typename outerProduct<vector, Type>::type,
 140         elementPatchField,
 141         elementMesh
 142     >
 143 >
 144 tetFec::elementGrad
 145 (
 146     const GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi
 147 )
 148 {
 149     typedef typename outerProduct<vector, Type>::type GradType;
 150
 151     const tetPolyMesh& tetMesh = psi.mesh();
 152
 153     const polyMesh& mesh = tetMesh();
 154
 155
 156     tmp<GeometricField<GradType, elementPatchField, elementMesh> > tElemGrad
 157     (
 158         new GeometricField<GradType, elementPatchField, elementMesh>
 159         (
 160             IOobject
 161             (
 162                 "grad("+psi.name()+')',
 163                 psi.instance(),
 164                 mesh,
 165                 IOobject::NO_READ,
 166                 IOobject::NO_WRITE
 167             ),
 168             tetMesh,
 169             dimensioned<GradType>
 170             (
 171                 "zero",
 172                 psi.dimensions()/dimLength,
 173                 pTraits<GradType>::zero
 174             )
 175         )
 176     );
 177
 178     GeometricField<GradType, elementPatchField, elementMesh>& elemGrad =
 179         tElemGrad();
 180
 181     pointField points = tetMesh.points();
 182
 183     scalarField weights (tetMesh.nCells(), 0.0);
 184
 185     const vectorField& C = mesh.cellCentres();
 186
 187
 188     for (label cellI = 0; cellI < tetMesh.nCells(); cellI++)
 189     {
 190         tetCellList tets = tetMesh.tets(cellI);
 191
 192         forAll (tets, tetI)
 193         {
 194             tetPointRef curTetrahedron = tets[tetI].tet(points);
 195
 196             cellShape curShape = tets[tetI].tetCellShape();
 197
 198             GradType tetGrad =
 199               - (1.0/3.0)*
 200                 (
 201                     curTetrahedron.Sa()*psi.internalField()[curShape[0]]
 202                   + curTetrahedron.Sb()*psi.internalField()[curShape[1]]
 203                   + curTetrahedron.Sc()*psi.internalField()[curShape[2]]
 204                   + curTetrahedron.Sd()*psi.internalField()[curShape[3]]
 205                 )/curTetrahedron.mag();
 206
 207             scalar weight = mag(C[cellI] - curShape.centre(points));
 208
 209             elemGrad.internalField()[cellI] += weight*tetGrad;
 210             weights[cellI] += weight;
 211         }
 212     }
 213
 214     elemGrad.internalField() /= weights;
 215
 216     return tElemGrad;
 217 }
 218
 219
 220
 221 template<class Type>
 222 tmp<GeometricField<Type, elementPatchField, elementMesh> > tetFec::ddt
 223 (
 224     const GeometricField<Type, elementPatchField, elementMesh>& ef
 225 )
 226 {
 227     const polyMesh& mesh = ef.mesh()();
 228
 229     dimensionedScalar rDeltaT = 1.0/mesh.time().deltaT();
 230
 231     IOobject ddtIOobject
 232     (
 233         "ddt("+ef.name()+')',
 234         mesh.time().timeName(),
 235         mesh,
 236         IOobject::NO_READ,
 237         IOobject::NO_WRITE
 238     );
 239
 240     return tmp<GeometricField<Type, elementPatchField, elementMesh> >
 241     (
 242         new GeometricField<Type, elementPatchField, elementMesh>
 243         (
 244             ddtIOobject,
 245             rDeltaT*(ef - ef.oldTime())
 246         )
 247     );
 248 }
 249
 250 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 251
 252 } // End namespace Foam
 253
 254 // ************************************************************************* //