surfaceDisplacementPointPatchVectorField.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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"
30 #include "fvMesh.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 = NULL;
85 
86  if (frozenPointsZone_.size() > 0)
87  {
88  const pointZoneMesh& 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 starting locations from the motionSolver
98  const pointField& points0 = mesh.lookupObject<displacementMotionSolver>
99  (
100  "dynamicMeshDict"
101  ).points0();
102 
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<< "surfaceDisplacement :"
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  fixedValuePointPatchVectorField(p, iF),
311  velocity_(Zero),
312  projectMode_(NEAREST),
313  projectDir_(Zero),
314  wedgePlane_(-1)
315 {}
316 
317 
320 (
321  const pointPatch& p,
323  const dictionary& dict
324 )
325 :
326  fixedValuePointPatchVectorField(p, iF, dict),
327  velocity_(dict.lookup("velocity")),
328  surfacesDict_(dict.subDict("geometry")),
329  projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
330  projectDir_(dict.lookup("projectDirection")),
331  wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
332  frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
333 {
334  if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
335  {
337  << "All components of velocity have to be positive : "
338  << velocity_ << nl
339  << "Set velocity components to a great value if no clipping"
340  << " necessary." << exit(FatalError);
341  }
342 }
343 
344 
347 (
349  const pointPatch& p,
351  const pointPatchFieldMapper& mapper
352 )
353 :
354  fixedValuePointPatchVectorField(ppf, p, iF, mapper),
355  velocity_(ppf.velocity_),
356  surfacesDict_(ppf.surfacesDict_),
357  projectMode_(ppf.projectMode_),
358  projectDir_(ppf.projectDir_),
359  wedgePlane_(ppf.wedgePlane_),
360  frozenPointsZone_(ppf.frozenPointsZone_)
361 {}
362 
363 
366 (
368 )
369 :
370  fixedValuePointPatchVectorField(ppf),
371  velocity_(ppf.velocity_),
372  surfacesDict_(ppf.surfacesDict_),
373  projectMode_(ppf.projectMode_),
374  projectDir_(ppf.projectDir_),
375  wedgePlane_(ppf.wedgePlane_),
376  frozenPointsZone_(ppf.frozenPointsZone_)
377 {}
378 
379 
382 (
385 )
386 :
387  fixedValuePointPatchVectorField(ppf, iF),
388  velocity_(ppf.velocity_),
389  surfacesDict_(ppf.surfacesDict_),
390  projectMode_(ppf.projectMode_),
391  projectDir_(ppf.projectDir_),
392  wedgePlane_(ppf.wedgePlane_),
393  frozenPointsZone_(ppf.frozenPointsZone_)
394 {}
395 
396 
397 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
398 
399 const searchableSurfaces&
401 {
402  if (surfacesPtr_.empty())
403  {
404  surfacesPtr_.reset
405  (
407  (
408  IOobject
409  (
410  "abc", // dummy name
411  db().time().constant(), // directory
412  "triSurface", // instance
413  db().time(), // registry
416  ),
417  surfacesDict_,
418  true // allow single-region shortcut
419  )
420  );
421  }
422  return surfacesPtr_();
423 }
424 
425 
427 {
428  if (this->updated())
429  {
430  return;
431  }
432 
433  const polyMesh& mesh = patch().boundaryMesh().mesh()();
434 
435  vectorField currentDisplacement(this->patchInternalField());
436 
437  // Calculate intersections with surface w.r.t points0.
438  vectorField displacement(currentDisplacement);
439  calcProjection(displacement);
440 
441  // offset wrt current displacement
442  vectorField offset(displacement-currentDisplacement);
443 
444  // Clip offset to maximum displacement possible: velocity*timestep
445 
446  const scalar deltaT = mesh.time().deltaTValue();
447  const vector clipVelocity = velocity_*deltaT;
448 
449  forAll(displacement, i)
450  {
451  vector& d = offset[i];
452 
453  for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
454  {
455  if (d[cmpt] < 0)
456  {
457  d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
458  }
459  else
460  {
461  d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
462  }
463  }
464  }
465 
466  this->operator==(currentDisplacement+offset);
467 
468  fixedValuePointPatchVectorField::updateCoeffs();
469 }
470 
471 
473 {
475  os.writeKeyword("velocity") << velocity_
476  << token::END_STATEMENT << nl;
477  os.writeKeyword("geometry") << surfacesDict_
478  << token::END_STATEMENT << nl;
479  os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
480  << token::END_STATEMENT << nl;
481  os.writeKeyword("projectDirection") << projectDir_
482  << token::END_STATEMENT << nl;
483  os.writeKeyword("wedgePlane") << wedgePlane_
484  << token::END_STATEMENT << nl;
485  if (frozenPointsZone_ != word::null)
486  {
487  os.writeKeyword("frozenPointsZone") << frozenPointsZone_
488  << token::END_STATEMENT << nl;
489  }
490 }
491 
492 
493 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
494 
496 (
497  fixedValuePointPatchVectorField,
499 );
500 
501 
502 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
503 
504 } // End namespace Foam
505 
506 // ************************************************************************* //
const Time & time() const
Return time.
const Cmpt & z() const
Definition: VectorI.H:87
Virtual base class for displacement motion solver.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
uint8_t direction
Definition: direction.H:46
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:427
const Cmpt & x() const
Definition: VectorI.H:75
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
bool hit() const
Is there a hit.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const searchableSurfaces & surfaces() const
Surface to follow. Demand loads surfaceNames.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:633
Foam::pointPatchFieldMapper.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
surfaceDisplacementPointPatchVectorField(const pointPatch &, const DimensionedField< vector, pointMesh > &)
Construct from patch and internal field.
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:52
Macros for easy insertion into run-time selection tables.
Spatial transformation functions for primitive fields.
const word & name() const
Return name.
Definition: zone.H:150
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:96
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.
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Cmpt & y() const
Definition: VectorI.H:81
static const word null
An empty word.
Definition: word.H:77
Container for searchableSurfaces.
static const zero Zero
Definition: zero.H:91
const Point & rawPoint() const
Return point with no checking.
Displacement fixed by projection onto triSurface. Use in a displacementMotionSolver as a bc on the po...
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
void findAnyIntersection(const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &) const
Find any intersection. Return hit point information and.
static const char nl
Definition: Ostream.H:262
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:54
label whichPoint(const label globalPointID) const
Helper function to re-direct to zone::localID(...)
Definition: pointZone.C:126
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:41
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:457
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
A subset of mesh points. The labels of points in the zone can be obtained from the addressing() list...
Definition: pointZone.H:62
Constant dispersed-phase particle diameter model.
Basic pointPatch represents a set of points from the mesh.
Definition: pointPatch.H:56
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
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
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Enum read(Istream &) const
Read a word from Istream and return the corresponding.
Definition: NamedEnum.C:61
runTime write()
Namespace for OpenFOAM.
IOporosityModelList pZones(mesh)
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451