motionSmootherAlgoTemplates.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) 2011-2019 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 "motionSmootherAlgo.H"
27 #include "meshTools.H"
29 #include "pointConstraint.H"
30 #include "pointConstraints.H"
31 #include "syncTools.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 template<class Type>
36 void Foam::motionSmootherAlgo::checkConstraints
37 (
38  GeometricField<Type, pointPatchField, pointMesh>& pf
39 )
40 {
41  typedef GeometricField<Type, pointPatchField, pointMesh> FldType;
42 
43  const polyMesh& mesh = pf.mesh();
44 
45  const polyBoundaryMesh& bm = mesh.boundaryMesh();
46 
47  // first count the total number of patch-patch points
48 
49  label nPatchPatchPoints = 0;
50 
51  forAll(bm, patchi)
52  {
53  if (!isA<emptyPolyPatch>(bm[patchi]))
54  {
55  nPatchPatchPoints += bm[patchi].boundaryPoints().size();
56  }
57  }
58 
59 
60  typename FldType::Boundary& bFld = pf.boundaryField();
61 
62 
63  // Evaluate in reverse order
64 
65  forAllReverse(bFld, patchi)
66  {
67  bFld[patchi].initEvaluate(Pstream::commsTypes::blocking); // buffered
68  }
69 
70  forAllReverse(bFld, patchi)
71  {
72  bFld[patchi].evaluate(Pstream::commsTypes::blocking);
73  }
74 
75 
76  // Save the values
77 
78  Field<Type> boundaryPointValues(nPatchPatchPoints);
79  nPatchPatchPoints = 0;
80 
81  forAll(bm, patchi)
82  {
83  if (!isA<emptyPolyPatch>(bm[patchi]))
84  {
85  const labelList& bp = bm[patchi].boundaryPoints();
86  const labelList& meshPoints = bm[patchi].meshPoints();
87 
88  forAll(bp, pointi)
89  {
90  label ppp = meshPoints[bp[pointi]];
91  boundaryPointValues[nPatchPatchPoints++] = pf[ppp];
92  }
93  }
94  }
95 
96 
97  // Forward evaluation
98 
99  bFld.evaluate();
100 
101 
102  // Check
103 
104  nPatchPatchPoints = 0;
105 
106  forAll(bm, patchi)
107  {
108  if (!isA<emptyPolyPatch>(bm[patchi]))
109  {
110  const labelList& bp = bm[patchi].boundaryPoints();
111  const labelList& meshPoints = bm[patchi].meshPoints();
112 
113  forAll(bp, pointi)
114  {
115  label ppp = meshPoints[bp[pointi]];
116 
117  const Type& savedVal = boundaryPointValues[nPatchPatchPoints++];
118 
119  if (savedVal != pf[ppp])
120  {
122  << "Patch fields are not consistent on mesh point "
123  << ppp << " coordinate " << mesh.points()[ppp]
124  << " at patch " << bm[patchi].name() << '.'
125  << endl
126  << "Reverse evaluation gives value " << savedVal
127  << " , forward evaluation gives value " << pf[ppp]
128  << abort(FatalError);
129  }
130  }
131  }
132  }
133 }
134 
135 
136 template<class Type>
138 Foam::motionSmootherAlgo::avg
139 (
141  const scalarField& edgeWeight
142 ) const
143 {
145  (
147  (
148  "avg("+fld.name()+')',
149  fld.mesh(),
150  dimensioned<Type>("zero", fld.dimensions(), Zero)
151  )
152  );
154 
155  const polyMesh& mesh = fld.mesh()();
156 
157 
158  // Sum local weighted values and weights
159  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
160 
161  // Note: on coupled edges use only one edge (through isMasterEdge)
162  // This is done so coupled edges do not get counted double.
163 
164  scalarField sumWeight(mesh.nPoints(), 0.0);
165 
166  const edgeList& edges = mesh.edges();
167 
168  forAll(edges, edgeI)
169  {
170  if (isMasterEdge_.get(edgeI) == 1)
171  {
172  const edge& e = edges[edgeI];
173  const scalar w = edgeWeight[edgeI];
174 
175  res[e[0]] += w*fld[e[1]];
176  sumWeight[e[0]] += w;
177 
178  res[e[1]] += w*fld[e[0]];
179  sumWeight[e[1]] += w;
180  }
181  }
182 
183 
184  // Add coupled contributions
185  // ~~~~~~~~~~~~~~~~~~~~~~~~~
186 
188  (
189  mesh,
190  res,
191  plusEqOp<Type>(),
192  Type(Zero) // null value
193  );
195  (
196  mesh,
197  sumWeight,
199  scalar(0) // null value
200  );
201 
202 
203  // Average
204  // ~~~~~~~
205 
206  forAll(res, pointi)
207  {
208  if (mag(sumWeight[pointi]) < vSmall)
209  {
210  // Unconnected point. Take over original value
211  res[pointi] = fld[pointi];
212  }
213  else
214  {
215  res[pointi] /= sumWeight[pointi];
216  }
217  }
218 
219  // Single and multi-patch constraints
220  pointConstraints::New(fld.mesh()).constrain(res, false);
221 
222  return tres;
223 }
224 
225 
226 template<class Type>
228 (
230  const scalarField& edgeWeight,
232 ) const
233 {
234  tmp<pointVectorField> tavgFld = avg(fld, edgeWeight);
235  const pointVectorField& avgFld = tavgFld();
236 
237  forAll(fld, pointi)
238  {
239  if (isInternalPoint(pointi))
240  {
241  newFld[pointi] = 0.5*fld[pointi] + 0.5*avgFld[pointi];
242  }
243  }
244 
245  // Single and multi-patch constraints
246  pointConstraints::New(fld.mesh()).constrain(newFld, false);
247 }
248 
249 
250 template<class Type, class CombineOp>
251 void Foam::motionSmootherAlgo::testSyncField
252 (
253  const Field<Type>& fld,
254  const CombineOp& cop,
255  const Type& zero,
256  const scalar maxMag
257 ) const
258 {
259  if (debug)
260  {
261  Pout<< "testSyncField : testing synchronisation of Field<Type>."
262  << endl;
263  }
264 
265  Field<Type> syncedFld(fld);
266 
268  (
269  mesh_,
270  syncedFld,
271  cop,
272  zero
273  );
274 
275  forAll(syncedFld, i)
276  {
277  if (mag(syncedFld[i] - fld[i]) > maxMag)
278  {
280  << "On element " << i << " value:" << fld[i]
281  << " synchronised value:" << syncedFld[i]
282  << abort(FatalError);
283  }
284  }
285 }
286 
287 
288 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
Definition: IOobject.H:315
unsigned int get(const label) const
Get value at index I.
Definition: PackedListI.H:954
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
Generic dimensioned Type class.
const dimensionSet & dimensions() const
Return dimensions.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void smooth(const GeometricField< Type, pointPatchField, pointMesh > &fld, const scalarField &edgeWeight, GeometricField< Type, pointPatchField, pointMesh > &newFld) const
Fully explicit smoothing of fields (not positions)
const Mesh & mesh() const
Return mesh.
Internal & ref()
Return a reference to the dimensioned internal field.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
fvConstraints constrain(rhoEqn)
label patchi
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A class for managing temporary objects.
Definition: PtrList.H:53