src/randomProcesses/Kmesh/Kmesh.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 "Kmesh.H"
  27 #include "polyMesh.H"
  28 #include "volFields.H"
  29 #include "mathematicalConstants.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35    //! \cond fileScope
  36    inline label rep
  37    (
  38        const label i,
  39        const label j,
  40        const label k,
  41        const labelList& nn
  42    )
  43    {
  44        return (k + j*nn[2] + i*nn[1]*nn[2]);
  45    }
  46    //! \endcond
  47
  48 } // End namespace Foam
  49
  50
  51 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  52
  53 // from fvMesh
  54 Foam::Kmesh::Kmesh(const fvMesh& mesh)
  55 :
  56     vectorField(mesh.V().size()),
  57     NN(vector::dim)
  58 {
  59     boundBox box = mesh.bounds();
  60     L = box.span();
  61
  62     vector cornerCellCentre = ::Foam::max(mesh.C().internalField());
  63     vector cellL = 2*(box.max() - cornerCellCentre);
  64
  65     vector rdeltaByL;
  66     label nTot = 1;
  67
  68     label i;
  69     forAll(NN, i)
  70     {
  71         NN[i] = label(L[i]/cellL[i] + 0.5);
  72         nTot *= NN[i];
  73
  74         if (NN[i] > 1)
  75         {
  76             L[i] -= cellL[i];
  77         }
  78
  79         rdeltaByL[i] = NN[i]/(L[i]*L[i]);
  80     }
  81
  82     if (nTot != mesh.nCells())
  83     {
  84         FatalErrorIn("Kmesh::Kmesh(const fvMesh& mesh)")
  85             << "calculated number of cells is incorrect"
  86             << abort(FatalError);
  87     }
  88
  89     for (i=0; i<NN[0]; i++)
  90     {
  91         scalar k1 = (i - NN[0]/2)*constant::mathematical::twoPi/L[0];
  92
  93         for (label j=0; j<NN[1]; j++)
  94         {
  95             scalar k2 = (j - NN[1]/2)*constant::mathematical::twoPi/L[1];
  96
  97             for (label k=0; k<NN[2]; k++)
  98             {
  99                 scalar k3 = (k - NN[2]/2)*constant::mathematical::twoPi/L[2];
 100
 101                 (*this)[rep(i, j, k, NN)] = vector(k1, k2, k3);
 102             }
 103         }
 104     }
 105
 106     Kmax = mag((*this)[rep(NN[0]-1, NN[1]-1, NN[2]-1, NN)]);
 107 }
 108
 109
 110 // ************************************************************************* //