FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2007 OpenCFD Ltd.
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 <OpenFOAM/Time.H>
30 #include <finiteVolume/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().deltaT().value();
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;
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[pZones.findZoneID(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<displacementFvMotionSolver>
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_(vector::zero),
312  projectMode_(NEAREST),
313  projectDir_(vector::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  (
338  "surfaceDisplacementPointPatchVectorField::\n"
339  "surfaceDisplacementPointPatchVectorField\n"
340  "(\n"
341  " const pointPatch& p,\n"
342  " const DimensionedField<vector, pointMesh>& iF,\n"
343  " const dictionary& dict\n"
344  ")\n"
345  ) << "All components of velocity have to be positive : "
346  << velocity_ << nl
347  << "Set velocity components to a great value if no clipping"
348  << " necessary." << exit(FatalError);
349  }
350 }
351 
352 
355 (
357  const pointPatch& p,
359  const pointPatchFieldMapper& mapper
360 )
361 :
362  fixedValuePointPatchVectorField(ppf, p, iF, mapper),
363  velocity_(ppf.velocity_),
364  surfacesDict_(ppf.surfacesDict_),
365  projectMode_(ppf.projectMode_),
366  projectDir_(ppf.projectDir_),
367  wedgePlane_(ppf.wedgePlane_),
368  frozenPointsZone_(ppf.frozenPointsZone_)
369 {}
370 
371 
374 (
376 )
377 :
378  fixedValuePointPatchVectorField(ppf),
379  velocity_(ppf.velocity_),
380  surfacesDict_(ppf.surfacesDict_),
381  projectMode_(ppf.projectMode_),
382  projectDir_(ppf.projectDir_),
383  wedgePlane_(ppf.wedgePlane_),
384  frozenPointsZone_(ppf.frozenPointsZone_)
385 {}
386 
387 
390 (
393 )
394 :
395  fixedValuePointPatchVectorField(ppf, iF),
396  velocity_(ppf.velocity_),
397  surfacesDict_(ppf.surfacesDict_),
398  projectMode_(ppf.projectMode_),
399  projectDir_(ppf.projectDir_),
400  wedgePlane_(ppf.wedgePlane_),
401  frozenPointsZone_(ppf.frozenPointsZone_)
402 {}
403 
404 
405 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
406 
407 const searchableSurfaces&
409 {
410  if (surfacesPtr_.empty())
411  {
412  surfacesPtr_.reset
413  (
415  (
416  IOobject
417  (
418  "abc", // dummy name
419  db().time().constant(), // directory
420  "triSurface", // instance
421  db().time(), // registry
424  ),
425  surfacesDict_
426  )
427  );
428  }
429  return surfacesPtr_();
430 }
431 
432 
434 {
435  if (this->updated())
436  {
437  return;
438  }
439 
440  const polyMesh& mesh = patch().boundaryMesh().mesh()();
441 
442  vectorField currentDisplacement = this->patchInternalField();
443 
444  // Calculate intersections with surface w.r.t points0.
445  vectorField displacement(currentDisplacement);
446  calcProjection(displacement);
447 
448  // offset wrt current displacement
449  vectorField offset = displacement-currentDisplacement;
450 
451  // Clip offset to maximum displacement possible: velocity*timestep
452 
453  const scalar deltaT = mesh.time().deltaT().value();
454  const vector clipVelocity = velocity_*deltaT;
455 
456  forAll(displacement, i)
457  {
458  vector& d = offset[i];
459 
460  for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
461  {
462  if (d[cmpt] < 0)
463  {
464  d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
465  }
466  else
467  {
468  d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
469  }
470  }
471  }
472 
473  this->operator==(currentDisplacement+offset);
474 
476 }
477 
478 
480 {
482  os.writeKeyword("velocity") << velocity_
483  << token::END_STATEMENT << nl;
484  os.writeKeyword("geometry") << surfacesDict_
485  << token::END_STATEMENT << nl;
486  os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
487  << token::END_STATEMENT << nl;
488  os.writeKeyword("projectDirection") << projectDir_
489  << token::END_STATEMENT << nl;
490  os.writeKeyword("wedgePlane") << wedgePlane_
491  << token::END_STATEMENT << nl;
492  if (frozenPointsZone_ != word::null)
493  {
494  os.writeKeyword("frozenPointsZone") << frozenPointsZone_
495  << token::END_STATEMENT << nl;
496  }
497 }
498 
499 
500 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
501 
503 (
504  fixedValuePointPatchVectorField,
506 );
507 
508 
509 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
510 
511 } // End namespace Foam
512 
513 // ************************ vim: set sw=4 sts=4 et: ************************ //