interpolationCellPointWallModified.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 
27 #include "syncTools.H"
28 #include "volPointInterpolation.H"
29 #include "wallPolyPatch.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 template<class Type>
35 template<class TYPE>
38 (
40 ) const
41 {
43  << typeName << " interpolation is only defined for vector fields"
44  << exit(FatalError);
45 
47 }
48 
49 
50 template<class Type>
53 (
54  const volVectorField& psi
55 ) const
56 {
57  using namespace constant::mathematical;
58 
59  const fvMesh& mesh = psi.mesh();
60 
61 
62  // Create the point field
64  (
66  (
67  "volPointInterpolateWallModified(" + psi.name() + ')',
68  pointMesh::New(mesh),
69  dimensioned<Type>("zero", psi.dimensions(), Zero)
70  )
71  );
72  pointVectorField& psip = tPsip.ref();
73 
74 
75  // Interpolate to the points with wall patches extrapolated
76  {
77  wordList patchTypes(psi.boundaryField().size());
79  {
80  if (isA<wallPolyPatch>(mesh.boundaryMesh()[patchi]))
81  {
82  patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName;
83  }
84  else
85  {
86  patchTypes[patchi] = calculatedFvPatchVectorField::typeName;
87  }
88  }
89  volVectorField psiExtrapolated
90  (
91  IOobject
92  (
93  psi.name() + "Extrapolated",
94  mesh.time().timeName(),
95  mesh
96  ),
97  psi,
99  );
100  psiExtrapolated.correctBoundaryConditions();
101  volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip);
102  }
103 
104 
105  // Generate point normals across the entire boundary
106  pointField pointNormals(mesh.nPoints(), vector::zero);
107  {
108  scalarField pointCount(mesh.nPoints(), 0);
109 
110  forAll(mesh.boundaryMesh(), patchi)
111  {
112  const polyPatch& patch = mesh.boundaryMesh()[patchi];
113 
114  forAll(patch, patchFacei)
115  {
116  const face& f = patch[patchFacei];
117  const vector& n = patch.faceNormals()[patchFacei];
118 
119  forAll(f, i)
120  {
121  pointNormals[f[i]] += n;
122  pointCount[f[i]] += 1;
123  }
124  }
125  }
126 
127  syncTools::syncPointList
128  (
129  mesh,
130  pointNormals,
132  vector::zero
133  );
134 
135  syncTools::syncPointList
136  (
137  mesh,
138  pointCount,
140  scalar(0)
141  );
142 
143  pointNormals /= max(pointCount, small);
144  }
145 
146 
147  // Calculate the rotation necessary from the vector to the (negative) point
148  // normal needed to make the interpolated field point into the mesh
149  scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat);
150  scalar maxVHatDotN = - vGreat;
151  forAll(mesh.boundaryMesh(), patchi)
152  {
153  const polyPatch& patch = mesh.boundaryMesh()[patchi];
154  if (isA<wallPolyPatch>(patch))
155  {
156  forAll(patch, patchFacei)
157  {
158  const label facei = patch.start() + patchFacei;
159  const face& f = patch[patchFacei];
160 
161  for (label i = 1; i < f.size() - 1; ++ i)
162  {
163  const triFace tri
164  (
165  tetIndices(mesh.faceOwner()[facei], facei, i)
166  .faceTriIs(mesh)
167  );
168 
169  const vector n = tri.normal(mesh.points());
170 
171  forAll(tri, triPointI)
172  {
173  const label pointi = tri[triPointI];
174 
175  const vector& v = psip[pointi];
176 
177  const scalar vHatDotN = normalised(v) & n;
178  maxVHatDotN = max(maxVHatDotN, vHatDotN);
179 
180  const vector a =
181  normalised(v ^ pointNormals[pointi]);
182 
183  const scalar C = v & n, S = (v ^ a) & n;
184 
185  const scalar theta = atan2(C, - S);
186 
187  theta0[pointi] = max(theta0[pointi], theta);
188  theta1[pointi] = min(theta1[pointi], theta + pi);
189  }
190  }
191  }
192  }
193  }
194 
195  if (debug)
196  {
197  reduce(maxVHatDotN, maxOp<scalar>());
198  Info<< typeName << ": Maximum in-to-wall dot product before = "
199  << maxVHatDotN << endl;
200  }
201 
202  syncTools::syncPointList
203  (
204  mesh,
205  theta0,
206  maxEqOp<scalar>(),
207  scalar(0)
208  );
209 
210  syncTools::syncPointList
211  (
212  mesh,
213  theta1,
214  minEqOp<scalar>(),
215  scalar(0)
216  );
217 
218 
219  if (debug > 1)
220  {
221  pointVectorField(psip.name() + "Before", psip).write();
222  }
223 
224 
225  // Apply the rotations so that the interpolated field points into the mesh
226  forAll(mesh.boundaryMesh(), patchi)
227  {
228  const polyPatch& patch = mesh.boundaryMesh()[patchi];
229  if (isA<wallPolyPatch>(patch))
230  {
231  forAll(patch.meshPoints(), patchPointi)
232  {
233  const label pointi = patch.meshPoints()[patchPointi];
234 
235  vector& v = psip[pointi];
236 
237  if (theta0[pointi] <= 0 && theta1[pointi] >= 0)
238  {
239  continue;
240  }
241 
242  if (theta0[pointi] >= theta1[pointi])
243  {
244  v = Zero;
245  continue;
246  }
247 
248  const scalar theta =
249  theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi];
250 
251  const scalar c = cos(theta), s = sin(theta);
252 
253  const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ...
254 
255  const vector a = normalised(v ^ pointNormals[pointi]);
256 
257  v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v;
258 
259  theta0[pointi] = theta1[pointi] = 0;
260  }
261  }
262  }
263 
264 
265  // Report the field-normal dot products
266  if (debug)
267  {
268  maxVHatDotN = - vGreat;
269 
270  forAll(mesh.boundaryMesh(), patchi)
271  {
272  const polyPatch& patch = mesh.boundaryMesh()[patchi];
273  if (isA<wallPolyPatch>(patch))
274  {
275  forAll(patch, patchFacei)
276  {
277  const label facei = patch.start() + patchFacei;
278  const face& f = patch[patchFacei];
279 
280  for (label i = 1; i < f.size() - 1; ++ i)
281  {
282  const triFace tri
283  (
284  tetIndices(mesh.faceOwner()[facei], facei, i)
285  .faceTriIs(mesh)
286  );
287 
288  const vector n = tri.normal(mesh.points());
289 
290  forAll(tri, triPointI)
291  {
292  const label pointi = tri[triPointI];
293 
294  const vector& v = psip[pointi];
295 
296  const scalar vHatDotN = normalised(v) & n;
297  maxVHatDotN = max(maxVHatDotN, vHatDotN);
298  }
299  }
300  }
301  }
302  }
303 
304  reduce(maxVHatDotN, maxOp<scalar>());
305  Info<< typeName << ": Maximum in-to-wall dot product after = "
306  << maxVHatDotN << endl;
307  }
308 
309 
310  // Write the interpolated field
311  if (debug > 1)
312  {
313  psip.write();
314  }
315 
316 
317  return tPsip;
318 }
319 
320 
321 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
322 
323 template<class Type>
326 (
328 )
329 :
330  interpolationCellPoint<Type>(psi, calcPointField(psi))
331 {}
332 
333 
334 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
Graphite solid properties.
Definition: C.H:48
#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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
const Boundary & boundaryField() const
Return const-reference to the boundary field.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
interpolationCellPointWallModified(const GeometricField< Type, fvPatchField, volMesh > &psi)
Construct from components.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
wordList patchTypes(nPatches)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Generic GeometricField class.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
Generic dimensioned Type class.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
const dimensionedScalar c
Speed of light in a vacuum.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
const dimensionSet & dimensions() const
Return dimensions.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
dynamicFvMesh & mesh
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
dimensionedScalar cos(const dimensionedScalar &ds)
static const Identity< scalar > I
Definition: Identity.H:93
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const GeometricField< Type, fvPatchField, volMesh > &)
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
As interpolationCellPoint, but with the point field modified on wall faces.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
static const zero Zero
Definition: zero.H:97
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
const Field< PointType > & faceNormals() const
Return face normals for patch.
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
const Mesh & mesh() const
Return mesh.
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Internal & ref()
Return a reference to the dimensioned internal field.
label patchi
Given cell centre values and point (vertex) values decompose into tetrahedra and linear interpolate w...
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
Definition: pointFields.H:50
void correctBoundaryConditions()
Correct boundary field.
messageStream Info
label n
label nPoints() const
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92