surfaceSlipDisplacementPointPatchVectorField.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 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 template<>
41 const char*
42 NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>::
43 names[] =
44 {
45  "nearest",
46  "pointNormal",
47  "fixedNormal"
48 };
49 
51  surfaceSlipDisplacementPointPatchVectorField::projectModeNames_;
52 
53 
54 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
55 
56 void surfaceSlipDisplacementPointPatchVectorField::calcProjection
57 (
58  vectorField& displacement
59 ) const
60 {
61  const polyMesh& mesh = patch().boundaryMesh().mesh()();
62  const pointField& localPoints = patch().localPoints();
63  const labelList& meshPoints = patch().meshPoints();
64 
65  // const scalar deltaT = mesh.time().deltaTValue();
66 
67  // Construct large enough vector in direction of projectDir so
68  // we're guaranteed to hit something.
69 
70  //- Per point projection vector:
71  const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
72 
73  // For case of fixed projection vector:
74  vector projectVec(0, 0, 0);
75  if (projectMode_ == FIXEDNORMAL)
76  {
77  vector n = projectDir_/mag(projectDir_);
78  projectVec = projectLen*n;
79  }
80 
81 
82  // Get fixed points (bit of a hack)
83  const pointZone* zonePtr = nullptr;
84 
85  if (frozenPointsZone_.size() > 0)
86  {
87  const meshPointZones& pZones = mesh.pointZones();
88 
89  zonePtr = &pZones[frozenPointsZone_];
90 
91  Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
92  << zonePtr->size() << " points in pointZone " << zonePtr->name()
93  << endl;
94  }
95 
96  // Get the motionSolver from the dynamic mesh
97  const motionSolver& motion =
98  refCast<const dynamicMotionSolverFvMesh>(mesh).motion();
99 
100  // Get the starting locations from the motionSolver
101  const pointField& points0 =
102  refCast<const displacementMotionSolver>(motion).points0();
103 
104  pointField start(meshPoints.size());
105  forAll(start, i)
106  {
107  start[i] = points0[meshPoints[i]] + displacement[i];
108  }
109 
110  label nNotProjected = 0;
111 
112  if (projectMode_ == NEAREST)
113  {
114  List<pointIndexHit> nearest;
115  labelList hitSurfaces;
117  (
118  start,
119  scalarField(start.size(), sqr(projectLen)),
120  hitSurfaces,
121  nearest
122  );
123 
124  forAll(nearest, i)
125  {
126  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
127  {
128  // Fixed point. Reset to point0 location.
129  displacement[i] = points0[meshPoints[i]] - localPoints[i];
130  }
131  else if (nearest[i].hit())
132  {
133  displacement[i] =
134  nearest[i].hitPoint()
135  - points0[meshPoints[i]];
136  }
137  else
138  {
139  nNotProjected++;
140 
141  if (debug)
142  {
143  Pout<< " point:" << meshPoints[i]
144  << " coord:" << localPoints[i]
145  << " did not find any surface within " << projectLen
146  << endl;
147  }
148  }
149  }
150  }
151  else
152  {
153  // Do tests on all points. Combine later on.
154 
155  // 1. Check if already on surface
156  List<pointIndexHit> nearest;
157  {
158  labelList nearestSurface;
160  (
161  start,
162  scalarField(start.size(), sqr(small)),
163  nearestSurface,
164  nearest
165  );
166  }
167 
168  // 2. intersection. (combined later on with information from nearest
169  // above)
170  vectorField projectVecs(start.size(), projectVec);
171 
172  if (projectMode_ == POINTNORMAL)
173  {
174  projectVecs = projectLen*patch().pointNormals();
175  }
176 
177  // Knock out any wedge component
178  scalarField offset(start.size(), 0.0);
179  if (wedgePlane_ >= 0 && wedgePlane_ < vector::nComponents)
180  {
181  forAll(offset, i)
182  {
183  offset[i] = start[i][wedgePlane_];
184  start[i][wedgePlane_] = 0;
185  projectVecs[i][wedgePlane_] = 0;
186  }
187  }
188 
189  List<pointIndexHit> rightHit;
190  {
191  labelList rightSurf;
193  (
194  start,
195  start+projectVecs,
196  rightSurf,
197  rightHit
198  );
199  }
200 
201  List<pointIndexHit> leftHit;
202  {
203  labelList leftSurf;
205  (
206  start,
207  start-projectVecs,
208  leftSurf,
209  leftHit
210  );
211  }
212 
213  // 3. Choose either -fixed, nearest, right, left.
214  forAll(displacement, i)
215  {
216  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
217  {
218  // Fixed point. Reset to point0 location.
219  displacement[i] = points0[meshPoints[i]] - localPoints[i];
220  }
221  else if (nearest[i].hit())
222  {
223  // Found nearest.
224  displacement[i] =
225  nearest[i].hitPoint()
226  - points0[meshPoints[i]];
227  }
228  else
229  {
230  pointIndexHit interPt;
231 
232  if (rightHit[i].hit())
233  {
234  if (leftHit[i].hit())
235  {
236  if
237  (
238  magSqr(rightHit[i].hitPoint()-start[i])
239  < magSqr(leftHit[i].hitPoint()-start[i])
240  )
241  {
242  interPt = rightHit[i];
243  }
244  else
245  {
246  interPt = leftHit[i];
247  }
248  }
249  else
250  {
251  interPt = rightHit[i];
252  }
253  }
254  else
255  {
256  if (leftHit[i].hit())
257  {
258  interPt = leftHit[i];
259  }
260  }
261 
262 
263  if (interPt.hit())
264  {
265  if (wedgePlane_ >= 0 && wedgePlane_ < vector::nComponents)
266  {
267  interPt.rawPoint()[wedgePlane_] += offset[i];
268  }
269  displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
270  }
271  else
272  {
273  nNotProjected++;
274 
275  if (debug)
276  {
277  Pout<< " point:" << meshPoints[i]
278  << " coord:" << localPoints[i]
279  << " did not find any intersection between"
280  << " ray from " << start[i]-projectVecs[i]
281  << " to " << start[i]+projectVecs[i] << endl;
282  }
283  }
284  }
285  }
286  }
287 
288  reduce(nNotProjected, sumOp<label>());
289 
290  if (nNotProjected > 0)
291  {
292  Info<< "surfaceSlipDisplacement :"
293  << " on patch " << patch().name()
294  << " did not project " << nNotProjected
295  << " out of " << returnReduce(localPoints.size(), sumOp<label>())
296  << " points." << endl;
297  }
298 }
299 
300 
301 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
302 
305 (
306  const pointPatch& p,
308 )
309 :
310  pointPatchVectorField(p, iF),
311  projectMode_(NEAREST),
312  projectDir_(Zero),
313  wedgePlane_(-1)
314 {}
315 
316 
319 (
320  const pointPatch& p,
322  const dictionary& dict
323 )
324 :
325  pointPatchVectorField(p, iF, dict),
326  surfacesDict_(dict.subDict("geometry")),
327  projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
328  projectDir_(dict.lookup("projectDirection")),
329  wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
330  frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
331 {}
332 
333 
336 (
338  const pointPatch& p,
340  const pointPatchFieldMapper&
341 )
342 :
343  pointPatchVectorField(p, iF),
344  surfacesDict_(ppf.surfacesDict_),
345  projectMode_(ppf.projectMode_),
346  projectDir_(ppf.projectDir_),
347  wedgePlane_(ppf.wedgePlane_),
348  frozenPointsZone_(ppf.frozenPointsZone_)
349 {}
350 
351 
354 (
357 )
358 :
359  pointPatchVectorField(ppf, iF),
360  surfacesDict_(ppf.surfacesDict_),
361  projectMode_(ppf.projectMode_),
362  projectDir_(ppf.projectDir_),
363  wedgePlane_(ppf.wedgePlane_),
364  frozenPointsZone_(ppf.frozenPointsZone_)
365 {}
366 
367 
368 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
369 
370 const searchableSurfaces&
372 {
373  if (surfacesPtr_.empty())
374  {
375  surfacesPtr_.reset
376  (
378  (
379  IOobject
380  (
381  "abc", // dummy name
382  db().time().constant(),
384  db().time(),
387  ),
388  surfacesDict_,
389  true // use single region naming shortcut
390  )
391  );
392  }
393  return surfacesPtr_();
394 }
395 
396 
398 (
399  const Pstream::commsTypes commsType
400 )
401 {
402  vectorField displacement(this->patchInternalField());
403 
404  // Calculate displacement to project points onto surface
405  calcProjection(displacement);
406 
407  // Get internal field to insert values into
408  Field<vector>& iF = const_cast<Field<vector>&>(this->primitiveField());
409 
410  // setInInternalField(iF, motionU);
411  setInInternalField(iF, displacement);
412 
414 }
415 
416 
418 {
420  writeEntry(os, "geometry", surfacesDict_);
421  writeEntry(os, "projectMode", projectModeNames_[projectMode_]);
422  writeEntry(os, "projectDirection", projectDir_);
423  writeEntry(os, "wedgePlane", wedgePlane_);
424  if (frozenPointsZone_ != word::null)
425  {
426  writeEntry(os, "frozenPointsZone", frozenPointsZone_);
427  }
428 }
429 
430 
431 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
432 
434 (
437 );
438 
439 
440 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
441 
442 } // End namespace Foam
443 
444 // ************************************************************************* //
virtual void evaluate(const Pstream::commsTypes commsType=Pstream::commsTypes::blocking)
Evaluate the patch field.
#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
Displacement follows a triSurface. Use in a displacementMotionSolver as a bc on the pointDisplacement...
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
virtual void evaluate(const Pstream::commsTypes commsType=Pstream::commsTypes::blocking)
Update the patch field.
commsTypes
Types of communications.
Definition: UPstream.H:64
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)
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
Virtual base class for mesh motion solver.
Definition: motionSolver.H:55
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
const Field< vector > & primitiveField() const
Return internal field reference.
Macros for easy insertion into run-time selection tables.
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.
const objectRegistry & db() const
Return local objectRegistry.
static const word & geometryDir()
Return the geometry directory name.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const pointPatch & patch() const
Return patch.
bool hit() const
Is there a hit.
static const word null
An empty word.
Definition: word.H:77
Container for searchableSurfaces.
static const zero Zero
Definition: zero.H:97
virtual const word & name() const =0
Return name.
const word & name() const
Return name.
Definition: zone.H:147
void setInInternalField(Field< Type1 > &iF, const Field< Type1 > &pF, const labelList &meshPoints) const
Given the internal field and a patch field,.
label whichPoint(const label globalPointID) const
Helper function to re-direct to zone::localID(...)
Definition: pointZone.C:126
tmp< Field< vector > > patchInternalField() const
Return field created from appropriate internal field values.
const Point & rawPoint() const
Return point with no checking.
const pointMesh & mesh() const
Return the mesh reference.
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
const searchableSurfaces & surfaces() const
Surface to follow. Demand loads surfaceNames.
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
pointPatchField< vector > pointPatchVectorField
virtual const vectorField & pointNormals() const =0
Return point normals.
const pointBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
Definition: pointPatch.H:133
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual const vectorField & localPoints() const =0
Return mesh points.
virtual void write(Ostream &) const
Write.
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
surfaceSlipDisplacementPointPatchVectorField(const pointPatch &, const DimensionedField< vector, pointMesh > &)
Construct from patch and internal field.
Namespace for OpenFOAM.
virtual const labelList & meshPoints() const =0
Return mesh points.
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