surfaceDisplacementPointPatchVectorField.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-2024 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 
28 #include "Time.H"
29 #include "transformField.H"
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 template<>
42 const char*
44 names[] =
45 {
46  "nearest",
47  "pointNormal",
48  "fixedNormal"
49 };
50 
52  surfaceDisplacementPointPatchVectorField::projectModeNames_;
53 
54 
55 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
56 
57 void surfaceDisplacementPointPatchVectorField::calcProjection
58 (
59  vectorField& displacement
60 ) const
61 {
62  const polyMesh& mesh = patch().boundaryMesh().mesh()();
63  const pointField& localPoints = patch().localPoints();
64  const labelList& meshPoints = patch().meshPoints();
65 
66  // const scalar deltaT = mesh.time().deltaTValue();
67 
68  // Construct large enough vector in direction of projectDir so
69  // we're guaranteed to hit something.
70 
71  //- Per point projection vector:
72  const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
73 
74  // For case of fixed projection vector:
75  vector projectVec(Zero);
76  if (projectMode_ == FIXEDNORMAL)
77  {
78  vector n = projectDir_/mag(projectDir_);
79  projectVec = projectLen*n;
80  }
81 
82 
83  // Get fixed points (bit of a hack)
84  const pointZone* zonePtr = nullptr;
85 
86  if (frozenPointsZone_.size() > 0)
87  {
88  const pointZoneList& pZones = mesh.pointZones();
89 
90  zonePtr = &pZones[frozenPointsZone_];
91 
92  Pout<< "surfaceDisplacementPointPatchVectorField : Fixing all "
93  << zonePtr->size() << " points in pointZone " << zonePtr->name()
94  << endl;
95  }
96 
97  // Get the motionSolver from the mesh
98  const motionSolver& motion =
99  refCast<const fvMeshMovers::motionSolver>
100  (
101  refCast<const fvMesh>(mesh).mover()
102  ).motion();
103 
104  // Get the starting locations from the motionSolver
105  const pointField& points0 =
106  refCast<const displacementMotionSolver>(motion).points0();
107 
108  pointField start(meshPoints.size());
109  forAll(start, i)
110  {
111  start[i] = points0[meshPoints[i]] + displacement[i];
112  }
113 
114  label nNotProjected = 0;
115 
116  if (projectMode_ == NEAREST)
117  {
118  List<pointIndexHit> nearest;
119  labelList hitSurfaces;
121  (
122  start,
123  scalarField(start.size(), sqr(projectLen)),
124  hitSurfaces,
125  nearest
126  );
127 
128  forAll(nearest, i)
129  {
130  if (zonePtr && (zonePtr->localIndex(meshPoints[i]) >= 0))
131  {
132  // Fixed point. Reset to point0 location.
133  displacement[i] = points0[meshPoints[i]] - localPoints[i];
134  }
135  else if (nearest[i].hit())
136  {
137  displacement[i] =
138  nearest[i].hitPoint()
139  - points0[meshPoints[i]];
140  }
141  else
142  {
143  nNotProjected++;
144 
145  if (debug)
146  {
147  Pout<< " point:" << meshPoints[i]
148  << " coord:" << localPoints[i]
149  << " did not find any surface within " << projectLen
150  << endl;
151  }
152  }
153  }
154  }
155  else
156  {
157  // Do tests on all points. Combine later on.
158 
159  // 1. Check if already on surface
160  List<pointIndexHit> nearest;
161  {
162  labelList nearestSurface;
164  (
165  start,
166  scalarField(start.size(), sqr(small)),
167  nearestSurface,
168  nearest
169  );
170  }
171 
172  // 2. intersection. (combined later on with information from nearest
173  // above)
174  vectorField projectVecs(start.size(), projectVec);
175 
176  if (projectMode_ == POINTNORMAL)
177  {
178  projectVecs = projectLen*patch().pointNormals();
179  }
180 
181  // Knock out any wedge component
182  scalarField offset(start.size(), 0.0);
183  if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
184  {
185  forAll(offset, i)
186  {
187  offset[i] = start[i][wedgePlane_];
188  start[i][wedgePlane_] = 0;
189  projectVecs[i][wedgePlane_] = 0;
190  }
191  }
192 
193  List<pointIndexHit> rightHit;
194  {
195  labelList rightSurf;
197  (
198  start,
199  start+projectVecs,
200  rightSurf,
201  rightHit
202  );
203  }
204 
205  List<pointIndexHit> leftHit;
206  {
207  labelList leftSurf;
209  (
210  start,
211  start-projectVecs,
212  leftSurf,
213  leftHit
214  );
215  }
216 
217  // 3. Choose either -fixed, nearest, right, left.
218  forAll(displacement, i)
219  {
220  if (zonePtr && (zonePtr->localIndex(meshPoints[i]) >= 0))
221  {
222  // Fixed point. Reset to point0 location.
223  displacement[i] = points0[meshPoints[i]] - localPoints[i];
224  }
225  else if (nearest[i].hit())
226  {
227  // Found nearest.
228  displacement[i] =
229  nearest[i].hitPoint()
230  - points0[meshPoints[i]];
231  }
232  else
233  {
234  pointIndexHit interPt;
235 
236  if (rightHit[i].hit())
237  {
238  if (leftHit[i].hit())
239  {
240  if
241  (
242  magSqr(rightHit[i].hitPoint()-start[i])
243  < magSqr(leftHit[i].hitPoint()-start[i])
244  )
245  {
246  interPt = rightHit[i];
247  }
248  else
249  {
250  interPt = leftHit[i];
251  }
252  }
253  else
254  {
255  interPt = rightHit[i];
256  }
257  }
258  else
259  {
260  if (leftHit[i].hit())
261  {
262  interPt = leftHit[i];
263  }
264  }
265 
266 
267  if (interPt.hit())
268  {
269  if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
270  {
271  interPt.rawPoint()[wedgePlane_] += offset[i];
272  }
273  displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
274  }
275  else
276  {
277  nNotProjected++;
278 
279  if (debug)
280  {
281  Pout<< " point:" << meshPoints[i]
282  << " coord:" << localPoints[i]
283  << " did not find any intersection between"
284  << " ray from " << start[i]-projectVecs[i]
285  << " to " << start[i]+projectVecs[i] << endl;
286  }
287  }
288  }
289  }
290  }
291 
292  reduce(nNotProjected, sumOp<label>());
293 
294  if (nNotProjected > 0)
295  {
296  Info<< "surfaceDisplacement :"
297  << " on patch " << patch().name()
298  << " did not project " << nNotProjected
299  << " out of " << returnReduce(localPoints.size(), sumOp<label>())
300  << " points." << endl;
301  }
302 }
303 
304 
305 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
306 
309 (
310  const pointPatch& p,
312  const dictionary& dict
313 )
314 :
315  fixedValuePointPatchVectorField(p, iF, dict),
316  velocity_(dict.lookup<vector>("velocity", dimVelocity)),
317  surfacesDict_(dict.subDict("geometry")),
318  projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
319  projectDir_(dict.lookup<vector>("projectDirection", dimless)),
320  wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
321  frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
322 {
323  if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
324  {
326  << "All components of velocity have to be positive : "
327  << velocity_ << nl
328  << "Set velocity components to a great value if no clipping"
329  << " necessary." << exit(FatalError);
330  }
331 }
332 
333 
336 (
338  const pointPatch& p,
340  const fieldMapper& mapper
341 )
342 :
343  fixedValuePointPatchVectorField(ppf, p, iF, mapper),
344  velocity_(ppf.velocity_),
345  surfacesDict_(ppf.surfacesDict_),
346  projectMode_(ppf.projectMode_),
347  projectDir_(ppf.projectDir_),
348  wedgePlane_(ppf.wedgePlane_),
349  frozenPointsZone_(ppf.frozenPointsZone_)
350 {}
351 
352 
355 (
358 )
359 :
360  fixedValuePointPatchVectorField(ppf, iF),
361  velocity_(ppf.velocity_),
362  surfacesDict_(ppf.surfacesDict_),
363  projectMode_(ppf.projectMode_),
364  projectDir_(ppf.projectDir_),
365  wedgePlane_(ppf.wedgePlane_),
366  frozenPointsZone_(ppf.frozenPointsZone_)
367 {}
368 
369 
370 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
371 
372 const searchableSurfaces&
374 {
375  if (surfacesPtr_.empty())
376  {
377  surfacesPtr_.reset
378  (
380  (
381  IOobject
382  (
383  "abc", // dummy name
384  db().time().constant(),
385  searchableSurface::geometryDir(db().time()),
386  db().time(),
389  ),
390  surfacesDict_,
391  true // allow single-region shortcut
392  )
393  );
394  }
395  return surfacesPtr_();
396 }
397 
398 
400 {
401  if (this->updated())
402  {
403  return;
404  }
405 
406  const polyMesh& mesh = patch().boundaryMesh().mesh()();
407 
408  vectorField currentDisplacement(this->patchInternalField());
409 
410  // Calculate intersections with surface w.r.t points0.
411  vectorField displacement(currentDisplacement);
412  calcProjection(displacement);
413 
414  // offset wrt current displacement
415  vectorField offset(displacement-currentDisplacement);
416 
417  // Clip offset to maximum displacement possible: velocity*timestep
418 
419  const scalar deltaT = mesh.time().deltaTValue();
420  const vector clipVelocity = velocity_*deltaT;
421 
422  forAll(displacement, i)
423  {
424  vector& d = offset[i];
425 
426  for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
427  {
428  if (d[cmpt] < 0)
429  {
430  d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
431  }
432  else
433  {
434  d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
435  }
436  }
437  }
438 
439  this->operator==(currentDisplacement+offset);
440 
441  fixedValuePointPatchVectorField::updateCoeffs();
442 }
443 
444 
446 {
448  writeEntry(os, "velocity", velocity_);
449  writeEntry(os, "geometry", surfacesDict_);
450  writeEntry(os, "projectMode", projectModeNames_[projectMode_]);
451  writeEntry(os, "projectDirection", projectDir_);
452  writeEntry(os, "wedgePlane", wedgePlane_);
453  if (frozenPointsZone_ != word::null)
454  {
455  writeEntry(os, "frozenPointsZone", frozenPointsZone_);
456  }
457 }
458 
459 
460 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
461 
463 (
464  fixedValuePointPatchVectorField,
466 );
467 
468 
469 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
470 
471 } // End namespace Foam
472 
473 // ************************************************************************* //
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:54
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:85
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:34
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:105
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:60
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:66
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Abstract base class for field mapping.
Definition: fieldMapper.H:48
const Time & time() const
Return time.
Basic pointPatch represents a set of points from the mesh.
Definition: pointPatch.H:57
const polyMesh & mesh() const
Return the mesh reference.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const pointZoneList & pointZones() const
Return point zones.
Definition: polyMesh.H:440
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:410
static const word & geometryDir()
Return the geometry directory name.
Container for searchableSurfaces.
void findAnyIntersection(const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &) const
Find any intersection. Return hit point information and.
void findNearest(const pointField &, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest. Return -1 (and a miss()) or surface and nearest.
Displacement fixed by projection onto triSurface. Use in a displacementMotionSolver as a bc on the po...
const searchableSurfaces & surfaces() const
Surface to follow. Demand loads surfaceNames.
surfaceDisplacementPointPatchVectorField(const pointPatch &, const DimensionedField< vector, pointMesh > &, const dictionary &)
Construct from patch, internal field and dictionary.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
A class for handling words, derived from string.
Definition: word.H:62
static const word null
An empty word.
Definition: word.H:77
IOporosityModelList pZones(mesh)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
makePointPatchTypeField(pointPatchVectorField, angularOscillatingDisplacementPointPatchVectorField)
bool read(const char *, int32_t &)
Definition: int32IO.C:85
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
const dimensionSet dimless
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Field< vector > vectorField
Specialisation of Field<T> for vector.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const dimensionSet dimVelocity
error FatalError
void offset(label &lst, const label o)
static const char nl
Definition: Ostream.H:266
dimensioned< scalar > magSqr(const dimensioned< Type > &)
uint8_t direction
Definition: direction.H:45
dictionary dict
volScalarField & p
Spatial transformation functions for primitive fields.