FreeFOAM The Cross-Platform CFD Toolkit
surfaceSlipDisplacementPointPatchVectorField.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-2010 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 // * * * * * * * * * * * * * * 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().deltaT().value();
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;
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 = NULL;
84 
85  if (frozenPointsZone_.size() > 0)
86  {
87  const pointZoneMesh& pZones = mesh.pointZones();
88  label zoneID = pZones.findZoneID(frozenPointsZone_);
89  if (zoneID == -1)
90  {
92  (
93  "surfaceSlipDisplacementPointPatchVectorField::calcProjection()"
94  ) << "Cannot find zone " << frozenPointsZone_ << endl
95  << "Valid zones are " << pZones.name() << exit(FatalError);
96  }
97 
98  zonePtr = &pZones[zoneID];
99 
100  Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
101  << zonePtr->size() << " points in pointZone " << zonePtr->name()
102  << endl;
103  }
104 
105  // Get the starting locations from the motionSolver
106  const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
107  (
108  "dynamicMeshDict"
109  ).points0();
110 
111 
112  pointField start(meshPoints.size());
113  forAll(start, i)
114  {
115  start[i] = points0[meshPoints[i]] + displacement[i];
116  }
117 
118  label nNotProjected = 0;
119 
120  if (projectMode_ == NEAREST)
121  {
122  List<pointIndexHit> nearest;
123  labelList hitSurfaces;
125  (
126  start,
127  scalarField(start.size(), sqr(projectLen)),
128  hitSurfaces,
129  nearest
130  );
131 
132  forAll(nearest, i)
133  {
134  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
135  {
136  // Fixed point. Reset to point0 location.
137  displacement[i] = points0[meshPoints[i]] - localPoints[i];
138  }
139  else if (nearest[i].hit())
140  {
141  displacement[i] =
142  nearest[i].hitPoint()
143  - points0[meshPoints[i]];
144  }
145  else
146  {
147  nNotProjected++;
148 
149  if (debug)
150  {
151  Pout<< " point:" << meshPoints[i]
152  << " coord:" << localPoints[i]
153  << " did not find any surface within " << projectLen
154  << endl;
155  }
156  }
157  }
158  }
159  else
160  {
161  // Do tests on all points. Combine later on.
162 
163  // 1. Check if already on surface
164  List<pointIndexHit> nearest;
165  {
166  labelList nearestSurface;
168  (
169  start,
170  scalarField(start.size(), sqr(SMALL)),
171  nearestSurface,
172  nearest
173  );
174  }
175 
176  // 2. intersection. (combined later on with information from nearest
177  // above)
178  vectorField projectVecs(start.size(), projectVec);
179 
180  if (projectMode_ == POINTNORMAL)
181  {
182  projectVecs = projectLen*patch().pointNormals();
183  }
184 
185  // Knock out any wedge component
186  scalarField offset(start.size(), 0.0);
187  if (wedgePlane_ >= 0 && wedgePlane_ < vector::nComponents)
188  {
189  forAll(offset, i)
190  {
191  offset[i] = start[i][wedgePlane_];
192  start[i][wedgePlane_] = 0;
193  projectVecs[i][wedgePlane_] = 0;
194  }
195  }
196 
197  List<pointIndexHit> rightHit;
198  {
199  labelList rightSurf;
201  (
202  start,
203  start+projectVecs,
204  rightSurf,
205  rightHit
206  );
207  }
208 
209  List<pointIndexHit> leftHit;
210  {
211  labelList leftSurf;
213  (
214  start,
215  start-projectVecs,
216  leftSurf,
217  leftHit
218  );
219  }
220 
221  // 3. Choose either -fixed, nearest, right, left.
222  forAll(displacement, i)
223  {
224  if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
225  {
226  // Fixed point. Reset to point0 location.
227  displacement[i] = points0[meshPoints[i]] - localPoints[i];
228  }
229  else if (nearest[i].hit())
230  {
231  // Found nearest.
232  displacement[i] =
233  nearest[i].hitPoint()
234  - points0[meshPoints[i]];
235  }
236  else
237  {
238  pointIndexHit interPt;
239 
240  if (rightHit[i].hit())
241  {
242  if (leftHit[i].hit())
243  {
244  if
245  (
246  magSqr(rightHit[i].hitPoint()-start[i])
247  < magSqr(leftHit[i].hitPoint()-start[i])
248  )
249  {
250  interPt = rightHit[i];
251  }
252  else
253  {
254  interPt = leftHit[i];
255  }
256  }
257  else
258  {
259  interPt = rightHit[i];
260  }
261  }
262  else
263  {
264  if (leftHit[i].hit())
265  {
266  interPt = leftHit[i];
267  }
268  }
269 
270 
271  if (interPt.hit())
272  {
273  if (wedgePlane_ >= 0 && wedgePlane_ < vector::nComponents)
274  {
275  interPt.rawPoint()[wedgePlane_] += offset[i];
276  }
277  displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
278  }
279  else
280  {
281  nNotProjected++;
282 
283  if (debug)
284  {
285  Pout<< " point:" << meshPoints[i]
286  << " coord:" << localPoints[i]
287  << " did not find any intersection between"
288  << " ray from " << start[i]-projectVecs[i]
289  << " to " << start[i]+projectVecs[i] << endl;
290  }
291  }
292  }
293  }
294  }
295 
296  reduce(nNotProjected, sumOp<label>());
297 
298  if (nNotProjected > 0)
299  {
300  Info<< "surfaceSlipDisplacement :"
301  << " on patch " << patch().name()
302  << " did not project " << nNotProjected
303  << " out of " << returnReduce(localPoints.size(), sumOp<label>())
304  << " points." << endl;
305  }
306 }
307 
308 
309 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
310 
313 (
314  const pointPatch& p,
316 )
317 :
318  pointPatchVectorField(p, iF),
319  projectMode_(NEAREST),
320  projectDir_(vector::zero),
321  wedgePlane_(-1)
322 {}
323 
324 
327 (
328  const pointPatch& p,
330  const dictionary& dict
331 )
332 :
333  pointPatchVectorField(p, iF, dict),
334  surfacesDict_(dict.subDict("geometry")),
335  projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
336  projectDir_(dict.lookup("projectDirection")),
337  wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
338  frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
339 {}
340 
341 
344 (
346  const pointPatch& p,
348  const pointPatchFieldMapper&
349 )
350 :
351  pointPatchVectorField(p, iF),
352  surfacesDict_(ppf.surfacesDict_),
353  projectMode_(ppf.projectMode_),
354  projectDir_(ppf.projectDir_),
355  wedgePlane_(ppf.wedgePlane_),
356  frozenPointsZone_(ppf.frozenPointsZone_)
357 {}
358 
359 
362 (
364 )
365 :
367  surfacesDict_(ppf.surfacesDict_),
368  projectMode_(ppf.projectMode_),
369  projectDir_(ppf.projectDir_),
370  wedgePlane_(ppf.wedgePlane_),
371  frozenPointsZone_(ppf.frozenPointsZone_)
372 {}
373 
374 
377 (
380 )
381 :
382  pointPatchVectorField(ppf, iF),
383  surfacesDict_(ppf.surfacesDict_),
384  projectMode_(ppf.projectMode_),
385  projectDir_(ppf.projectDir_),
386  wedgePlane_(ppf.wedgePlane_),
387  frozenPointsZone_(ppf.frozenPointsZone_)
388 {}
389 
390 
391 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
392 
393 const searchableSurfaces&
395 {
396  if (surfacesPtr_.empty())
397  {
398  surfacesPtr_.reset
399  (
401  (
402  IOobject
403  (
404  "abc", // dummy name
405  db().time().constant(), // directory
406  "triSurface", // instance
407  db().time(), // registry
410  ),
411  surfacesDict_
412  )
413  );
414  }
415  return surfacesPtr_();
416 }
417 
418 
420 (
421  const Pstream::commsTypes commsType
422 )
423 {
424  vectorField displacement(this->patchInternalField());
425 
426  // Calculate displacement to project points onto surface
427  calcProjection(displacement);
428 
429  // Get internal field to insert values into
430  Field<vector>& iF = const_cast<Field<vector>&>(this->internalField());
431 
432  //setInInternalField(iF, motionU);
433  setInInternalField(iF, displacement);
434 
436 }
437 
438 
440 {
442  os.writeKeyword("geometry") << surfacesDict_
443  << token::END_STATEMENT << nl;
444  os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
445  << token::END_STATEMENT << nl;
446  os.writeKeyword("projectDirection") << projectDir_
447  << token::END_STATEMENT << nl;
448  os.writeKeyword("wedgePlane") << wedgePlane_
449  << token::END_STATEMENT << nl;
450  if (frozenPointsZone_ != word::null)
451  {
452  os.writeKeyword("frozenPointsZone") << frozenPointsZone_
453  << token::END_STATEMENT << nl;
454  }
455 }
456 
457 
458 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
459 
461 (
464 );
465 
466 
467 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
468 
469 } // End namespace Foam
470 
471 // ************************ vim: set sw=4 sts=4 et: ************************ //