Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "Kmesh.H"
00027 #include <OpenFOAM/polyMesh.H>
00028 #include <finiteVolume/volFields.H>
00029 #include <OpenFOAM/mathematicalConstants.H>
00030
00031
00032
00033 namespace Foam
00034 {
00036 inline label rep
00037 (
00038 const label i,
00039 const label j,
00040 const label k,
00041 const labelList& nn
00042 )
00043 {
00044 return (k + j*nn[2] + i*nn[1]*nn[2]);
00045 }
00047
00048 }
00049
00050
00051
00052
00053
00054 Foam::Kmesh::Kmesh(const fvMesh& mesh)
00055 :
00056 vectorField(mesh.V().size()),
00057 NN(vector::dim)
00058 {
00059 const scalar pi = mathematicalConstant::pi;
00060 const scalar twoPi = 2.0*pi;
00061
00062 boundBox box = mesh.bounds();
00063 L = box.span();
00064
00065 vector cornerCellCentre = ::Foam::max(mesh.C().internalField());
00066 vector cellL = 2 * (box.max() - cornerCellCentre);
00067
00068 vector rdeltaByL;
00069 label nTot = 1;
00070
00071 label i;
00072 forAll(NN, i)
00073 {
00074 NN[i] = label(L[i]/cellL[i] + 0.5);
00075 nTot *= NN[i];
00076
00077 if (NN[i] > 1)
00078 {
00079 L[i] -= cellL[i];
00080 }
00081
00082 rdeltaByL[i] = NN[i]/(L[i]*L[i]);
00083 }
00084
00085 if (nTot != mesh.nCells())
00086 {
00087 FatalErrorIn("Kmesh::Kmesh(const fvMesh& mesh)")
00088 << "calculated number of cells is incorrect"
00089 << abort(FatalError);
00090 }
00091
00092 for (i=0; i<NN[0]; i++)
00093 {
00094 scalar k1 = (i - NN[0]/2)*twoPi/L[0];
00095
00096 for (label j=0; j<NN[1]; j++)
00097 {
00098 scalar k2 = (j - NN[1]/2)*twoPi/L[1];
00099
00100 for (label k=0; k<NN[2]; k++)
00101 {
00102 scalar k3 = (k - NN[2]/2)*twoPi/L[2];
00103
00104 (*this)[rep(i, j, k, NN)] = vector(k1, k2, k3);
00105 }
00106 }
00107 }
00108
00109 Kmax = mag((*this)[rep(NN[0]-1, NN[1]-1, NN[2]-1, NN)]);
00110 }
00111
00112
00113