FreeFOAM The Cross-Platform CFD Toolkit
sampledCuttingPlane.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 
26 #include "sampledCuttingPlane.H"
27 #include <OpenFOAM/dictionary.H>
28 #include <finiteVolume/volFields.H>
31 #include <finiteVolume/fvMesh.H>
32 #include <sampling/isoSurface.H>
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(sampledCuttingPlane, 0);
40  (
41  sampledSurface,
42  sampledCuttingPlane,
43  word,
44  cuttingPlane
45  );
46 }
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void Foam::sampledCuttingPlane::createGeometry()
51 {
52  if (debug)
53  {
54  Pout<< "sampledCuttingPlane::createGeometry :updating geometry."
55  << endl;
56  }
57 
58  // Clear any stored topologies
59  facesPtr_.clear();
60  isoSurfPtr_.ptr();
61  pointDistance_.clear();
62  cellDistancePtr_.clear();
63 
64  // Clear derived data
65  clearGeom();
66 
67  // Get any subMesh
68  if (zoneID_.index() != -1 && !subMeshPtr_.valid())
69  {
70  const polyBoundaryMesh& patches = mesh().boundaryMesh();
71 
72  // Patch to put exposed internal faces into
73  label exposedPatchI = patches.findPatchID(exposedPatchName_);
74 
75  if (debug)
76  {
77  Info<< "Allocating subset of size "
78  << mesh().cellZones()[zoneID_.index()].size()
79  << " with exposed faces into patch "
80  << patches[exposedPatchI].name() << endl;
81  }
82 
83  subMeshPtr_.reset
84  (
85  new fvMeshSubset(static_cast<const fvMesh&>(mesh()))
86  );
87  subMeshPtr_().setLargeCellSubset
88  (
89  labelHashSet(mesh().cellZones()[zoneID_.index()]),
90  exposedPatchI
91  );
92  }
93 
94 
95  // Select either the submesh or the underlying mesh
96  const fvMesh& fvm =
97  (
98  subMeshPtr_.valid()
99  ? subMeshPtr_().subMesh()
100  : static_cast<const fvMesh&>(mesh())
101  );
102 
103 
104  // Distance to cell centres
105  // ~~~~~~~~~~~~~~~~~~~~~~~~
106 
107  cellDistancePtr_.reset
108  (
109  new volScalarField
110  (
111  IOobject
112  (
113  "cellDistance",
114  fvm.time().timeName(),
115  fvm.time(),
118  false
119  ),
120  fvm,
121  dimensionedScalar("zero", dimLength, 0)
122  )
123  );
124  volScalarField& cellDistance = cellDistancePtr_();
125 
126  // Internal field
127  {
128  const pointField& cc = fvm.cellCentres();
129  scalarField& fld = cellDistance.internalField();
130 
131  forAll(cc, i)
132  {
133  // Signed distance
134  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
135  }
136  }
137 
138  // Patch fields
139  {
140  forAll(cellDistance.boundaryField(), patchI)
141  {
142  if
143  (
144  isA<emptyFvPatchScalarField>
145  (
146  cellDistance.boundaryField()[patchI]
147  )
148  )
149  {
150  cellDistance.boundaryField().set
151  (
152  patchI,
153  new calculatedFvPatchScalarField
154  (
155  fvm.boundary()[patchI],
156  cellDistance
157  )
158  );
159 
160  const polyPatch& pp = fvm.boundary()[patchI].patch();
161  pointField::subField cc = pp.patchSlice(fvm.faceCentres());
162 
163  fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
164  fld.setSize(pp.size());
165  forAll(fld, i)
166  {
167  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
168  }
169  }
170  else
171  {
172  const pointField& cc = fvm.C().boundaryField()[patchI];
173  fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
174 
175  forAll(fld, i)
176  {
177  fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
178  }
179  }
180  }
181  }
182 
183 
184  // On processor patches the mesh.C() will already be the cell centre
185  // on the opposite side so no need to swap cellDistance.
186 
187 
188  // Distance to points
189  pointDistance_.setSize(fvm.nPoints());
190  {
191  const pointField& pts = fvm.points();
192 
193  forAll(pointDistance_, i)
194  {
195  pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal();
196  }
197  }
198 
199 
200  if (debug)
201  {
202  Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
203  cellDistance.write();
204  pointScalarField pDist
205  (
206  IOobject
207  (
208  "pointDistance",
209  fvm.time().timeName(),
210  fvm.time(),
212  IOobject::NO_WRITE,
213  false
214  ),
215  pointMesh::New(fvm),
216  dimensionedScalar("zero", dimLength, 0)
217  );
218  pDist.internalField() = pointDistance_;
219 
220  Pout<< "Writing point distance:" << pDist.objectPath() << endl;
221  pDist.write();
222  }
223 
224 
225  //- Direct from cell field and point field.
226  isoSurfPtr_.reset
227  (
228  new isoSurface
229  (
230  cellDistance,
231  pointDistance_,
232  0.0,
233  regularise_
234  )
235  );
236 
237  if (debug)
238  {
239  print(Pout);
240  Pout<< endl;
241  }
242 }
243 
244 
245 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
246 
248 (
249  const word& name,
250  const polyMesh& mesh,
251  const dictionary& dict
252 )
253 :
254  sampledSurface(name, mesh, dict),
255  plane_(dict),
256  mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
257  regularise_(dict.lookupOrDefault("regularise", true)),
258  zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
259  exposedPatchName_(word::null),
260  needsUpdate_(true),
261  subMeshPtr_(NULL),
262  cellDistancePtr_(NULL),
263  isoSurfPtr_(NULL),
264  facesPtr_(NULL)
265 {
266  if (zoneID_.index() != -1)
267  {
268  dict.lookup("exposedPatchName") >> exposedPatchName_;
269 
270  if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
271  {
273  (
274  "sampledCuttingPlane::sampledCuttingPlane"
275  "(const word&, const polyMesh&, const dictionary&)"
276  ) << "Cannot find patch " << exposedPatchName_
277  << " in which to put exposed faces." << endl
278  << "Valid patches are " << mesh.boundaryMesh().names()
279  << exit(FatalError);
280  }
281 
282  if (debug && zoneID_.index() != -1)
283  {
284  Info<< "Restricting to cellZone " << zoneID_.name()
285  << " with exposed internal faces into patch "
286  << exposedPatchName_ << endl;
287  }
288  }
289 }
290 
291 
292 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
293 
295 {}
296 
297 
298 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
299 
301 {
302  return needsUpdate_;
303 }
304 
305 
307 {
308  if (debug)
309  {
310  Pout<< "sampledCuttingPlane::expire :"
311  << " have-facesPtr_:" << facesPtr_.valid()
312  << " needsUpdate_:" << needsUpdate_ << endl;
313  }
314 
315  // Clear any stored topologies
316  facesPtr_.clear();
317 
318  // Clear derived data
319  clearGeom();
320 
321  // already marked as expired
322  if (needsUpdate_)
323  {
324  return false;
325  }
326 
327  needsUpdate_ = true;
328  return true;
329 }
330 
331 
333 {
334  if (debug)
335  {
336  Pout<< "sampledCuttingPlane::update :"
337  << " have-facesPtr_:" << facesPtr_.valid()
338  << " needsUpdate_:" << needsUpdate_ << endl;
339  }
340 
341  if (!needsUpdate_)
342  {
343  return false;
344  }
345 
346  createGeometry();
347 
348  needsUpdate_ = false;
349  return true;
350 }
351 
352 
355 (
356  const volScalarField& vField
357 ) const
358 {
359  return sampleField(vField);
360 }
361 
362 
365 (
366  const volVectorField& vField
367 ) const
368 {
369  return sampleField(vField);
370 }
371 
372 
375 (
376  const volSphericalTensorField& vField
377 ) const
378 {
379  return sampleField(vField);
380 }
381 
382 
385 (
386  const volSymmTensorField& vField
387 ) const
388 {
389  return sampleField(vField);
390 }
391 
392 
395 (
396  const volTensorField& vField
397 ) const
398 {
399  return sampleField(vField);
400 }
401 
402 
405 (
406  const interpolation<scalar>& interpolator
407 ) const
408 {
409  return interpolateField(interpolator);
410 }
411 
412 
415 (
416  const interpolation<vector>& interpolator
417 ) const
418 {
419  return interpolateField(interpolator);
420 }
421 
424 (
425  const interpolation<sphericalTensor>& interpolator
426 ) const
427 {
428  return interpolateField(interpolator);
429 }
430 
431 
434 (
435  const interpolation<symmTensor>& interpolator
436 ) const
437 {
438  return interpolateField(interpolator);
439 }
440 
441 
444 (
445  const interpolation<tensor>& interpolator
446 ) const
447 {
448  return interpolateField(interpolator);
449 }
450 
451 
453 {
454  os << "sampledCuttingPlane: " << name() << " :"
455  << " plane:" << plane_
456  << " faces:" << faces().size()
457  << " points:" << points().size();
458 }
459 
460 
461 // ************************ vim: set sw=4 sts=4 et: ************************ //