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-2021 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 mapDistribute& 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>
128 (
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 (
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 
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  {
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>
212 (
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
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;
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(),
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 
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>
330 (
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
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:303
void interpolateBoundaryField(const GeometricField< Type, fvPatchField, volMesh > &vf, GeometricField< Type, pointPatchField, pointMesh > &pf) const
Interpolate boundary field without applying constraints/boundary.
label nInternalFaces() const
const Boundary & boundaryField() const
Return const-reference to the boundary field.
label nFaces() const
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
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.
A calculated boundary condition for pointField.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
Application of (multi-)patch point constraints.
const dimensionSet & dimensions() const
Return dimensions.
A list of faces which address into the list of points.
Pre-declare SubField and related Field type.
Definition: Field.H:56
A class for handling words, derived from string.
Definition: word.H:59
static void syncUntransformedData(const polyMesh &mesh, List< Type > &pointData, const CombineOp &cop)
Helper: sync data on collocated points only.
const objectRegistry & thisDb() const
Return database. For now is its polyMesh.
Definition: pointMesh.H:115
static const word null
An empty word.
Definition: word.H:77
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1394
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
void interpolateInternalField(const GeometricField< Type, fvPatchField, volMesh > &, GeometricField< Type, pointPatchField, pointMesh > &) const
Interpolate internal field from volField to pointField.
faceListList boundary(nPatches)
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:34
const labelListList & pointFaces() const
Return point-face addressing.
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
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
const fileName & instance() const
Definition: IOobject.H:390
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Foam::fvBoundaryMesh.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
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.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
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.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
void constrain(GeometricField< Type, pointPatchField, pointMesh > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:540