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-2021 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*
43 NamedEnum<surfaceDisplacementPointPatchVectorField::projectMode, 3>::
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 meshPointZones& 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 dynamic mesh
98  const motionSolver& motion =
99  refCast<const dynamicMotionSolverFvMesh>(mesh).motion();
100 
101  // Get the starting locations from the motionSolver
102  const pointField& points0 =
103  refCast<const displacementMotionSolver>(motion).points0();
104 
105  pointField start(meshPoints.size());
106  forAll(start, i)
107  {
108  start[i] = points0[meshPoints[i]] + displacement[i];
109  }
110 
111  label nNotProjected = 0;
112 
113  if (projectMode_ == NEAREST)
114  {
115  List<pointIndexHit> nearest;
116  labelList hitSurfaces;
118  (
119  start,
120  scalarField(start.size(), sqr(projectLen)),
121  hitSurfaces,
122  nearest
123  );
124 
125  forAll(nearest, i)
126  {
127  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
128  {
129  // Fixed point. Reset to point0 location.
130  displacement[i] = points0[meshPoints[i]] - localPoints[i];
131  }
132  else if (nearest[i].hit())
133  {
134  displacement[i] =
135  nearest[i].hitPoint()
136  - points0[meshPoints[i]];
137  }
138  else
139  {
140  nNotProjected++;
141 
142  if (debug)
143  {
144  Pout<< " point:" << meshPoints[i]
145  << " coord:" << localPoints[i]
146  << " did not find any surface within " << projectLen
147  << endl;
148  }
149  }
150  }
151  }
152  else
153  {
154  // Do tests on all points. Combine later on.
155 
156  // 1. Check if already on surface
157  List<pointIndexHit> nearest;
158  {
159  labelList nearestSurface;
161  (
162  start,
163  scalarField(start.size(), sqr(small)),
164  nearestSurface,
165  nearest
166  );
167  }
168 
169  // 2. intersection. (combined later on with information from nearest
170  // above)
171  vectorField projectVecs(start.size(), projectVec);
172 
173  if (projectMode_ == POINTNORMAL)
174  {
175  projectVecs = projectLen*patch().pointNormals();
176  }
177 
178  // Knock out any wedge component
179  scalarField offset(start.size(), 0.0);
180  if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
181  {
182  forAll(offset, i)
183  {
184  offset[i] = start[i][wedgePlane_];
185  start[i][wedgePlane_] = 0;
186  projectVecs[i][wedgePlane_] = 0;
187  }
188  }
189 
190  List<pointIndexHit> rightHit;
191  {
192  labelList rightSurf;
194  (
195  start,
196  start+projectVecs,
197  rightSurf,
198  rightHit
199  );
200  }
201 
202  List<pointIndexHit> leftHit;
203  {
204  labelList leftSurf;
206  (
207  start,
208  start-projectVecs,
209  leftSurf,
210  leftHit
211  );
212  }
213 
214  // 3. Choose either -fixed, nearest, right, left.
215  forAll(displacement, i)
216  {
217  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
218  {
219  // Fixed point. Reset to point0 location.
220  displacement[i] = points0[meshPoints[i]] - localPoints[i];
221  }
222  else if (nearest[i].hit())
223  {
224  // Found nearest.
225  displacement[i] =
226  nearest[i].hitPoint()
227  - points0[meshPoints[i]];
228  }
229  else
230  {
231  pointIndexHit interPt;
232 
233  if (rightHit[i].hit())
234  {
235  if (leftHit[i].hit())
236  {
237  if
238  (
239  magSqr(rightHit[i].hitPoint()-start[i])
240  < magSqr(leftHit[i].hitPoint()-start[i])
241  )
242  {
243  interPt = rightHit[i];
244  }
245  else
246  {
247  interPt = leftHit[i];
248  }
249  }
250  else
251  {
252  interPt = rightHit[i];
253  }
254  }
255  else
256  {
257  if (leftHit[i].hit())
258  {
259  interPt = leftHit[i];
260  }
261  }
262 
263 
264  if (interPt.hit())
265  {
266  if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
267  {
268  interPt.rawPoint()[wedgePlane_] += offset[i];
269  }
270  displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
271  }
272  else
273  {
274  nNotProjected++;
275 
276  if (debug)
277  {
278  Pout<< " point:" << meshPoints[i]
279  << " coord:" << localPoints[i]
280  << " did not find any intersection between"
281  << " ray from " << start[i]-projectVecs[i]
282  << " to " << start[i]+projectVecs[i] << endl;
283  }
284  }
285  }
286  }
287  }
288 
289  reduce(nNotProjected, sumOp<label>());
290 
291  if (nNotProjected > 0)
292  {
293  Info<< "surfaceDisplacement :"
294  << " on patch " << patch().name()
295  << " did not project " << nNotProjected
296  << " out of " << returnReduce(localPoints.size(), sumOp<label>())
297  << " points." << endl;
298  }
299 }
300 
301 
302 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
303 
306 (
307  const pointPatch& p,
309 )
310 :
311  fixedValuePointPatchVectorField(p, iF),
312  velocity_(Zero),
313  projectMode_(NEAREST),
314  projectDir_(Zero),
315  wedgePlane_(-1)
316 {}
317 
318 
321 (
322  const pointPatch& p,
324  const dictionary& dict
325 )
326 :
327  fixedValuePointPatchVectorField(p, iF, dict),
328  velocity_(dict.lookup("velocity")),
329  surfacesDict_(dict.subDict("geometry")),
330  projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
331  projectDir_(dict.lookup("projectDirection")),
332  wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
333  frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
334 {
335  if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
336  {
338  << "All components of velocity have to be positive : "
339  << velocity_ << nl
340  << "Set velocity components to a great value if no clipping"
341  << " necessary." << exit(FatalError);
342  }
343 }
344 
345 
348 (
350  const pointPatch& p,
352  const pointPatchFieldMapper& mapper
353 )
354 :
355  fixedValuePointPatchVectorField(ppf, p, iF, mapper),
356  velocity_(ppf.velocity_),
357  surfacesDict_(ppf.surfacesDict_),
358  projectMode_(ppf.projectMode_),
359  projectDir_(ppf.projectDir_),
360  wedgePlane_(ppf.wedgePlane_),
361  frozenPointsZone_(ppf.frozenPointsZone_)
362 {}
363 
364 
367 (
370 )
371 :
372  fixedValuePointPatchVectorField(ppf, iF),
373  velocity_(ppf.velocity_),
374  surfacesDict_(ppf.surfacesDict_),
375  projectMode_(ppf.projectMode_),
376  projectDir_(ppf.projectDir_),
377  wedgePlane_(ppf.wedgePlane_),
378  frozenPointsZone_(ppf.frozenPointsZone_)
379 {}
380 
381 
382 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
383 
384 const searchableSurfaces&
386 {
387  if (surfacesPtr_.empty())
388  {
389  surfacesPtr_.reset
390  (
392  (
393  IOobject
394  (
395  "abc", // dummy name
396  db().time().constant(),
397  searchableSurface::geometryDir(db().time()),
398  db().time(),
401  ),
402  surfacesDict_,
403  true // allow single-region shortcut
404  )
405  );
406  }
407  return surfacesPtr_();
408 }
409 
410 
412 {
413  if (this->updated())
414  {
415  return;
416  }
417 
418  const polyMesh& mesh = patch().boundaryMesh().mesh()();
419 
420  vectorField currentDisplacement(this->patchInternalField());
421 
422  // Calculate intersections with surface w.r.t points0.
423  vectorField displacement(currentDisplacement);
424  calcProjection(displacement);
425 
426  // offset wrt current displacement
427  vectorField offset(displacement-currentDisplacement);
428 
429  // Clip offset to maximum displacement possible: velocity*timestep
430 
431  const scalar deltaT = mesh.time().deltaTValue();
432  const vector clipVelocity = velocity_*deltaT;
433 
434  forAll(displacement, i)
435  {
436  vector& d = offset[i];
437 
438  for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
439  {
440  if (d[cmpt] < 0)
441  {
442  d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
443  }
444  else
445  {
446  d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
447  }
448  }
449  }
450 
451  this->operator==(currentDisplacement+offset);
452 
453  fixedValuePointPatchVectorField::updateCoeffs();
454 }
455 
456 
458 {
460  writeEntry(os, "velocity", velocity_);
461  writeEntry(os, "geometry", surfacesDict_);
462  writeEntry(os, "projectMode", projectModeNames_[projectMode_]);
463  writeEntry(os, "projectDirection", projectDir_);
464  writeEntry(os, "wedgePlane", wedgePlane_);
465  if (frozenPointsZone_ != word::null)
466  {
467  writeEntry(os, "frozenPointsZone", frozenPointsZone_);
468  }
469 }
470 
471 
472 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
473 
475 (
476  fixedValuePointPatchVectorField,
478 );
479 
480 
481 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
482 
483 } // End namespace Foam
484 
485 // ************************************************************************* //
#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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
void findAnyIntersection(const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &) const
Find any intersection. Return hit point information and.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
uint8_t direction
Definition: direction.H:45
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Foam::pointPatchFieldMapper.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
surfaceDisplacementPointPatchVectorField(const pointPatch &, const DimensionedField< vector, pointMesh > &)
Construct from patch and internal field.
Virtual base class for mesh motion solver.
Definition: motionSolver.H:55
const Cmpt & z() const
Definition: VectorI.H:87
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
Macros for easy insertion into run-time selection tables.
const Cmpt & y() const
Definition: VectorI.H:81
const searchableSurfaces & surfaces() const
Surface to follow. Demand loads surfaceNames.
const meshPointZones & pointZones() const
Return point zones.
Definition: polyMesh.H:470
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:982
Spatial transformation functions for primitive fields.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:99
dynamicFvMesh & mesh
void findNearest(const pointField &, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest. Return -1 (and a miss()) or surface and nearest.
static const word & geometryDir()
Return the geometry directory name.
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:41
bool hit() const
Is there a hit.
static const word null
An empty word.
Definition: word.H:77
Container for searchableSurfaces.
void write(std::ostream &os, const bool binary, List< floatScalar > &fField)
Write floats ascii or binary.
static const zero Zero
Definition: zero.H:97
const word & name() const
Return name.
Definition: zone.H:147
label whichPoint(const label globalPointID) const
Helper function to re-direct to zone::localID(...)
Definition: pointZone.C:126
const Cmpt & x() const
Definition: VectorI.H:75
Displacement fixed by projection onto triSurface. Use in a displacementMotionSolver as a bc on the po...
const Point & rawPoint() const
Return point with no checking.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
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)
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:440
A subset of mesh points. The labels of points in the zone can be obtained from the addressing() list...
Definition: pointZone.H:62
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Basic pointPatch represents a set of points from the mesh.
Definition: pointPatch.H:56
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
messageStream Info
makePointPatchTypeField(pointPatchVectorField, solidBodyMotionDisplacementPointPatchVectorField)
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Enum read(Istream &) const
Read a word from Istream and return the corresponding.
Definition: NamedEnum.C:61
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.
IOporosityModelList pZones(mesh)
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
virtual void updateCoeffs()
Update the coefficients associated with the patch field.