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-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 "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  const word patchType = pf.boundaryField()[patchi].patchType();
267  if (patchType != word::null)
268  {
269  havePatchTypes = true;
270  patchTypes[patchi] = patchType;
271  }
272  }
273  if (!havePatchTypes)
274  {
275  return;
276  }
277 
278  // If the patch types have been overridden than we need to re-normalise the
279  // boundary points weights. Re-calculate the weight sum.
280  pointScalarField psw
281  (
282  IOobject
283  (
284  "volPointSumWeights",
286  mesh()
287  ),
288  pointMesh::New(mesh()),
290  wordList
291  (
292  pf.boundaryField().size(),
294  ),
295  patchTypes
296  );
297 
298  scalarField& pswi = psw.primitiveFieldRef();
299 
300  forAll(boundary.meshPoints(), i)
301  {
302  const label pointi = boundary.meshPoints()[i];
303 
304  if (isPatchPoint_[pointi])
305  {
306  pswi[pointi] = sum(boundaryPointWeights_[i]);
307  }
308  }
309 
311 
312  addSeparated(psw);
313 
314  pushUntransformedData(pswi);
315 
316  // Apply the new weight sum to the result to re-normalise it
317  forAll(boundary.meshPoints(), i)
318  {
319  const label pointi = boundary.meshPoints()[i];
320 
321  if (isPatchPoint_[pointi])
322  {
323  pfi[pointi] /= pswi[pointi];
324  }
325  }
326 }
327 
328 
329 template<class Type>
331 (
334  const bool overrideFixedValue
335 ) const
336 {
337  interpolateBoundaryField(vf, pf);
338 
339  // Apply constraints
340  const pointConstraints& pcs = pointConstraints::New(pf.mesh());
341 
342  pcs.constrain(pf, overrideFixedValue);
343 }
344 
345 
346 template<class Type>
348 (
351 ) const
352 {
353  if (debug)
354  {
355  Pout<< "volPointInterpolation::interpolate("
356  << "const GeometricField<Type, fvPatchField, volMesh>&, "
357  << "GeometricField<Type, pointPatchField, pointMesh>&) : "
358  << "interpolating field from cells to points"
359  << endl;
360  }
361 
362  interpolateInternalField(vf, pf);
363 
364  // Interpolate to the patches preserving fixed value BCs
365  interpolateBoundaryField(vf, pf, false);
366 }
367 
368 
369 template<class Type>
372 (
374  const wordList& patchFieldTypes
375 ) const
376 {
377  const pointMesh& pm = pointMesh::New(vf.mesh());
378 
379  // Construct tmp<pointField>
381  (
383  (
384  IOobject
385  (
386  "volPointInterpolate(" + vf.name() + ')',
387  vf.instance(),
388  pm.thisDb()
389  ),
390  pm,
391  vf.dimensions(),
392  patchFieldTypes
393  )
394  );
395 
396  interpolateInternalField(vf, tpf());
397 
398  // Interpolate to the patches overriding fixed value BCs
399  interpolateBoundaryField(vf, tpf(), true);
400 
401  return tpf;
402 }
403 
404 
405 template<class Type>
408 (
410  const wordList& patchFieldTypes
411 ) const
412 {
413  // Construct tmp<pointField>
415  interpolate(tvf(), patchFieldTypes);
416  tvf.clear();
417  return tpf;
418 }
419 
420 
421 template<class Type>
424 (
426  const word& name,
427  const bool cache
428 ) const
429 {
431 
432  const pointMesh& pm = pointMesh::New(vf.mesh());
433  const objectRegistry& db = pm.thisDb();
434 
435  if (!cache || vf.mesh().changing())
436  {
437  // Delete any old occurrences to avoid double registration
438  if (db.objectRegistry::template foundObject<PointFieldType>(name))
439  {
440  PointFieldType& pf =
441  db.objectRegistry::template lookupObjectRef<PointFieldType>
442  (
443  name
444  );
445 
446  if (pf.ownedByRegistry())
447  {
448  solution::cachePrintMessage("Deleting", name, vf);
449  pf.release();
450  delete &pf;
451  }
452  }
453 
454 
456  (
458  (
459  IOobject
460  (
461  name,
462  vf.instance(),
463  pm.thisDb()
464  ),
465  pm,
466  vf.dimensions()
467  )
468  );
469 
470  interpolate(vf, tpf.ref());
471 
472  return tpf;
473  }
474  else
475  {
476  if (!db.objectRegistry::template foundObject<PointFieldType>(name))
477  {
478  solution::cachePrintMessage("Calculating and caching", name, vf);
479  tmp<PointFieldType> tpf = interpolate(vf, name, false);
480  PointFieldType* pfPtr = tpf.ptr();
481  regIOobject::store(pfPtr);
482  return *pfPtr;
483  }
484  else
485  {
486  PointFieldType& pf =
487  db.objectRegistry::template lookupObjectRef<PointFieldType>
488  (
489  name
490  );
491 
492  if (pf.upToDate(vf)) // TBD: , vf.mesh().points()))
493  {
494  solution::cachePrintMessage("Reusing", name, vf);
495  return pf;
496  }
497  else
498  {
499  solution::cachePrintMessage("Deleting", name, vf);
500  pf.release();
501  delete &pf;
502 
503  solution::cachePrintMessage("Recalculating", name, vf);
504  tmp<PointFieldType> tpf = interpolate(vf, name, false);
505 
506  solution::cachePrintMessage("Storing", name, vf);
507  PointFieldType* pfPtr = tpf.ptr();
508  regIOobject::store(pfPtr);
509 
510  // Note: return reference, not pointer
511  return *pfPtr;
512  }
513  }
514  }
515 }
516 
517 
518 template<class Type>
521 (
523 ) const
524 {
525  return interpolate(vf, "volPointInterpolate(" + vf.name() + ')', false);
526 }
527 
528 
529 template<class Type>
532 (
534 ) const
535 {
536  // Construct tmp<pointField>
538  interpolate(tvf());
539  tvf.clear();
540  return tpf;
541 }
542 
543 
544 // ************************************************************************* //
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:295
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:174
wordList patchTypes(nPatches)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Generic GeometricField class.
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)
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
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)
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:378
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
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:198
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:354
static void cachePrintMessage(const char *message, const word &name, const FieldType &vf)
Helper for printing cache message.
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, &oldCyclicPolyPatch::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
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