volPointInterpolate.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-2022 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 "volPointInterpolation.H"
27 #include "volFields.H"
28 #include "pointFields.H"
29 #include "emptyFvPatch.H"
30 #include "coupledPointPatchField.H"
31 #include "pointConstraints.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 template<class Type>
36 void Foam::volPointInterpolation::pushUntransformedData
37 (
38  List<Type>& pointData
39 ) const
40 {
41  // Transfer onto coupled patch
42  const globalMeshData& gmd = mesh().globalData();
43  const indirectPrimitivePatch& cpp = gmd.coupledPatch();
44  const labelList& meshPoints = cpp.meshPoints();
45 
46  const distributionMap& slavesMap = gmd.globalCoPointSlavesMap();
47  const labelListList& slaves = gmd.globalCoPointSlaves();
48 
49  List<Type> elems(slavesMap.constructSize());
50  forAll(meshPoints, i)
51  {
52  elems[i] = pointData[meshPoints[i]];
53  }
54 
55  // Combine master data with slave data
56  forAll(slaves, i)
57  {
58  const labelList& slavePoints = slaves[i];
59 
60  // Copy master data to slave slots
61  forAll(slavePoints, j)
62  {
63  elems[slavePoints[j]] = elems[i];
64  }
65  }
66 
67  // Push slave-slot data back to slaves
68  slavesMap.reverseDistribute(elems.size(), elems, false);
69 
70  // Extract back onto mesh
71  forAll(meshPoints, i)
72  {
73  pointData[meshPoints[i]] = elems[i];
74  }
75 }
76 
77 
78 template<class Type>
79 void Foam::volPointInterpolation::addSeparated
80 (
81  GeometricField<Type, pointPatchField, pointMesh>& pf
82 ) const
83 {
84  if (debug)
85  {
86  Pout<< "volPointInterpolation::addSeparated" << endl;
87  }
88 
89  typename GeometricField<Type, pointPatchField, pointMesh>::
90  Internal& pfi = pf.ref();
91 
92  typename GeometricField<Type, pointPatchField, pointMesh>::
93  Boundary& pfbf = pf.boundaryFieldRef();
94 
95  forAll(pfbf, patchi)
96  {
97  if (pfbf[patchi].coupled())
98  {
99  refCast<coupledPointPatchField<Type>>
100  (pfbf[patchi]).initSwapAddSeparated
101  (
103  pfi
104  );
105  }
106  }
107 
108  // Block for any outstanding requests
110 
111  forAll(pfbf, patchi)
112  {
113  if (pfbf[patchi].coupled())
114  {
115  refCast<coupledPointPatchField<Type>>
116  (pfbf[patchi]).swapAddSeparated
117  (
119  pfi
120  );
121  }
122  }
123 }
124 
125 
126 template<class Type>
127 void Foam::volPointInterpolation::interpolateInternalField
128 (
129  const GeometricField<Type, fvPatchField, volMesh>& vf,
130  GeometricField<Type, pointPatchField, pointMesh>& pf
131 ) const
132 {
133  if (debug)
134  {
135  Pout<< "volPointInterpolation::interpolateInternalField("
136  << "const GeometricField<Type, fvPatchField, volMesh>&, "
137  << "GeometricField<Type, pointPatchField, pointMesh>&) : "
138  << "interpolating field from cells to points"
139  << endl;
140  }
141 
142  const labelListList& pointCells = vf.mesh().pointCells();
143 
144  // Multiply volField by weighting factor matrix to create pointField
145  forAll(pointCells, pointi)
146  {
147  if (!isPatchPoint_[pointi])
148  {
149  const scalarList& pw = pointWeights_[pointi];
150  const labelList& ppc = pointCells[pointi];
151 
152  pf[pointi] = Zero;
153 
154  forAll(ppc, pointCelli)
155  {
156  pf[pointi] += pw[pointCelli]*vf[ppc[pointCelli]];
157  }
158  }
159  }
160 }
161 
162 
163 template<class Type>
164 Foam::tmp<Foam::Field<Type>> Foam::volPointInterpolation::flatBoundaryField
165 (
166  const GeometricField<Type, fvPatchField, volMesh>& vf
167 ) const
168 {
169  const fvMesh& mesh = vf.mesh();
170  const fvBoundaryMesh& bm = mesh.boundary();
171 
172  tmp<Field<Type>> tboundaryVals
173  (
174  new Field<Type>(mesh.nFaces()-mesh.nInternalFaces())
175  );
176  Field<Type>& boundaryVals = tboundaryVals.ref();
177 
178  forAll(vf.boundaryField(), patchi)
179  {
180  label bFacei = bm[patchi].patch().start() - mesh.nInternalFaces();
181 
182  if
183  (
184  !isA<emptyFvPatch>(bm[patchi])
185  && !vf.boundaryField()[patchi].coupled()
186  )
187  {
188  SubList<Type>
189  (
190  boundaryVals,
191  vf.boundaryField()[patchi].size(),
192  bFacei
193  ) = vf.boundaryField()[patchi];
194  }
195  else
196  {
197  const polyPatch& pp = bm[patchi].patch();
198 
199  forAll(pp, i)
200  {
201  boundaryVals[bFacei++] = Zero;
202  }
203  }
204  }
205 
206  return tboundaryVals;
207 }
208 
209 
210 template<class Type>
211 void Foam::volPointInterpolation::interpolateBoundaryField
212 (
213  const GeometricField<Type, fvPatchField, volMesh>& vf,
214  GeometricField<Type, pointPatchField, pointMesh>& pf
215 ) const
216 {
217  const primitivePatch& boundary = boundaryPtr_();
218 
219  Field<Type>& pfi = pf.primitiveFieldRef();
220 
221  // Get face data in flat list
222  tmp<Field<Type>> tboundaryVals(flatBoundaryField(vf));
223  const Field<Type>& boundaryVals = tboundaryVals();
224 
225  // Do points on 'normal' patches from the surrounding patch faces
226  forAll(boundary.meshPoints(), i)
227  {
228  const label pointi = boundary.meshPoints()[i];
229 
230  if (isPatchPoint_[pointi])
231  {
232  const labelList& pFaces = boundary.pointFaces()[i];
233  const scalarList& pWeights = boundaryPointWeights_[i];
234 
235  Type& val = pfi[pointi];
236 
237  val = Zero;
238  forAll(pFaces, j)
239  {
240  if (boundaryIsPatchFace_[pFaces[j]])
241  {
242  val += pWeights[j]*boundaryVals[pFaces[j]];
243  }
244  }
245  }
246  }
247 
248  // Sum collocated contributions
249  pointConstraints::syncUntransformedData(mesh(), pfi, plusEqOp<Type>());
250 
251  // And add separated contributions
252  addSeparated(pf);
253 
254  // Push master data to slaves. It is possible (not sure how often) for
255  // a coupled point to have its master on a different patch so
256  // to make sure just push master data to slaves.
257  pushUntransformedData(pfi);
258 
259 
260  // Detect whether the field has overridden constraint patch types. If not,
261  // we are done, so return.
262  bool havePatchTypes = false;
263  wordList patchTypes(pf.boundaryField().size(), word::null);
264  forAll(pf.boundaryField(), patchi)
265  {
266  if (pf.boundaryField()[patchi].overridesConstraint())
267  {
268  havePatchTypes = true;
269  patchTypes[patchi] = pf.boundaryField()[patchi].patch().type();
270  }
271  }
272  if (!havePatchTypes)
273  {
274  return;
275  }
276 
277  // If the patch types have been overridden than we need to re-normalise the
278  // boundary points weights. Re-calculate the weight sum.
279  pointScalarField psw
280  (
281  IOobject
282  (
283  "volPointSumWeights",
285  mesh()
286  ),
287  pointMesh::New(mesh()),
289  wordList
290  (
291  pf.boundaryField().size(),
292  calculatedPointPatchField<scalar>::typeName
293  ),
294  patchTypes
295  );
296 
297  scalarField& pswi = psw.primitiveFieldRef();
298 
299  forAll(boundary.meshPoints(), i)
300  {
301  const label pointi = boundary.meshPoints()[i];
302 
303  if (isPatchPoint_[pointi])
304  {
305  pswi[pointi] = sum(boundaryPointWeights_[i]);
306  }
307  }
308 
309  pointConstraints::syncUntransformedData(mesh(), pswi, plusEqOp<scalar>());
310 
311  addSeparated(psw);
312 
313  pushUntransformedData(pswi);
314 
315  // Apply the new weight sum to the result to re-normalise it
316  forAll(boundary.meshPoints(), i)
317  {
318  const label pointi = boundary.meshPoints()[i];
319 
320  if (isPatchPoint_[pointi])
321  {
322  pfi[pointi] /= pswi[pointi];
323  }
324  }
325 }
326 
327 
328 template<class Type>
329 void Foam::volPointInterpolation::interpolateBoundaryField
330 (
331  const GeometricField<Type, fvPatchField, volMesh>& vf,
332  GeometricField<Type, pointPatchField, pointMesh>& pf,
333  const bool overrideFixedValue
334 ) const
335 {
336  interpolateBoundaryField(vf, pf);
337 
338  // Apply constraints
339  const pointConstraints& pcs = pointConstraints::New(pf.mesh());
340 
341  pcs.constrain(pf, overrideFixedValue);
342 }
343 
344 
345 template<class Type>
347 (
350 ) const
351 {
352  if (debug)
353  {
354  Pout<< "volPointInterpolation::interpolate("
355  << "const GeometricField<Type, fvPatchField, volMesh>&, "
356  << "GeometricField<Type, pointPatchField, pointMesh>&) : "
357  << "interpolating field from cells to points"
358  << endl;
359  }
360 
361  interpolateInternalField(vf, pf);
362 
363  // Interpolate to the patches preserving fixed value BCs
364  interpolateBoundaryField(vf, pf, false);
365 }
366 
367 
368 template<class Type>
371 (
373  const wordList& patchFieldTypes
374 ) const
375 {
376  const pointMesh& pm = pointMesh::New(vf.mesh());
377 
378  // Construct tmp<pointField>
380  (
382  (
383  "volPointInterpolate(" + vf.name() + ')',
384  pm,
385  vf.dimensions(),
386  patchFieldTypes
387  )
388  );
389 
390  interpolateInternalField(vf, tpf());
391 
392  // Interpolate to the patches overriding fixed value BCs
393  interpolateBoundaryField(vf, tpf(), true);
394 
395  return tpf;
396 }
397 
398 
399 template<class Type>
402 (
404  const wordList& patchFieldTypes
405 ) const
406 {
407  // Construct tmp<pointField>
409  interpolate(tvf(), patchFieldTypes);
410  tvf.clear();
411  return tpf;
412 }
413 
414 
415 template<class Type>
418 (
420  const word& name,
421  const bool cache
422 ) const
423 {
425 
426  const pointMesh& pm = pointMesh::New(vf.mesh());
427  const objectRegistry& db = pm.thisDb();
428 
429  if (!cache || vf.mesh().changing())
430  {
431  // Delete any old occurrences to avoid double registration
432  if (db.objectRegistry::template foundObject<PointFieldType>(name))
433  {
434  PointFieldType& pf =
435  db.objectRegistry::template lookupObjectRef<PointFieldType>
436  (
437  name
438  );
439 
440  if (pf.ownedByRegistry())
441  {
442  solution::cachePrintMessage("Deleting", name, vf);
443  pf.release();
444  delete &pf;
445  }
446  }
447 
448 
450  (
452  (
453  name,
454  pm,
455  vf.dimensions()
456  )
457  );
458 
459  interpolate(vf, tpf.ref());
460 
461  return tpf;
462  }
463  else
464  {
465  if (!db.objectRegistry::template foundObject<PointFieldType>(name))
466  {
467  solution::cachePrintMessage("Calculating and caching", name, vf);
468  tmp<PointFieldType> tpf = interpolate(vf, name, false);
469  PointFieldType* pfPtr = tpf.ptr();
470  regIOobject::store(pfPtr);
471  return *pfPtr;
472  }
473  else
474  {
475  PointFieldType& pf =
476  db.objectRegistry::template lookupObjectRef<PointFieldType>
477  (
478  name
479  );
480 
481  if (pf.upToDate(vf)) // TBD: , vf.mesh().points()))
482  {
483  solution::cachePrintMessage("Reusing", name, vf);
484  return pf;
485  }
486  else
487  {
488  solution::cachePrintMessage("Deleting", name, vf);
489  pf.release();
490  delete &pf;
491 
492  solution::cachePrintMessage("Recalculating", name, vf);
493  tmp<PointFieldType> tpf = interpolate(vf, name, false);
494 
495  solution::cachePrintMessage("Storing", name, vf);
496  PointFieldType* pfPtr = tpf.ptr();
497  regIOobject::store(pfPtr);
498 
499  // Note: return reference, not pointer
500  return *pfPtr;
501  }
502  }
503  }
504 }
505 
506 
507 template<class Type>
510 (
512 ) const
513 {
514  return interpolate(vf, "volPointInterpolate(" + vf.name() + ')', false);
515 }
516 
517 
518 template<class Type>
521 (
523 ) const
524 {
525  // Construct tmp<pointField>
527  interpolate(tvf());
528  tvf.clear();
529  return tpf;
530 }
531 
532 
533 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
tmp< GeometricField< Type, pointPatchField, pointMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &) const
Interpolate volField using inverse distance weighting.
#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
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
IOobject(const word &name, const fileName &instance, const objectRegistry &registry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
Definition: IOobject.C:167
wordList patchTypes(nPatches)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Generic GeometricField class.
const dimensionSet dimless
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
const dimensionSet & dimensions() const
Return dimensions.
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static void syncUntransformedData(const polyMesh &mesh, List< Type > &pointData, const CombineOp &cop)
Helper: sync data on collocated points only.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
const objectRegistry & thisDb() const
Return database. For now is its polyMesh.
Definition: pointMesh.H:109
static const word null
An empty word.
Definition: word.H:77
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1435
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:40
const Mesh & mesh() const
Return mesh.
static void waitRequests(const label start=0)
Wait until all requests (from start onwards) have finished.
Definition: UPstream.C:147
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
List< word > wordList
A List of words.
Definition: fileName.H:54
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
fileName & instance() const
Return the instance directory, constant, system, <time> etc.
Definition: IOobject.C:355
T * ptr() const
Return tmp pointer for reuse.
Definition: tmpI.H:205
A class for managing temporary objects.
Definition: PtrList.H:53
Registry of regIOobjects.
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:312
static void cachePrintMessage(const char *message, const word &name, const FieldType &vf)
Helper for printing cache message.
void constrain(GeometricField< Type, pointPatchField, pointMesh > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.