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-2023 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  PointField<Type>& pf
39 )
40 {
41  const polyMesh& mesh = pf.mesh();
42 
43  const polyBoundaryMesh& bm = mesh.boundaryMesh();
44 
45  // first count the total number of patch-patch points
46 
47  label nPatchPatchPoints = 0;
48 
49  forAll(bm, patchi)
50  {
51  if (!isA<emptyPolyPatch>(bm[patchi]))
52  {
53  nPatchPatchPoints += bm[patchi].boundaryPoints().size();
54  }
55  }
56 
57 
58  typename PointField<Type>::Boundary& bFld = pf.boundaryField();
59 
60 
61  // Evaluate in reverse order
62 
63  forAllReverse(bFld, patchi)
64  {
65  bFld[patchi].initEvaluate(Pstream::commsTypes::blocking); // buffered
66  }
67 
68  forAllReverse(bFld, patchi)
69  {
70  bFld[patchi].evaluate(Pstream::commsTypes::blocking);
71  }
72 
73 
74  // Save the values
75 
76  Field<Type> boundaryPointValues(nPatchPatchPoints);
77  nPatchPatchPoints = 0;
78 
79  forAll(bm, patchi)
80  {
81  if (!isA<emptyPolyPatch>(bm[patchi]))
82  {
83  const labelList& bp = bm[patchi].boundaryPoints();
84  const labelList& meshPoints = bm[patchi].meshPoints();
85 
86  forAll(bp, pointi)
87  {
88  label ppp = meshPoints[bp[pointi]];
89  boundaryPointValues[nPatchPatchPoints++] = pf[ppp];
90  }
91  }
92  }
93 
94 
95  // Forward evaluation
96 
97  bFld.evaluate();
98 
99 
100  // Check
101 
102  nPatchPatchPoints = 0;
103 
104  forAll(bm, patchi)
105  {
106  if (!isA<emptyPolyPatch>(bm[patchi]))
107  {
108  const labelList& bp = bm[patchi].boundaryPoints();
109  const labelList& meshPoints = bm[patchi].meshPoints();
110 
111  forAll(bp, pointi)
112  {
113  label ppp = meshPoints[bp[pointi]];
114 
115  const Type& savedVal = boundaryPointValues[nPatchPatchPoints++];
116 
117  if (savedVal != pf[ppp])
118  {
120  << "Patch fields are not consistent on mesh point "
121  << ppp << " coordinate " << mesh.points()[ppp]
122  << " at patch " << bm[patchi].name() << '.'
123  << endl
124  << "Reverse evaluation gives value " << savedVal
125  << " , forward evaluation gives value " << pf[ppp]
126  << abort(FatalError);
127  }
128  }
129  }
130  }
131 }
132 
133 
134 template<class Type>
136 Foam::motionSmootherAlgo::avg
137 (
138  const PointField<Type>& fld,
139  const scalarField& edgeWeight
140 ) const
141 {
143  (
145  (
146  "avg("+fld.name()+')',
147  fld.mesh(),
148  dimensioned<Type>("zero", fld.dimensions(), Zero)
149  )
150  );
151  PointField<Type>& res = tres.ref();
152 
153  const polyMesh& mesh = fld.mesh()();
154 
155 
156  // Sum local weighted values and weights
157  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
158 
159  // Note: on coupled edges use only one edge (through isMasterEdge)
160  // This is done so coupled edges do not get counted double.
161 
162  scalarField sumWeight(mesh.nPoints(), 0.0);
163 
164  const edgeList& edges = mesh.edges();
165 
166  forAll(edges, edgeI)
167  {
168  if (isMasterEdge_.get(edgeI) == 1)
169  {
170  const edge& e = edges[edgeI];
171  const scalar w = edgeWeight[edgeI];
172 
173  res[e[0]] += w*fld[e[1]];
174  sumWeight[e[0]] += w;
175 
176  res[e[1]] += w*fld[e[0]];
177  sumWeight[e[1]] += w;
178  }
179  }
180 
181 
182  // Add coupled contributions
183  // ~~~~~~~~~~~~~~~~~~~~~~~~~
184 
186  (
187  mesh,
188  res,
189  plusEqOp<Type>(),
190  Type(Zero) // null value
191  );
193  (
194  mesh,
195  sumWeight,
197  scalar(0) // null value
198  );
199 
200 
201  // Average
202  // ~~~~~~~
203 
204  forAll(res, pointi)
205  {
206  if (mag(sumWeight[pointi]) < vSmall)
207  {
208  // Unconnected point. Take over original value
209  res[pointi] = fld[pointi];
210  }
211  else
212  {
213  res[pointi] /= sumWeight[pointi];
214  }
215  }
216 
217  // Single and multi-patch constraints
218  pointConstraints::New(fld.mesh()).constrain(res, false);
219 
220  return tres;
221 }
222 
223 
224 template<class Type>
226 (
227  const PointField<Type>& fld,
228  const scalarField& edgeWeight,
229  PointField<Type>& newFld
230 ) const
231 {
232  tmp<pointVectorField> tavgFld = avg(fld, edgeWeight);
233  const pointVectorField& avgFld = tavgFld();
234 
235  forAll(fld, pointi)
236  {
237  if (isInternalPoint(pointi))
238  {
239  newFld[pointi] = 0.5*fld[pointi] + 0.5*avgFld[pointi];
240  }
241  }
242 
243  // Single and multi-patch constraints
244  pointConstraints::New(fld.mesh()).constrain(newFld, false);
245 }
246 
247 
248 template<class Type, class CombineOp>
249 void Foam::motionSmootherAlgo::testSyncField
250 (
251  const Field<Type>& fld,
252  const CombineOp& cop,
253  const Type& zero,
254  const scalar maxMag
255 ) const
256 {
257  if (debug)
258  {
259  Pout<< "testSyncField : testing synchronisation of Field<Type>."
260  << endl;
261  }
262 
263  Field<Type> syncedFld(fld);
264 
266  (
267  mesh_,
268  syncedFld,
269  cop,
270  zero
271  );
272 
273  forAll(syncedFld, i)
274  {
275  if (mag(syncedFld[i] - fld[i]) > maxMag)
276  {
278  << "On element " << i << " value:" << fld[i]
279  << " synchronised value:" << syncedFld[i]
280  << abort(FatalError);
281  }
282  }
283 }
284 
285 
286 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
static pointConstraints & New(const word &name, const pointMesh &mesh)
Construct and return the named DemandDrivenMeshObject.
GeometricBoundaryField< Type, pointPatchField, pointMesh > Boundary
Type of the boundary field.
Generic dimensioned Type class.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:61
const polyMesh & mesh() const
Reference to mesh.
void smooth(const PointField< Type > &fld, const scalarField &edgeWeight, PointField< Type > &newFld) const
Fully explicit smoothing of fields (not positions)
void constrain(PointField< Type > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
label nPoints() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
A class for managing temporary objects.
Definition: tmp.H:55
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:50
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
gmvFile<< "tracers "<< particles.size()<< nl;{ pointField positions(particles.size());label particlei=0;forAllConstIter(Cloud< passiveParticle >, particles, iter) { positions[particlei++]=iter().position(mesh);} for(i=0;i< pTraits< point >::nComponents;i++) { forAll(positions, particlei) { gmvFile<< component(positions[particlei], i)<< ' ';} gmvFile<< nl;}}forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static const zero Zero
Definition: zero.H:97
const doubleScalar e
Definition: doubleScalar.H:106
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensioned< scalar > mag(const dimensioned< Type > &)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
error FatalError