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-2018 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  IOobject
149  (
150  "avg("+fld.name()+')',
151  fld.time().timeName(),
152  fld.db(),
155  false
156  ),
157  fld.mesh(),
158  dimensioned<Type>("zero", fld.dimensions(), Zero)
159  )
160  );
162 
163  const polyMesh& mesh = fld.mesh()();
164 
165 
166  // Sum local weighted values and weights
167  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
168 
169  // Note: on coupled edges use only one edge (through isMasterEdge)
170  // This is done so coupled edges do not get counted double.
171 
172  scalarField sumWeight(mesh.nPoints(), 0.0);
173 
174  const edgeList& edges = mesh.edges();
175 
176  forAll(edges, edgeI)
177  {
178  if (isMasterEdge_.get(edgeI) == 1)
179  {
180  const edge& e = edges[edgeI];
181  const scalar w = edgeWeight[edgeI];
182 
183  res[e[0]] += w*fld[e[1]];
184  sumWeight[e[0]] += w;
185 
186  res[e[1]] += w*fld[e[0]];
187  sumWeight[e[1]] += w;
188  }
189  }
190 
191 
192  // Add coupled contributions
193  // ~~~~~~~~~~~~~~~~~~~~~~~~~
194 
196  (
197  mesh,
198  res,
199  plusEqOp<Type>(),
200  Type(Zero) // null value
201  );
203  (
204  mesh,
205  sumWeight,
207  scalar(0) // null value
208  );
209 
210 
211  // Average
212  // ~~~~~~~
213 
214  forAll(res, pointi)
215  {
216  if (mag(sumWeight[pointi]) < vSmall)
217  {
218  // Unconnected point. Take over original value
219  res[pointi] = fld[pointi];
220  }
221  else
222  {
223  res[pointi] /= sumWeight[pointi];
224  }
225  }
226 
227  // Single and multi-patch constraints
228  pointConstraints::New(fld.mesh()).constrain(res, false);
229 
230  return tres;
231 }
232 
233 
234 template<class Type>
236 (
238  const scalarField& edgeWeight,
240 ) const
241 {
242  tmp<pointVectorField> tavgFld = avg(fld, edgeWeight);
243  const pointVectorField& avgFld = tavgFld();
244 
245  forAll(fld, pointi)
246  {
247  if (isInternalPoint(pointi))
248  {
249  newFld[pointi] = 0.5*fld[pointi] + 0.5*avgFld[pointi];
250  }
251  }
252 
253  // Single and multi-patch constraints
254  pointConstraints::New(fld.mesh()).constrain(newFld, false);
255 }
256 
257 
258 template<class Type, class CombineOp>
259 void Foam::motionSmootherAlgo::testSyncField
260 (
261  const Field<Type>& fld,
262  const CombineOp& cop,
263  const Type& zero,
264  const scalar maxMag
265 ) const
266 {
267  if (debug)
268  {
269  Pout<< "testSyncField : testing synchronisation of Field<Type>."
270  << endl;
271  }
272 
273  Field<Type> syncedFld(fld);
274 
276  (
277  mesh_,
278  syncedFld,
279  cop,
280  zero
281  );
282 
283  forAll(syncedFld, i)
284  {
285  if (mag(syncedFld[i] - fld[i]) > maxMag)
286  {
288  << "On element " << i << " value:" << fld[i]
289  << " synchronised value:" << syncedFld[i]
290  << abort(FatalError);
291  }
292  }
293 }
294 
295 
296 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
Definition: IOobject.H:295
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:319
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
Generic GeometricField class.
Generic dimensioned Type class.
static const pointConstraints & New(const pointMesh &mesh)
Definition: MeshObject.C:44
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.
fvOptions constrain(rhoEqn)
label patchi
const Time & time() const
Return time.
Definition: IOobject.C:360
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A class for managing temporary objects.
Definition: PtrList.H:53
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:354
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92