Moment.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2020 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 "Moment.H"
27 
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 template<class Type>
33 (
34  const IOobject& io,
35  const dictionary& dict,
36  const fvMesh& mesh
37 )
38 :
39  AveragingMethod<Type>(io, dict, mesh, labelList(4, mesh.nCells())),
44  transform_(mesh.nCells(), Zero),
45  scale_(0.5*pow(mesh.V(), 1.0/3.0))
46 {
47  scalar a = 1.0/24.0;
48  scalar b = 0.5854101966249685;
49  scalar c = 0.1381966011250105;
50 
51  scalarField wQ(4);
52  wQ[0] = a;
53  wQ[1] = a;
54  wQ[2] = a;
55  wQ[3] = a;
56 
57  vectorField xQ(4);
58  xQ[0] = vector(b, c, c);
59  xQ[1] = vector(c, b, c);
60  xQ[2] = vector(c, c, b);
61  xQ[3] = vector(c, c, c);
62 
63  forAll(mesh.C(), celli)
64  {
65  const List<tetIndices> cellTets =
66  polyMeshTetDecomposition::cellTetIndices(mesh, celli);
67 
68  symmTensor A(Zero);
69 
70  forAll(cellTets, tetI)
71  {
72  const tetIndices& tetIs = cellTets[tetI];
73  const triFace triIs = tetIs.faceTriIs(mesh);
74 
75  const tensor T
76  (
77  tensor
78  (
79  mesh.points()[triIs[0]] - mesh.C()[celli],
80  mesh.points()[triIs[1]] - mesh.C()[celli],
81  mesh.points()[triIs[2]] - mesh.C()[celli]
82  ).T()
83  );
84 
85  const vectorField d((T & xQ)/scale_[celli]);
86 
87  const scalar v(6.0*tetIs.tet(mesh).mag()/mesh.V()[celli]);
88 
89  A += v*sum(wQ*sqr(d));
90  }
91 
92  transform_[celli] = inv(A);
93  }
94 }
95 
96 
97 template<class Type>
99 (
100  const Moment<Type>& am
101 )
102 :
108  transform_(am.transform_)
109 {}
110 
111 
112 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
113 
114 template<class Type>
116 {}
117 
118 
119 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
120 
121 template<class Type>
123 {}
124 
125 
126 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
127 
128 template<class Type>
130 (
131  const barycentric& coordinates,
132  const tetIndices& tetIs,
133  const Type& value
134 )
135 {
136  const label celli = tetIs.cell();
137  const triFace triIs = tetIs.faceTriIs(this->mesh_);
138 
139  const point delta =
140  (coordinates[0] - 1)*this->mesh_.C()[celli]
141  + coordinates[1]*this->mesh_.points()[triIs[0]]
142  + coordinates[2]*this->mesh_.points()[triIs[1]]
143  + coordinates[3]*this->mesh_.points()[triIs[2]];
144 
145  const Type v = value/this->mesh_.V()[celli];
146  const TypeGrad dv = transform_[celli] & (v*delta/scale_[celli]);
147 
148  data_[celli] += v;
149  dataX_[celli] += v + dv.x();
150  dataY_[celli] += v + dv.y();
151  dataZ_[celli] += v + dv.z();
152 }
153 
154 
155 template<class Type>
157 (
158  const barycentric& coordinates,
159  const tetIndices& tetIs
160 ) const
161 {
162  const label celli = tetIs.cell();
163  const triFace triIs = tetIs.faceTriIs(this->mesh_);
164 
165  const point delta =
166  (coordinates[0] - 1)*this->mesh_.C()[celli]
167  + coordinates[1]*this->mesh_.points()[triIs[0]]
168  + coordinates[2]*this->mesh_.points()[triIs[1]]
169  + coordinates[3]*this->mesh_.points()[triIs[2]];
170 
171  return
172  data_[celli]
173  + (
174  TypeGrad
175  (
176  dataX_[celli] - data_[celli],
177  dataY_[celli] - data_[celli],
178  dataZ_[celli] - data_[celli]
179  )
180  & delta/scale_[celli]
181  );
182 }
183 
184 
185 template<class Type>
188 (
189  const barycentric& coordinates,
190  const tetIndices& tetIs
191 ) const
192 {
193  const label celli(tetIs.cell());
194 
195  return
196  TypeGrad
197  (
198  dataX_[celli] - data_[celli],
199  dataY_[celli] - data_[celli],
200  dataZ_[celli] - data_[celli]
201  )/scale_[celli];
202 }
203 
204 
205 template<class Type>
208 {
209  return tmp<Field<Type>>(data_);
210 }
211 
212 
213 // ************************************************************************* //
scalar delta
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
AveragingMethod< Type >::TypeGrad TypeGrad
Gradient type.
Definition: Moment.H:70
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
void add(const barycentric &coordinates, const tetIndices &tetIs, const Type &value)
Member Functions.
Definition: Moment.C:130
Type interpolate(const barycentric &coordinates, const tetIndices &tetIs) const
Interpolate.
Definition: Moment.C:157
outerProduct< vector, Type >::type TypeGrad
Protected typedefs.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
triFace faceTriIs(const polyMesh &mesh) const
Return the indices corresponding to the tri on the face for.
Definition: tetIndicesI.H:67
fvMesh & mesh
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
Definition: tetIndicesI.H:108
const dimensionedScalar c
Speed of light in a vacuum.
Generic field type.
Definition: FieldField.H:51
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Moment lagrangian averaging procedure.
Definition: Moment.H:61
TypeGrad interpolateGrad(const barycentric &coordinates, const tetIndices &tetIs) const
Interpolate gradient.
Definition: Moment.C:188
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
label cell() const
Return the cell.
Definition: tetIndicesI.H:31
scalar mag() const
Return volume.
Definition: tetrahedronI.H:170
Base class for lagrangian averaging methods.
Definition: MPPICParcel.H:55
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
tmp< Field< Type > > primitiveField() const
Return an internal field of the average.
Definition: Moment.C:207
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const volScalarField & T
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
Moment(const IOobject &io, const dictionary &dict, const fvMesh &mesh)
Constructors.
Definition: Moment.C:33
virtual ~Moment()
Destructor.
Definition: Moment.C:115
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98