volPointInterpolate.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) 2011-2013 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 namespace Foam
36 {
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 template<class Type>
41 void volPointInterpolation::pushUntransformedData
42 (
43  List<Type>& pointData
44 ) const
45 {
46  // Transfer onto coupled patch
47  const globalMeshData& gmd = mesh().globalData();
48  const indirectPrimitivePatch& cpp = gmd.coupledPatch();
49  const labelList& meshPoints = cpp.meshPoints();
50 
51  const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
52  const labelListList& slaves = gmd.globalCoPointSlaves();
53 
54  List<Type> elems(slavesMap.constructSize());
55  forAll(meshPoints, i)
56  {
57  elems[i] = pointData[meshPoints[i]];
58  }
59 
60  // Combine master data with slave data
61  forAll(slaves, i)
62  {
63  const labelList& slavePoints = slaves[i];
64 
65  // Copy master data to slave slots
66  forAll(slavePoints, j)
67  {
68  elems[slavePoints[j]] = elems[i];
69  }
70  }
71 
72  // Push slave-slot data back to slaves
73  slavesMap.reverseDistribute(elems.size(), elems, false);
74 
75  // Extract back onto mesh
76  forAll(meshPoints, i)
77  {
78  pointData[meshPoints[i]] = elems[i];
79  }
80 }
81 
82 
83 template<class Type>
84 void volPointInterpolation::addSeparated
85 (
86  GeometricField<Type, pointPatchField, pointMesh>& pf
87 ) const
88 {
89  if (debug)
90  {
91  Pout<< "volPointInterpolation::addSeparated" << endl;
92  }
93 
94  forAll(pf.boundaryField(), patchI)
95  {
96  if (pf.boundaryField()[patchI].coupled())
97  {
98  refCast<coupledPointPatchField<Type> >
99  (pf.boundaryField()[patchI]).initSwapAddSeparated
100  (
102  pf.internalField()
103  );
104  }
105  }
106 
107  // Block for any outstanding requests
109 
110  forAll(pf.boundaryField(), patchI)
111  {
112  if (pf.boundaryField()[patchI].coupled())
113  {
114  refCast<coupledPointPatchField<Type> >
115  (pf.boundaryField()[patchI]).swapAddSeparated
116  (
118  pf.internalField()
119  );
120  }
121  }
122 }
123 
124 
125 template<class Type>
127 (
130 ) const
131 {
132  if (debug)
133  {
134  Pout<< "volPointInterpolation::interpolateInternalField("
135  << "const GeometricField<Type, fvPatchField, volMesh>&, "
136  << "GeometricField<Type, pointPatchField, pointMesh>&) : "
137  << "interpolating field from cells to points"
138  << endl;
139  }
140 
141  const labelListList& pointCells = vf.mesh().pointCells();
142 
143  // Multiply volField by weighting factor matrix to create pointField
144  forAll(pointCells, pointi)
145  {
146  if (!isPatchPoint_[pointi])
147  {
148  const scalarList& pw = pointWeights_[pointi];
149  const labelList& ppc = pointCells[pointi];
150 
151  pf[pointi] = pTraits<Type>::zero;
152 
153  forAll(ppc, pointCelli)
154  {
155  pf[pointi] += pw[pointCelli]*vf[ppc[pointCelli]];
156  }
157  }
158  }
159 }
160 
161 
162 template<class Type>
163 tmp<Field<Type> > volPointInterpolation::flatBoundaryField
164 (
166 ) const
167 {
168  const fvMesh& mesh = vf.mesh();
169  const fvBoundaryMesh& bm = mesh.boundary();
170 
171  tmp<Field<Type> > tboundaryVals
172  (
173  new Field<Type>(mesh.nFaces()-mesh.nInternalFaces())
174  );
175  Field<Type>& boundaryVals = tboundaryVals();
176 
177  forAll(vf.boundaryField(), patchI)
178  {
179  label bFaceI = bm[patchI].patch().start() - mesh.nInternalFaces();
180 
181  if
182  (
183  !isA<emptyFvPatch>(bm[patchI])
184  && !vf.boundaryField()[patchI].coupled()
185  )
186  {
188  (
189  boundaryVals,
190  vf.boundaryField()[patchI].size(),
191  bFaceI
192  ).assign(vf.boundaryField()[patchI]);
193  }
194  else
195  {
196  const polyPatch& pp = bm[patchI].patch();
197 
198  forAll(pp, i)
199  {
200  boundaryVals[bFaceI++] = pTraits<Type>::zero;
201  }
202  }
203  }
204 
205  return tboundaryVals;
206 }
207 
208 
209 template<class Type>
211 (
214 ) const
215 {
216  const primitivePatch& boundary = boundaryPtr_();
217 
218  Field<Type>& pfi = pf.internalField();
219 
220  // Get face data in flat list
221  tmp<Field<Type> > tboundaryVals(flatBoundaryField(vf));
222  const Field<Type>& boundaryVals = tboundaryVals();
223 
224 
225  // Do points on 'normal' patches from the surrounding patch faces
226  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
227 
228  forAll(boundary.meshPoints(), i)
229  {
230  label pointI = boundary.meshPoints()[i];
231 
232  if (isPatchPoint_[pointI])
233  {
234  const labelList& pFaces = boundary.pointFaces()[i];
235  const scalarList& pWeights = boundaryPointWeights_[i];
236 
237  Type& val = pfi[pointI];
238 
239  val = pTraits<Type>::zero;
240  forAll(pFaces, j)
241  {
242  if (boundaryIsPatchFace_[pFaces[j]])
243  {
244  val += pWeights[j]*boundaryVals[pFaces[j]];
245  }
246  }
247  }
248  }
249 
250  // Sum collocated contributions
252 
253  // And add separated contributions
254  addSeparated(pf);
255 
256  // Push master data to slaves. It is possible (not sure how often) for
257  // a coupled point to have its master on a different patch so
258  // to make sure just push master data to slaves.
259  pushUntransformedData(pfi);
260 }
261 
262 
263 template<class Type>
265 (
268  const bool overrideFixedValue
269 ) const
270 {
271  interpolateBoundaryField(vf, pf);
272 
273  // Apply constraints
274  const pointConstraints& pcs = pointConstraints::New(pf.mesh());
275 
276  pcs.constrain(pf, overrideFixedValue);
277 }
278 
279 
280 template<class Type>
282 (
285 ) const
286 {
287  if (debug)
288  {
289  Pout<< "volPointInterpolation::interpolate("
290  << "const GeometricField<Type, fvPatchField, volMesh>&, "
291  << "GeometricField<Type, pointPatchField, pointMesh>&) : "
292  << "interpolating field from cells to points"
293  << endl;
294  }
295 
296  interpolateInternalField(vf, pf);
297 
298  // Interpolate to the patches preserving fixed value BCs
299  interpolateBoundaryField(vf, pf, false);
300 }
301 
302 
303 template<class Type>
306 (
308  const wordList& patchFieldTypes
309 ) const
310 {
311  const pointMesh& pm = pointMesh::New(vf.mesh());
312 
313  // Construct tmp<pointField>
315  (
317  (
318  IOobject
319  (
320  "volPointInterpolate(" + vf.name() + ')',
321  vf.instance(),
322  pm.thisDb()
323  ),
324  pm,
325  vf.dimensions(),
326  patchFieldTypes
327  )
328  );
329 
330  interpolateInternalField(vf, tpf());
331 
332  // Interpolate to the patches overriding fixed value BCs
333  interpolateBoundaryField(vf, tpf(), true);
334 
335  return tpf;
336 }
337 
338 
339 template<class Type>
342 (
344  const wordList& patchFieldTypes
345 ) const
346 {
347  // Construct tmp<pointField>
349  interpolate(tvf(), patchFieldTypes);
350  tvf.clear();
351  return tpf;
352 }
353 
354 
355 template<class Type>
358 (
360  const word& name,
361  const bool cache
362 ) const
363 {
365 
366  const pointMesh& pm = pointMesh::New(vf.mesh());
367  const objectRegistry& db = pm.thisDb();
368 
369  if (!cache || vf.mesh().changing())
370  {
371  // Delete any old occurences to avoid double registration
372  if (db.objectRegistry::template foundObject<PointFieldType>(name))
373  {
374  PointFieldType& pf = const_cast<PointFieldType&>
375  (
376  db.objectRegistry::template lookupObject<PointFieldType>(name)
377  );
378 
379  if (pf.ownedByRegistry())
380  {
381  solution::cachePrintMessage("Deleting", name, vf);
382  pf.release();
383  delete &pf;
384  }
385  }
386 
387 
389  (
391  (
392  IOobject
393  (
394  name,
395  vf.instance(),
396  pm.thisDb()
397  ),
398  pm,
399  vf.dimensions()
400  )
401  );
402 
403  interpolate(vf, tpf());
404 
405  return tpf;
406  }
407  else
408  {
409  if (!db.objectRegistry::template foundObject<PointFieldType>(name))
410  {
411  solution::cachePrintMessage("Calculating and caching", name, vf);
412  tmp<PointFieldType> tpf = interpolate(vf, name, false);
413  PointFieldType* pfPtr = tpf.ptr();
414  regIOobject::store(pfPtr);
415  return *pfPtr;
416  }
417  else
418  {
419  PointFieldType& pf = const_cast<PointFieldType&>
420  (
421  db.objectRegistry::template lookupObject<PointFieldType>(name)
422  );
423 
424  if (pf.upToDate(vf)) //TBD: , vf.mesh().points()))
425  {
426  solution::cachePrintMessage("Reusing", name, vf);
427  return pf;
428  }
429  else
430  {
431  solution::cachePrintMessage("Deleting", name, vf);
432  pf.release();
433  delete &pf;
434 
435  solution::cachePrintMessage("Recalculating", name, vf);
436  tmp<PointFieldType> tpf = interpolate(vf, name, false);
437 
438  solution::cachePrintMessage("Storing", name, vf);
439  PointFieldType* pfPtr = tpf.ptr();
440  regIOobject::store(pfPtr);
441 
442  // Note: return reference, not pointer
443  return *pfPtr;
444  }
445  }
446  }
447 }
448 
449 
450 template<class Type>
453 (
455 ) const
456 {
457  return interpolate(vf, "volPointInterpolate(" + vf.name() + ')', false);
458 }
459 
460 
461 template<class Type>
464 (
466 ) const
467 {
468  // Construct tmp<pointField>
470  interpolate(tvf());
471  tvf.clear();
472  return tpf;
473 }
474 
475 
476 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
477 
478 } // End namespace Foam
479 
480 // ************************************************************************* //
const objectRegistry & thisDb() const
Return database. For now is its polyMesh.
Definition: pointMesh.H:118
void constrain(GeometricField< Type, pointPatchField, pointMesh > &pf, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints) and.
T * ptr() const
Return tmp pointer for reuse.
Definition: tmpI.H:148
faceListList boundary(nPatches)
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
label nFaces() const
tmp< GeometricField< Type, pointPatchField, pointMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &) const
Interpolate volField using inverse distance weighting.
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
static const pointConstraints & New(const pointMesh &mesh)
const labelListList & pointFaces() const
Return point-face addressing.
Foam::fvBoundaryMesh.
A class for handling words, derived from string.
Definition: word.H:59
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 fileName & instance() const
Definition: IOobject.H:337
InternalField & internalField()
Return internal field.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
const Mesh & mesh() const
Return mesh.
Namespace for OpenFOAM.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
static void syncUntransformedData(const polyMesh &mesh, List< Type > &pointData, const CombineOp &cop)
Helper: sync data on collocated points only.
static void waitRequests(const label start=0)
Wait until all requests (from start onwards) have finished.
Definition: UPstream.C:106
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
void store()
Transfer ownership of this object to its registry.
Definition: regIOobjectI.H:34
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
#define forAll(list, i)
Definition: UList.H:421
const dimensionSet & dimensions() const
Return dimensions.
Pre-declare SubField and related Field type.
Definition: Field.H:57
const word & name() const
Return name.
Definition: IOobject.H:260
static void cachePrintMessage(const char *message, const word &name, const FieldType &vf)
Helper for printing cache message.
A list of faces which address into the list of points.
label nInternalFaces() const
Generic GeometricField class.
Traits class for primitives.
Definition: pTraits.H:50
Registry of regIOobjects.
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:245
Application of (multi-)patch point contraints.
void interpolateInternalField(const GeometricField< Type, fvPatchField, volMesh > &, GeometricField< Type, pointPatchField, pointMesh > &) const
Interpolate internal field from volField to pointField.
A class for managing temporary objects.
Definition: PtrList.H:118
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1200
void interpolateBoundaryField(const GeometricField< Type, fvPatchField, volMesh > &vf, GeometricField< Type, pointPatchField, pointMesh > &pf) const
Interpolate boundary field without applying constraints/boundary.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:552