FreeFOAM The Cross-Platform CFD Toolkit
motionSmootherTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "motionSmoother.H"
27 #include <meshTools/meshTools.H>
31 #include <OpenFOAM/syncTools.H>
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 template<class Type>
36 void Foam::motionSmoother::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::GeometricBoundaryField& bFld = pf.boundaryField();
61 
62 
63  // Evaluate in reverse order
64 
65  forAllReverse(bFld, patchi)
66  {
67  bFld[patchi].initEvaluate(Pstream::blocking); // buffered
68  }
69 
70  forAllReverse(bFld, patchi)
71  {
72  bFld[patchi].evaluate(Pstream::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  (
123  "motionSmoother::checkConstraints"
124  "(GeometricField<Type, pointPatchField, pointMesh>&)"
125  ) << "Patch fields are not consistent on mesh point "
126  << ppp << " coordinate " << mesh.points()[ppp]
127  << " at patch " << bm[patchi].name() << '.'
128  << endl
129  << "Reverse evaluation gives value " << savedVal
130  << " , forward evaluation gives value " << pf[ppp]
131  << abort(FatalError);
132  }
133  }
134  }
135  }
136 }
137 
138 
139 template<class Type>
140 void Foam::motionSmoother::applyCornerConstraints
141 (
142  GeometricField<Type, pointPatchField, pointMesh>& pf
143 ) const
144 {
145  forAll(patchPatchPointConstraintPoints_, pointi)
146  {
147  pf[patchPatchPointConstraintPoints_[pointi]] = transform
148  (
149  patchPatchPointConstraintTensors_[pointi],
150  pf[patchPatchPointConstraintPoints_[pointi]]
151  );
152  }
153 }
154 
155 
156 // Average of connected points.
157 template <class Type>
159  Foam::motionSmoother::avg
160 (
161  const GeometricField<Type, pointPatchField, pointMesh>& fld,
162  const scalarField& edgeWeight,
163  const bool separation
164 ) const
165 {
166  tmp<GeometricField<Type, pointPatchField, pointMesh> > tres
167  (
168  new GeometricField<Type, pointPatchField, pointMesh>
169  (
170  IOobject
171  (
172  "avg("+fld.name()+')',
173  fld.time().timeName(),
174  fld.db(),
177  ),
178  fld.mesh(),
179  dimensioned<Type>("zero", fld.dimensions(), pTraits<Type>::zero)
180  )
181  );
182  GeometricField<Type, pointPatchField, pointMesh>& res = tres();
183 
184  const polyMesh& mesh = fld.mesh()();
185 
186 
187  // Sum local weighted values and weights
188  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
189 
190  // Note: on coupled edges use only one edge (through isMasterEdge)
191  // This is done so coupled edges do not get counted double.
192 
193  scalarField sumWeight(mesh.nPoints(), 0.0);
194 
195  const edgeList& edges = mesh.edges();
196 
197  forAll(edges, edgeI)
198  {
199  if (isMasterEdge_.get(edgeI) == 1)
200  {
201  const edge& e = edges[edgeI];
202  const scalar w = edgeWeight[edgeI];
203 
204  res[e[0]] += w*fld[e[1]];
205  sumWeight[e[0]] += w;
206 
207  res[e[1]] += w*fld[e[0]];
208  sumWeight[e[1]] += w;
209  }
210  }
211 
212 
213  // Add coupled contributions
214  // ~~~~~~~~~~~~~~~~~~~~~~~~~
215 
217  (
218  mesh,
219  res,
220  plusEqOp<Type>(),
221  pTraits<Type>::zero, // null value
222  separation // separation
223  );
224 
226  (
227  mesh,
228  sumWeight,
229  plusEqOp<scalar>(),
230  scalar(0), // null value
231  false // separation
232  );
233 
234 
235  // Average
236  // ~~~~~~~
237 
238  forAll(res, pointI)
239  {
240  if (mag(sumWeight[pointI]) < VSMALL)
241  {
242  // Unconnected point. Take over original value
243  res[pointI] = fld[pointI];
244  }
245  else
246  {
247  res[pointI] /= sumWeight[pointI];
248  }
249  }
250 
251  res.correctBoundaryConditions();
252  applyCornerConstraints(res);
253 
254  return tres;
255 }
256 
257 
258 // smooth field (point-jacobi)
259 template <class Type>
261 (
263  const scalarField& edgeWeight,
264  const bool separation,
266 ) const
267 {
268  tmp<pointVectorField> tavgFld = avg
269  (
270  fld,
271  edgeWeight, // weighting
272  separation // whether to apply separation vector
273  );
274  const pointVectorField& avgFld = tavgFld();
275 
276  forAll(fld, pointI)
277  {
278  if (isInternalPoint(pointI))
279  {
280  newFld[pointI] = 0.5*fld[pointI] + 0.5*avgFld[pointI];
281  }
282  }
283 
284  newFld.correctBoundaryConditions();
285  applyCornerConstraints(newFld);
286 }
287 
288 
289 //- Test synchronisation of pointField
290 template<class Type, class CombineOp>
291 void Foam::motionSmoother::testSyncField
292 (
293  const Field<Type>& fld,
294  const CombineOp& cop,
295  const Type& zero,
296  const bool separation,
297  const scalar maxMag
298 ) const
299 {
300  if (debug)
301  {
302  Pout<< "testSyncField : testing synchronisation of Field<Type>."
303  << endl;
304  }
305 
306  Field<Type> syncedFld(fld);
307 
309  (
310  mesh_,
311  syncedFld,
312  cop,
313  zero, // null value
314  separation // separation
315  );
316 
317  forAll(syncedFld, i)
318  {
319  if (mag(syncedFld[i] - fld[i]) > maxMag)
320  {
322  (
323  "motionSmoother::testSyncField"
324  "(const Field<Type>&, const CombineOp&"
325  ", const Type&, const bool)"
326  ) << "On element " << i << " value:" << fld[i]
327  << " synchronised value:" << syncedFld[i]
328  << abort(FatalError);
329  }
330  }
331 }
332 
333 
334 // ************************ vim: set sw=4 sts=4 et: ************************ //