FreeFOAM The Cross-Platform CFD Toolkit
displacementComponentLaplacianFvMotionSolver.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 
30 #include <OpenFOAM/mapPolyMesh.H>
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(displacementComponentLaplacianFvMotionSolver, 0);
38 
40  (
41  fvMotionSolver,
42  displacementComponentLaplacianFvMotionSolver,
43  dictionary
44  );
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 Foam::direction Foam::displacementComponentLaplacianFvMotionSolver::cmpt
51 (
52  const word& cmptName
53 ) const
54 {
55  if (cmptName == "x")
56  {
57  return vector::X;
58  }
59  else if (cmptName == "y")
60  {
61  return vector::Y;
62  }
63  else if (cmptName == "z")
64  {
65  return vector::Z;
66  }
67  else
68  {
70  (
71  "displacementComponentLaplacianFvMotionSolver::"
72  "displacementComponentLaplacianFvMotionSolver"
73  "(const polyMesh& mesh, Istream& msData)"
74  ) << "Given component name " << cmptName << " should be x, y or z"
75  << exit(FatalError);
76 
77  return 0;
78  }
79 }
80 
81 
82 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
83 
84 Foam::displacementComponentLaplacianFvMotionSolver::
85 displacementComponentLaplacianFvMotionSolver
86 (
87  const polyMesh& mesh,
88  Istream& msData
89 )
90 :
91  fvMotionSolver(mesh),
92  cmptName_(msData),
93  cmpt_(cmpt(cmptName_)),
94  points0_
95  (
97  (
98  IOobject
99  (
100  "points",
101  time().constant(),
103  mesh,
106  false
107  )
108  ).component(cmpt_)
109  ),
110  pointDisplacement_
111  (
112  IOobject
113  (
114  "pointDisplacement" + cmptName_,
115  fvMesh_.time().timeName(),
116  fvMesh_,
119  ),
120  pointMesh::New(fvMesh_)
121  ),
122  cellDisplacement_
123  (
124  IOobject
125  (
126  "cellDisplacement" + cmptName_,
127  mesh.time().timeName(),
128  mesh,
131  ),
132  fvMesh_,
134  (
135  "cellDisplacement",
136  pointDisplacement_.dimensions(),
137  0
138  ),
139  cellMotionBoundaryTypes<scalar>(pointDisplacement_.boundaryField())
140  ),
141  pointLocation_(NULL),
142  diffusivityPtr_
143  (
144  motionDiffusivity::New(*this, lookup("diffusivity"))
145  ),
146  frozenPointsZone_
147  (
148  found("frozenPointsZone")
149  ? fvMesh_.pointZones().findZoneID(lookup("frozenPointsZone"))
150  : -1
151  )
152 {
153  IOobject io
154  (
155  "pointLocation",
156  fvMesh_.time().timeName(),
157  fvMesh_,
160  );
161 
162  if (debug)
163  {
164  Info<< "displacementComponentLaplacianFvMotionSolver:" << nl
165  << " diffusivity : " << diffusivityPtr_().type() << nl
166  << " frozenPoints zone : " << frozenPointsZone_ << endl;
167  }
168 
169  if (io.headerOk())
170  {
171  pointLocation_.reset
172  (
173  new pointVectorField
174  (
175  io,
176  pointMesh::New(fvMesh_)
177  )
178  );
179 
180  if (debug)
181  {
182  Info<< "displacementComponentLaplacianFvMotionSolver :"
183  << " Read pointVectorField "
184  << io.name() << " to be used for boundary conditions on points."
185  << nl
186  << "Boundary conditions:"
187  << pointLocation_().boundaryField().types() << endl;
188  }
189  }
190 }
191 
192 
193 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
194 
197 {}
198 
199 
200 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
201 
204 {
206  (
207  cellDisplacement_,
208  pointDisplacement_
209  );
210 
211  if (pointLocation_.valid())
212  {
213  if (debug)
214  {
215  Info<< "displacementComponentLaplacianFvMotionSolver : applying "
216  << " boundary conditions on " << pointLocation_().name()
217  << " to new point location."
218  << endl;
219  }
220 
221  // Apply pointLocation_ b.c. to mesh points.
222 
223  pointLocation_().internalField() = fvMesh_.points();
224 
225  pointLocation_().internalField().replace
226  (
227  cmpt_,
228  points0_ + pointDisplacement_.internalField()
229  );
230 
231  pointLocation_().correctBoundaryConditions();
232 
233  // Implement frozen points
234  if (frozenPointsZone_ != -1)
235  {
236  const pointZone& pz = fvMesh_.pointZones()[frozenPointsZone_];
237 
238  forAll(pz, i)
239  {
240  label pointI = pz[i];
241 
242  pointLocation_()[pointI][cmpt_] = points0_[pointI];
243  }
244  }
245 
246  twoDCorrectPoints(pointLocation_().internalField());
247 
248  return tmp<pointField>(pointLocation_().internalField());
249  }
250  else
251  {
252  tmp<pointField> tcurPoints(new pointField(fvMesh_.points()));
253 
254  tcurPoints().replace
255  (
256  cmpt_,
257  points0_ + pointDisplacement_.internalField()
258  );
259 
260  // Implement frozen points
261  if (frozenPointsZone_ != -1)
262  {
263  const pointZone& pz = fvMesh_.pointZones()[frozenPointsZone_];
264 
265  forAll(pz, i)
266  {
267  label pointI = pz[i];
268 
269  tcurPoints()[pointI][cmpt_] = points0_[pointI];
270  }
271  }
272 
273  twoDCorrectPoints(tcurPoints());
274 
275  return tcurPoints;
276  }
277 }
278 
279 
281 {
282  // The points have moved so before interpolation update
283  // the fvMotionSolver accordingly
284  movePoints(fvMesh_.points());
285 
286  diffusivityPtr_->correct();
287  pointDisplacement_.boundaryField().updateCoeffs();
288 
290  (
292  (
293  diffusivityPtr_->operator()(),
294  cellDisplacement_,
295  "laplacian(diffusivity,cellDisplacement)"
296  )
297  );
298 }
299 
300 
302 (
303  const mapPolyMesh& mpm
304 )
305 {
307 
308  // Map points0_. Bit special since we somehow have to come up with
309  // a sensible points0 position for introduced points.
310  // Find out scaling between points0 and current points
311 
312  // Get the new points either from the map or the mesh
313  const scalarField points
314  (
315  mpm.hasMotionPoints()
316  ? mpm.preMotionPoints().component(cmpt_)
317  : fvMesh_.points().component(cmpt_)
318  );
319 
320  // Get extents of points0 and points and determine scale
321  const scalar scale =
322  (gMax(points0_)-gMin(points0_))
323  /(gMax(points)-gMin(points));
324 
325  scalarField newPoints0(mpm.pointMap().size());
326 
327  forAll(newPoints0, pointI)
328  {
329  label oldPointI = mpm.pointMap()[pointI];
330 
331  if (oldPointI >= 0)
332  {
333  label masterPointI = mpm.reversePointMap()[oldPointI];
334 
335  if (masterPointI == pointI)
336  {
337  newPoints0[pointI] = points0_[oldPointI];
338  }
339  else
340  {
341  // New point. Assume motion is scaling.
342  newPoints0[pointI] =
343  points0_[oldPointI]
344  + scale*(points[pointI]-points[masterPointI]);
345  }
346  }
347  else
348  {
350  (
351  "displacementLaplacianFvMotionSolver::updateMesh"
352  "(const mapPolyMesh& mpm)"
353  ) << "Cannot work out coordinates of introduced vertices."
354  << " New vertex " << pointI << " at coordinate "
355  << points[pointI] << exit(FatalError);
356  }
357  }
358  points0_.transfer(newPoints0);
359 
360  // Update diffusivity. Note two stage to make sure old one is de-registered
361  // before creating/registering new one.
362  diffusivityPtr_.reset(NULL);
363  diffusivityPtr_ = motionDiffusivity::New(*this, lookup("diffusivity"));
364 }
365 
366 
367 // ************************ vim: set sw=4 sts=4 et: ************************ //