volPointInterpolation.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 "fvMesh.H"
28 #include "volFields.H"
29 #include "pointFields.H"
30 #include "demandDrivenData.H"
31 #include "pointConstraints.H"
32 #include "surfaceFields.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(volPointInterpolation, 0);
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 void Foam::volPointInterpolation::calcBoundaryAddressing()
45 {
46  if (debug)
47  {
48  Pout<< "volPointInterpolation::calcBoundaryAddressing() : "
49  << "constructing boundary addressing"
50  << endl;
51  }
52 
53  boundaryPtr_.reset
54  (
55  new primitivePatch
56  (
57  SubList<face>
58  (
59  mesh().faces(),
60  mesh().nFaces()-mesh().nInternalFaces(),
61  mesh().nInternalFaces()
62  ),
63  mesh().points()
64  )
65  );
66  const primitivePatch& boundary = boundaryPtr_();
67 
68  boundaryIsPatchFace_.setSize(boundary.size());
69  boundaryIsPatchFace_ = false;
70 
71  isPatchPoint_.setSize(mesh().nPoints());
72  isPatchPoint_ = false;
73 
74  const polyBoundaryMesh& pbm = mesh().boundaryMesh();
75 
76  // Get precalculated volField only so we can use coupled() tests for
77  // cyclicAMI
78  const surfaceScalarField& magSf = mesh().magSf();
79 
80  forAll(pbm, patchi)
81  {
82  const polyPatch& pp = pbm[patchi];
83 
84  if
85  (
86  !isA<emptyPolyPatch>(pp)
87  && !magSf.boundaryField()[patchi].coupled()
88  )
89  {
90  label bFacei = pp.start()-mesh().nInternalFaces();
91 
92  forAll(pp, i)
93  {
94  boundaryIsPatchFace_[bFacei] = true;
95 
96  const face& f = boundary[bFacei++];
97 
98  forAll(f, fp)
99  {
100  isPatchPoint_[f[fp]] = true;
101  }
102  }
103  }
104  }
105 
106  // Make sure point status is synchronised so even processor that holds
107  // no face of a certain patch still can have boundary points marked.
108  if (debug)
109  {
110  boolList oldData(isPatchPoint_);
111 
113  (
114  mesh(),
115  isPatchPoint_,
116  orEqOp<bool>()
117  );
118 
119  forAll(isPatchPoint_, pointi)
120  {
121  if (isPatchPoint_[pointi] != oldData[pointi])
122  {
123  Pout<< "volPointInterpolation::calcBoundaryAddressing():"
124  << " added dangling mesh point:" << pointi
125  << " at:" << mesh().points()[pointi]
126  << endl;
127  }
128  }
129 
130  label nPatchFace = 0;
131  forAll(boundaryIsPatchFace_, i)
132  {
133  if (boundaryIsPatchFace_[i])
134  {
135  nPatchFace++;
136  }
137  }
138  label nPatchPoint = 0;
139  forAll(isPatchPoint_, i)
140  {
141  if (isPatchPoint_[i])
142  {
143  nPatchPoint++;
144  }
145  }
146  Pout<< "boundary:" << nl
147  << " faces :" << boundary.size() << nl
148  << " of which on proper patch:" << nPatchFace << nl
149  << " points:" << boundary.nPoints() << nl
150  << " of which on proper patch:" << nPatchPoint << endl;
151  }
152 }
153 
154 
155 void Foam::volPointInterpolation::makeInternalWeights(scalarField& sumWeights)
156 {
157  if (debug)
158  {
159  Pout<< "volPointInterpolation::makeInternalWeights() : "
160  << "constructing weighting factors for internal and non-coupled"
161  << " points." << endl;
162  }
163 
164  const pointField& points = mesh().points();
165  const labelListList& pointCells = mesh().pointCells();
166  const vectorField& cellCentres = mesh().cellCentres();
167 
168  // Allocate storage for weighting factors
169  pointWeights_.clear();
170  pointWeights_.setSize(points.size());
171 
172  // Calculate inverse distances between cell centres and points
173  // and store in weighting factor array
174  forAll(points, pointi)
175  {
176  if (!isPatchPoint_[pointi])
177  {
178  const labelList& pcp = pointCells[pointi];
179 
180  scalarList& pw = pointWeights_[pointi];
181  pw.setSize(pcp.size());
182 
183  forAll(pcp, pointCelli)
184  {
185  pw[pointCelli] =
186  1.0/mag(points[pointi] - cellCentres[pcp[pointCelli]]);
187 
188  sumWeights[pointi] += pw[pointCelli];
189  }
190  }
191  }
192 }
193 
194 
195 void Foam::volPointInterpolation::makeBoundaryWeights(scalarField& sumWeights)
196 {
197  if (debug)
198  {
199  Pout<< "volPointInterpolation::makeBoundaryWeights() : "
200  << "constructing weighting factors for boundary points." << endl;
201  }
202 
203  const pointField& points = mesh().points();
204  const pointField& faceCentres = mesh().faceCentres();
205 
206  const primitivePatch& boundary = boundaryPtr_();
207 
208  boundaryPointWeights_.clear();
209  boundaryPointWeights_.setSize(boundary.meshPoints().size());
210 
211  forAll(boundary.meshPoints(), i)
212  {
213  label pointi = boundary.meshPoints()[i];
214 
215  if (isPatchPoint_[pointi])
216  {
217  const labelList& pFaces = boundary.pointFaces()[i];
218 
219  scalarList& pw = boundaryPointWeights_[i];
220  pw.setSize(pFaces.size());
221 
222  sumWeights[pointi] = 0.0;
223 
224  forAll(pFaces, i)
225  {
226  if (boundaryIsPatchFace_[pFaces[i]])
227  {
228  label facei = mesh().nInternalFaces() + pFaces[i];
229 
230  pw[i] = 1.0/mag(points[pointi] - faceCentres[facei]);
231  sumWeights[pointi] += pw[i];
232  }
233  else
234  {
235  pw[i] = 0.0;
236  }
237  }
238  }
239  }
240 }
241 
242 
243 void Foam::volPointInterpolation::makeWeights()
244 {
245  if (debug)
246  {
247  Pout<< "volPointInterpolation::makeWeights() : "
248  << "constructing weighting factors"
249  << endl;
250  }
251 
252  // Update addressing over all boundary faces
253  calcBoundaryAddressing();
254 
255 
256  // Running sum of weights
257  pointScalarField sumWeights
258  (
259  IOobject
260  (
261  "volPointSumWeights",
263  mesh()
264  ),
265  pointMesh::New(mesh()),
267  );
268 
269 
270  // Create internal weights; add to sumWeights
271  makeInternalWeights(sumWeights);
272 
273 
274  // Create boundary weights; override sumWeights
275  makeBoundaryWeights(sumWeights);
276 
277 
278  // forAll(boundary.meshPoints(), i)
279  //{
280  // label pointi = boundary.meshPoints()[i];
281  //
282  // if (isPatchPoint_[pointi])
283  // {
284  // Pout<< "Calculated Weight at boundary point:" << i
285  // << " at:" << mesh().points()[pointi]
286  // << " sumWeight:" << sumWeights[pointi]
287  // << " from:" << boundaryPointWeights_[i]
288  // << endl;
289  // }
290  //}
291 
292 
293  // Sum collocated contributions
295  (
296  mesh(),
297  sumWeights,
298  plusEqOp<scalar>()
299  );
300 
301  // And add separated contributions
302  addSeparated(sumWeights);
303 
304  // Push master data to slaves. It is possible (not sure how often) for
305  // a coupled point to have its master on a different patch so
306  // to make sure just push master data to slaves. Reuse the syncPointData
307  // structure.
308  pushUntransformedData(sumWeights);
309 
310 
311  // Normalise internal weights
312  forAll(pointWeights_, pointi)
313  {
314  scalarList& pw = pointWeights_[pointi];
315  // Note:pw only sized for !isPatchPoint
316  forAll(pw, i)
317  {
318  pw[i] /= sumWeights[pointi];
319  }
320  }
321 
322  // Normalise boundary weights
323  const primitivePatch& boundary = boundaryPtr_();
324 
325  forAll(boundary.meshPoints(), i)
326  {
327  label pointi = boundary.meshPoints()[i];
328 
329  scalarList& pw = boundaryPointWeights_[i];
330  // Note:pw only sized for isPatchPoint
331  forAll(pw, i)
332  {
333  pw[i] /= sumWeights[pointi];
334  }
335  }
336 
337 
338  if (debug)
339  {
340  Pout<< "volPointInterpolation::makeWeights() : "
341  << "finished constructing weighting factors"
342  << endl;
343  }
344 }
345 
346 
347 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
348 
350 :
352 {
353  makeWeights();
354 }
355 
356 
357 // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
358 
360 {}
361 
362 
363 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
364 
366 {
367  makeWeights();
368 }
369 
370 
372 {
373  makeWeights();
374 
375  return true;
376 }
377 
378 
380 (
381  const volVectorField& vf,
382  pointVectorField& pf
383 ) const
384 {
385  interpolateInternalField(vf, pf);
386 
387  // Interpolate to the patches but no constraints
388  interpolateBoundaryField(vf, pf);
389 
390  // Apply displacement constraints
391  const pointConstraints& pcs = pointConstraints::New(pf.mesh());
392 
393  pcs.constrainDisplacement(pf, false);
394 }
395 
396 
397 // ************************************************************************* //
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
Foam::surfaceFields.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void updateMesh(const mapPolyMesh &)
Update mesh topology using the morph engine.
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
void interpolateBoundaryField(const GeometricField< Type, fvPatchField, volMesh > &vf, GeometricField< Type, pointPatchField, pointMesh > &pf) const
Interpolate boundary field without applying constraints/boundary.
label nInternalFaces() const
void constrainDisplacement(pointVectorField &displacement, const bool overrideValue=false) const
Apply boundary conditions (single-patch constraints),.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Application of (multi-)patch point constraints.
Templated abstract base-class for optional mesh objects used to automate their allocation to the mesh...
Definition: MeshObject.H:86
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dynamicFvMesh & mesh
const pointField & points
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
label nPoints
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
List< label > labelList
A List of labels.
Definition: labelList.H:56
void interpolateInternalField(const GeometricField< Type, fvPatchField, volMesh > &, GeometricField< Type, pointPatchField, pointMesh > &) const
Interpolate internal field from volField to pointField.
const vectorField & cellCentres() const
Interpolate from cell centres to points (vertices) using inverse distance weighting.
const labelListList & pointCells() const
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
static const char nl
Definition: Ostream.H:265
const Mesh & mesh() const
Return mesh.
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
void interpolateDisplacement(const volVectorField &, pointVectorField &) const
Interpolate from volField to pointField.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
volPointInterpolation(const fvMesh &)
Constructor given fvMesh and pointMesh.
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
Template functions to aid in the implementation of demand driven data.
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.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
bool movePoints()
Correct weighting factors for moving mesh.
Namespace for OpenFOAM.