FreeFOAM The Cross-Platform CFD Toolkit
displacementInterpolationFvMotionSolver.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/SortableList.H>
29 #include <OpenFOAM/IOList.H>
30 #include <OpenFOAM/Tuple2.H>
31 #include <OpenFOAM/mapPolyMesh.H>
32 #include <OpenFOAM/interpolateXY.H>
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(displacementInterpolationFvMotionSolver, 0);
39 
41  (
42  fvMotionSolver,
43  displacementInterpolationFvMotionSolver,
44  dictionary
45  );
46 
47  template<>
48  const word IOList<Tuple2<scalar, vector> >::typeName("scalarVectorTable");
49 }
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
55 
56 Foam::displacementInterpolationFvMotionSolver::
57 displacementInterpolationFvMotionSolver
58 (
59  const polyMesh& mesh,
60  Istream& is
61 )
62 :
64  dynamicMeshCoeffs_
65  (
67  (
68  IOobject
69  (
70  "dynamicMeshDict",
71  mesh.time().constant(),
72  mesh,
73  IOobject::MUST_READ,
74  IOobject::NO_WRITE,
75  false
76  )
77  ).subDict(typeName + "Coeffs")
78  )
79 {
80  // Get zones and their interpolation tables for displacement
81  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
82 
83  List<Pair<word> > faceZoneToTable
84  (
85  dynamicMeshCoeffs_.lookup("interpolationTables")
86  );
87 
88  const faceZoneMesh& fZones = mesh.faceZones();
89 
90  times_.setSize(fZones.size());
91  displacements_.setSize(fZones.size());
92 
93  forAll(faceZoneToTable, i)
94  {
95  const word& zoneName = faceZoneToTable[i][0];
96  label zoneI = fZones.findZoneID(zoneName);
97 
98  if (zoneI == -1)
99  {
101  (
102  "displacementInterpolationFvMotionSolver::"
103  "displacementInterpolationFvMotionSolver(const polyMesh&,"
104  "Istream&)"
105  ) << "Cannot find zone " << zoneName << endl
106  << "Valid zones are " << mesh.faceZones().names()
107  << exit(FatalError);
108  }
109 
110  const word& tableName = faceZoneToTable[i][1];
111 
113  (
114  IOobject
115  (
116  tableName,
117  mesh.time().constant(),
118  "tables",
119  mesh,
120  IOobject::MUST_READ,
121  IOobject::NO_WRITE,
122  false
123  )
124  );
125 
126  // Copy table
127  times_[zoneI].setSize(table.size());
128  displacements_[zoneI].setSize(table.size());
129 
130  forAll(table, j)
131  {
132  times_[zoneI][j] = table[j].first();
133  displacements_[zoneI][j] = table[j].second();
134  }
135  }
136 
137 
138 
139  // Sort points into bins according to position relative to faceZones
140  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141  // Done in all three directions.
142 
143  for (direction dir = 0; dir < vector::nComponents; dir++)
144  {
145  // min and max coordinates of all faceZones
146  SortableList<scalar> zoneCoordinates(2*faceZoneToTable.size());
147 
148  forAll(faceZoneToTable, i)
149  {
150  const word& zoneName = faceZoneToTable[i][0];
151  label zoneI = fZones.findZoneID(zoneName);
152  const faceZone& fz = fZones[zoneI];
153 
154  scalar minCoord = VGREAT;
155  scalar maxCoord = -VGREAT;
156 
157  forAll(fz().meshPoints(), localI)
158  {
159  label pointI = fz().meshPoints()[localI];
160  const scalar coord = points0()[pointI][dir];
161  minCoord = min(minCoord, coord);
162  maxCoord = max(maxCoord, coord);
163  }
164 
165  zoneCoordinates[2*i] = returnReduce(minCoord, minOp<scalar>());
166  zoneCoordinates[2*i+1] = returnReduce(maxCoord, maxOp<scalar>());
167 
168  if (debug)
169  {
170  Pout<< "direction " << dir << " : "
171  << "zone " << zoneName
172  << " ranges from coordinate " << zoneCoordinates[2*i]
173  << " to " << zoneCoordinates[2*i+1]
174  << endl;
175  }
176  }
177  zoneCoordinates.sort();
178 
179  // Slightly tweak min and max face zone so points sort within
180  zoneCoordinates[0] -= SMALL;
181  zoneCoordinates[zoneCoordinates.size()-1] += SMALL;
182 
183  // Check if we have static min and max mesh bounds
184  const scalarField meshCoords = points0().component(dir);
185 
186  scalar minCoord = gMin(meshCoords);
187  scalar maxCoord = gMax(meshCoords);
188 
189  if (debug)
190  {
191  Pout<< "direction " << dir << " : "
192  << "mesh ranges from coordinate " << minCoord << " to "
193  << maxCoord << endl;
194  }
195 
196  // Make copy of zoneCoordinates; include min and max of mesh
197  // if necessary. Mark min and max with zoneI=-1.
198  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
199 
200  labelList& rangeZone = rangeToZone_[dir];
201  labelListList& rangePoints = rangeToPoints_[dir];
202  List<scalarField>& rangeWeights = rangeToWeights_[dir];
203 
204  scalarField rangeToCoord(zoneCoordinates.size());
205  rangeZone.setSize(zoneCoordinates.size());
206  label rangeI = 0;
207 
208  if (minCoord < zoneCoordinates[0])
209  {
210  label sz = rangeZone.size();
211  rangeToCoord.setSize(sz+1);
212  rangeZone.setSize(sz+1);
213  rangeToCoord[rangeI] = minCoord-SMALL;
214  rangeZone[rangeI] = -1;
215 
216  if (debug)
217  {
218  Pout<< "direction " << dir << " : "
219  << "range " << rangeI << " at coordinate "
220  << rangeToCoord[rangeI] << " from min of mesh "
221  << rangeZone[rangeI] << endl;
222  }
223  rangeI = 1;
224  }
225  forAll(zoneCoordinates, i)
226  {
227  rangeToCoord[rangeI] = zoneCoordinates[i];
228  rangeZone[rangeI] = zoneCoordinates.indices()[i]/2;
229 
230  if (debug)
231  {
232  Pout<< "direction " << dir << " : "
233  << "range " << rangeI << " at coordinate "
234  << rangeToCoord[rangeI]
235  << " from zone " << rangeZone[rangeI] << endl;
236  }
237  rangeI++;
238  }
239  if (maxCoord > zoneCoordinates[zoneCoordinates.size()-1])
240  {
241  label sz = rangeToCoord.size();
242  rangeToCoord.setSize(sz+1);
243  rangeZone.setSize(sz+1);
244  rangeToCoord[sz] = maxCoord+SMALL;
245  rangeZone[sz] = -1;
246 
247  if (debug)
248  {
249  Pout<< "direction " << dir << " : "
250  << "range " << rangeI << " at coordinate "
251  << rangeToCoord[sz] << " from max of mesh "
252  << rangeZone[sz] << endl;
253  }
254  }
255 
256 
257  // Sort the points
258  // ~~~~~~~~~~~~~~~
259 
260  // Count all the points inbetween rangeI and rangeI+1
261  labelList nRangePoints(rangeToCoord.size(), 0);
262 
263  forAll(meshCoords, pointI)
264  {
265  label rangeI = findLower(rangeToCoord, meshCoords[pointI]);
266 
267  if (rangeI == -1 || rangeI == rangeToCoord.size()-1)
268  {
270  (
271  "displacementInterpolationFvMotionSolver::"
272  "displacementInterpolationFvMotionSolver"
273  "(const polyMesh&, Istream&)"
274  ) << "Did not find point " << points0()[pointI]
275  << " coordinate " << meshCoords[pointI]
276  << " in ranges " << rangeToCoord
277  << abort(FatalError);
278  }
279  nRangePoints[rangeI]++;
280  }
281 
282  if (debug)
283  {
284  for (label rangeI = 0; rangeI < rangeToCoord.size()-1; rangeI++)
285  {
286  // Get the two zones bounding the range
287  Pout<< "direction " << dir << " : "
288  << "range from " << rangeToCoord[rangeI]
289  << " to " << rangeToCoord[rangeI+1]
290  << " contains " << nRangePoints[rangeI]
291  << " points." << endl;
292  }
293  }
294 
295  // Sort
296  rangePoints.setSize(nRangePoints.size());
297  rangeWeights.setSize(nRangePoints.size());
298  forAll(rangePoints, rangeI)
299  {
300  rangePoints[rangeI].setSize(nRangePoints[rangeI]);
301  rangeWeights[rangeI].setSize(nRangePoints[rangeI]);
302  }
303  nRangePoints = 0;
304  forAll(meshCoords, pointI)
305  {
306  label rangeI = findLower(rangeToCoord, meshCoords[pointI]);
307  label& nPoints = nRangePoints[rangeI];
308  rangePoints[rangeI][nPoints] = pointI;
309  rangeWeights[rangeI][nPoints] =
310  (meshCoords[pointI]-rangeToCoord[rangeI])
311  / (rangeToCoord[rangeI+1]-rangeToCoord[rangeI]);
312  nPoints++;
313  }
314  }
315 }
316 
317 
318 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
319 
322 {}
323 
324 
325 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
326 
329 {
330  if (mesh().nPoints() != points0().size())
331  {
333  (
334  "displacementInterpolationFvMotionSolver::curPoints() const"
335  ) << "The number of points in the mesh seems to have changed." << endl
336  << "In constant/polyMesh there are " << points0().size()
337  << " points; in the current mesh there are " << mesh().nPoints()
338  << " points." << exit(FatalError);
339  }
340 
341  tmp<pointField> tcurPoints(new pointField(points0()));
342  pointField& curPoints = tcurPoints();
343 
344  // Interpolate the diplacement of the face zones.
345  vectorField zoneDisp(displacements_.size(), vector::zero);
346  forAll(zoneDisp, zoneI)
347  {
348  if (times_[zoneI].size())
349  {
350  zoneDisp[zoneI] = interpolateXY
351  (
352  mesh().time().value(),
353  times_[zoneI],
354  displacements_[zoneI]
355  );
356  }
357  }
358  if (debug)
359  {
360  Pout<< "Zone displacements:" << zoneDisp << endl;
361  }
362 
363 
364  // Interpolate the point location
365  for (direction dir = 0; dir < vector::nComponents; dir++)
366  {
367  const labelList& rangeZone = rangeToZone_[dir];
368  const labelListList& rangePoints = rangeToPoints_[dir];
369  const List<scalarField>& rangeWeights = rangeToWeights_[dir];
370 
371  for (label rangeI = 0; rangeI < rangeZone.size()-1; rangeI++)
372  {
373  const labelList& rPoints = rangePoints[rangeI];
374  const scalarField& rWeights = rangeWeights[rangeI];
375 
376  // Get the two zones bounding the range
377  label minZoneI = rangeZone[rangeI];
378  //vector minDisp =
379  // (minZoneI == -1 ? vector::zero : zoneDisp[minZoneI]);
380  scalar minDisp = (minZoneI == -1 ? 0.0 : zoneDisp[minZoneI][dir]);
381  label maxZoneI = rangeZone[rangeI+1];
382  //vector maxDisp =
383  // (maxZoneI == -1 ? vector::zero : zoneDisp[maxZoneI]);
384  scalar maxDisp = (maxZoneI == -1 ? 0.0 : zoneDisp[maxZoneI][dir]);
385 
386  forAll(rPoints, i)
387  {
388  label pointI = rPoints[i];
389  scalar w = rWeights[i];
390  //curPoints[pointI] += (1.0-w)*minDisp+w*maxDisp;
391  curPoints[pointI][dir] += (1.0-w)*minDisp+w*maxDisp;
392  }
393  }
394  }
395  return tcurPoints;
396 }
397 
398 
399 // ************************ vim: set sw=4 sts=4 et: ************************ //