FreeFOAM The Cross-Platform CFD Toolkit
sampledIsoSurface.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 "sampledIsoSurface.H"
27 #include <OpenFOAM/dictionary.H>
28 #include <finiteVolume/volFields.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(sampledIsoSurface, 0);
38  (
39  sampledSurface,
40  sampledIsoSurface,
41  word,
42  isoSurface
43  );
44 }
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 void Foam::sampledIsoSurface::getIsoFields() const
49 {
50  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
51 
52  // Get volField
53  // ~~~~~~~~~~~~
54 
55  if (fvm.foundObject<volScalarField>(isoField_))
56  {
57  if (debug)
58  {
59  Info<< "sampledIsoSurface::getIsoField() : lookup volField "
60  << isoField_ << endl;
61  }
62  storedVolFieldPtr_.clear();
63  volFieldPtr_ = &fvm.lookupObject<volScalarField>(isoField_);
64  }
65  else
66  {
67  // Bit of a hack. Read field and store.
68 
69  if (debug)
70  {
71  Info<< "sampledIsoSurface::getIsoField() : checking "
72  << isoField_ << " for same time " << fvm.time().timeName()
73  << endl;
74  }
75 
76  if
77  (
78  storedVolFieldPtr_.empty()
79  || (fvm.time().timeName() != storedVolFieldPtr_().instance())
80  )
81  {
82  if (debug)
83  {
84  Info<< "sampledIsoSurface::getIsoField() : reading volField "
85  << isoField_ << " from time " << fvm.time().timeName()
86  << endl;
87  }
88 
89  storedVolFieldPtr_.reset
90  (
91  new volScalarField
92  (
93  IOobject
94  (
95  isoField_,
96  fvm.time().timeName(),
97  fvm,
100  false
101  ),
102  fvm
103  )
104  );
105  volFieldPtr_ = storedVolFieldPtr_.operator->();
106  }
107  }
108 
109 
110 
111  // Get pointField
112  // ~~~~~~~~~~~~~~
113 
114  if (!subMeshPtr_.valid())
115  {
116  word pointFldName = "volPointInterpolate(" + isoField_ + ')';
117 
118  if (fvm.foundObject<pointScalarField>(pointFldName))
119  {
120  if (debug)
121  {
122  Info<< "sampledIsoSurface::getIsoField() : lookup pointField "
123  << pointFldName << endl;
124  }
125  pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
126  }
127  else
128  {
129  // Not in registry. Interpolate.
130 
131  if (debug)
132  {
133  Info<< "sampledIsoSurface::getIsoField() : checking pointField "
134  << pointFldName << " for same time "
135  << fvm.time().timeName() << endl;
136  }
137 
138  if
139  (
140  storedPointFieldPtr_.empty()
141  || (fvm.time().timeName() != storedPointFieldPtr_().instance())
142  )
143  {
144  if (debug)
145  {
146  Info<< "sampledIsoSurface::getIsoField() :"
147  << " interpolating volField " << volFieldPtr_->name()
148  << " to get pointField " << pointFldName << endl;
149  }
150 
151  storedPointFieldPtr_.reset
152  (
154  .interpolate(*volFieldPtr_).ptr()
155  );
156  storedPointFieldPtr_->checkOut();
157  pointFieldPtr_ = storedPointFieldPtr_.operator->();
158  }
159  }
160 
161 
162  // If averaging redo the volField. Can only be done now since needs the
163  // point field.
164  if (average_)
165  {
166  storedVolFieldPtr_.reset(average(fvm, *pointFieldPtr_).ptr());
167  volFieldPtr_ = storedVolFieldPtr_.operator->();
168  }
169 
170 
171  if (debug)
172  {
173  Info<< "sampledIsoSurface::getIsoField() : volField "
174  << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
175  << " max:" << max(*volFieldPtr_).value() << endl;
176  Info<< "sampledIsoSurface::getIsoField() : pointField "
177  << pointFieldPtr_->name()
178  << " min:" << gMin(pointFieldPtr_->internalField())
179  << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
180  }
181  }
182  else
183  {
184  // Get subMesh variants
185  const fvMesh& subFvm = subMeshPtr_().subMesh();
186 
187  // Either lookup on the submesh or subset the whole-mesh volField
188 
189  if (subFvm.foundObject<volScalarField>(isoField_))
190  {
191  if (debug)
192  {
193  Info<< "sampledIsoSurface::getIsoField() :"
194  << " submesh lookup volField "
195  << isoField_ << endl;
196  }
197  storedVolSubFieldPtr_.clear();
198  volSubFieldPtr_ = &subFvm.lookupObject<volScalarField>(isoField_);
199  }
200  else
201  {
202  if (debug)
203  {
204  Info<< "sampledIsoSurface::getIsoField() : subsetting volField "
205  << isoField_ << endl;
206  }
207  storedVolSubFieldPtr_.reset
208  (
209  subMeshPtr_().interpolate
210  (
211  *volFieldPtr_
212  ).ptr()
213  );
214  storedVolSubFieldPtr_->checkOut();
215  volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
216  }
217 
218 
219  // Pointfield on submesh
220 
221  word pointFldName =
222  "volPointInterpolate("
223  + volSubFieldPtr_->name()
224  + ')';
225 
226  if (subFvm.foundObject<pointScalarField>(pointFldName))
227  {
228  if (debug)
229  {
230  Info<< "sampledIsoSurface::getIsoField() :"
231  << " submesh lookup pointField " << pointFldName << endl;
232  }
233  storedPointSubFieldPtr_.clear();
234  pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
235  (
236  pointFldName
237  );
238  }
239  else
240  {
241  if (debug)
242  {
243  Info<< "sampledIsoSurface::getIsoField() :"
244  << " interpolating submesh volField "
245  << volSubFieldPtr_->name()
246  << " to get submesh pointField " << pointFldName << endl;
247  }
248  storedPointSubFieldPtr_.reset
249  (
251  (
252  subFvm
253  ).interpolate(*volSubFieldPtr_).ptr()
254  );
255  storedPointSubFieldPtr_->checkOut();
256  pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
257  }
258 
259 
260 
261  // If averaging redo the volField. Can only be done now since needs the
262  // point field.
263  if (average_)
264  {
265  storedVolSubFieldPtr_.reset
266  (
267  average(subFvm, *pointSubFieldPtr_).ptr()
268  );
269  volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
270  }
271 
272 
273  if (debug)
274  {
275  Info<< "sampledIsoSurface::getIsoField() : volSubField "
276  << volSubFieldPtr_->name()
277  << " min:" << min(*volSubFieldPtr_).value()
278  << " max:" << max(*volSubFieldPtr_).value() << endl;
279  Info<< "sampledIsoSurface::getIsoField() : pointSubField "
280  << pointSubFieldPtr_->name()
281  << " min:" << gMin(pointSubFieldPtr_->internalField())
282  << " max:" << gMax(pointSubFieldPtr_->internalField()) << endl;
283  }
284  }
285 }
286 
287 
289 (
290  const fvMesh& mesh,
291  const pointScalarField& pfld
292 ) const
293 {
294  tmp<volScalarField> tcellAvg
295  (
296  new volScalarField
297  (
298  IOobject
299  (
300  "cellAvg",
301  mesh.time().timeName(),
302  mesh,
305  false
306  ),
307  mesh,
308  dimensionedScalar("zero", dimless, scalar(0.0))
309  )
310  );
311  volScalarField& cellAvg = tcellAvg();
312 
313  labelField nPointCells(mesh.nCells(), 0);
314  {
315  for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
316  {
317  const labelList& pCells = mesh.pointCells(pointI);
318 
319  forAll(pCells, i)
320  {
321  label cellI = pCells[i];
322 
323  cellAvg[cellI] += pfld[pointI];
324  nPointCells[cellI]++;
325  }
326  }
327  }
328  forAll(cellAvg, cellI)
329  {
330  cellAvg[cellI] /= nPointCells[cellI];
331  }
332  // Give value to calculatedFvPatchFields
333  cellAvg.correctBoundaryConditions();
334 
335  return tcellAvg;
336 }
337 
338 
340 (
341  const pointMesh& pMesh,
342  const volScalarField& fld
343 ) const
344 {
345  tmp<pointScalarField> tpointAvg
346  (
347  new pointScalarField
348  (
349  IOobject
350  (
351  "pointAvg",
352  fld.time().timeName(),
353  fld.db(),
356  false
357  ),
358  pMesh,
359  dimensionedScalar("zero", dimless, scalar(0.0))
360  )
361  );
362  pointScalarField& pointAvg = tpointAvg();
363 
364  for (label pointI = 0; pointI < fld.mesh().nPoints(); pointI++)
365  {
366  const labelList& pCells = fld.mesh().pointCells(pointI);
367 
368  forAll(pCells, i)
369  {
370  pointAvg[pointI] += fld[pCells[i]];
371  }
372  pointAvg[pointI] /= pCells.size();
373  }
374  // Give value to calculatedFvPatchFields
375  pointAvg.correctBoundaryConditions();
376 
377  return tpointAvg;
378 }
379 
380 
381 bool Foam::sampledIsoSurface::updateGeometry() const
382 {
383  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
384 
385  // no update needed
386  if (fvm.time().timeIndex() == prevTimeIndex_)
387  {
388  return false;
389  }
390 
391  // Get any subMesh
392  if (zoneID_.index() != -1 && !subMeshPtr_.valid())
393  {
394  const polyBoundaryMesh& patches = mesh().boundaryMesh();
395 
396  // Patch to put exposed internal faces into
397  label exposedPatchI = patches.findPatchID(exposedPatchName_);
398 
399  if (debug)
400  {
401  Info<< "Allocating subset of size "
402  << mesh().cellZones()[zoneID_.index()].size()
403  << " with exposed faces into patch "
404  << patches[exposedPatchI].name() << endl;
405  }
406 
407  subMeshPtr_.reset
408  (
409  new fvMeshSubset(fvm)
410  );
411  subMeshPtr_().setLargeCellSubset
412  (
413  labelHashSet(mesh().cellZones()[zoneID_.index()]),
414  exposedPatchI
415  );
416  }
417 
418 
419  prevTimeIndex_ = fvm.time().timeIndex();
420  getIsoFields();
421 
422  // Clear any stored topo
423  surfPtr_.clear();
424  facesPtr_.clear();
425 
426  // Clear derived data
427  clearGeom();
428 
429  if (subMeshPtr_.valid())
430  {
431  surfPtr_.reset
432  (
433  new isoSurface
434  (
435  *volSubFieldPtr_,
436  *pointSubFieldPtr_,
437  isoVal_,
438  regularise_,
439  mergeTol_
440  )
441  );
442  }
443  else
444  {
445  surfPtr_.reset
446  (
447  new isoSurface
448  (
449  *volFieldPtr_,
450  *pointFieldPtr_,
451  isoVal_,
452  regularise_,
453  mergeTol_
454  )
455  );
456  }
457 
458 
459  if (debug)
460  {
461  Pout<< "sampledIsoSurface::updateGeometry() : constructed iso:"
462  << nl
463  << " regularise : " << regularise_ << nl
464  << " average : " << average_ << nl
465  << " isoField : " << isoField_ << nl
466  << " isoValue : " << isoVal_ << nl;
467  if (subMeshPtr_.valid())
468  {
469  Pout<< " zone size : " << subMeshPtr_().subMesh().nCells()
470  << nl;
471  }
472  Pout<< " points : " << points().size() << nl
473  << " tris : " << surface().size() << nl
474  << " cut cells : " << surface().meshCells().size()
475  << endl;
476  }
477 
478  return true;
479 }
480 
481 
482 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
483 
485 (
486  const word& name,
487  const polyMesh& mesh,
488  const dictionary& dict
489 )
490 :
491  sampledSurface(name, mesh, dict),
492  isoField_(dict.lookup("isoField")),
493  isoVal_(readScalar(dict.lookup("isoValue"))),
494  mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
495  regularise_(dict.lookupOrDefault("regularise", true)),
496  average_(dict.lookupOrDefault("average", false)),
497  zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
498  exposedPatchName_(word::null),
499  surfPtr_(NULL),
500  facesPtr_(NULL),
501  prevTimeIndex_(-1),
502  storedVolFieldPtr_(NULL),
503  volFieldPtr_(NULL),
504  storedPointFieldPtr_(NULL),
505  pointFieldPtr_(NULL)
506 {
508  {
510  (
511  "sampledIsoSurface::sampledIsoSurface"
512  "(const word&, const polyMesh&, const dictionary&)",
513  dict
514  ) << "Non-interpolated iso surface not supported since triangles"
515  << " span across cells." << exit(FatalIOError);
516  }
517 
518  if (zoneID_.index() != -1)
519  {
520  dict.lookup("exposedPatchName") >> exposedPatchName_;
521 
522  if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
523  {
525  (
526  "sampledIsoSurface::sampledIsoSurface"
527  "(const word&, const polyMesh&, const dictionary&)",
528  dict
529  ) << "Cannot find patch " << exposedPatchName_
530  << " in which to put exposed faces." << endl
531  << "Valid patches are " << mesh.boundaryMesh().names()
532  << exit(FatalIOError);
533  }
534 
535  if (debug && zoneID_.index() != -1)
536  {
537  Info<< "Restricting to cellZone " << zoneID_.name()
538  << " with exposed internal faces into patch "
539  << exposedPatchName_ << endl;
540  }
541  }
542 }
543 
544 
545 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
546 
548 {}
549 
550 
551 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
552 
554 {
555  const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
556 
557  return fvm.time().timeIndex() != prevTimeIndex_;
558 }
559 
560 
562 {
563  surfPtr_.clear();
564  facesPtr_.clear();
565  subMeshPtr_.clear();
566 
567  // Clear derived data
568  clearGeom();
569 
570  // already marked as expired
571  if (prevTimeIndex_ == -1)
572  {
573  return false;
574  }
575 
576  // force update
577  prevTimeIndex_ = -1;
578  return true;
579 }
580 
581 
583 {
584  return updateGeometry();
585 }
586 
587 
589 (
590  const volScalarField& vField
591 ) const
592 {
593  return sampleField(vField);
594 }
595 
596 
598 (
599  const volVectorField& vField
600 ) const
601 {
602  return sampleField(vField);
603 }
604 
605 
607 (
608  const volSphericalTensorField& vField
609 ) const
610 {
611  return sampleField(vField);
612 }
613 
614 
616 (
617  const volSymmTensorField& vField
618 ) const
619 {
620  return sampleField(vField);
621 }
622 
623 
625 (
626  const volTensorField& vField
627 ) const
628 {
629  return sampleField(vField);
630 }
631 
632 
634 (
635  const interpolation<scalar>& interpolator
636 ) const
637 {
638  return interpolateField(interpolator);
639 }
640 
641 
643 (
644  const interpolation<vector>& interpolator
645 ) const
646 {
647  return interpolateField(interpolator);
648 }
649 
651 (
652  const interpolation<sphericalTensor>& interpolator
653 ) const
654 {
655  return interpolateField(interpolator);
656 }
657 
658 
660 (
661  const interpolation<symmTensor>& interpolator
662 ) const
663 {
664  return interpolateField(interpolator);
665 }
666 
667 
669 (
670  const interpolation<tensor>& interpolator
671 ) const
672 {
673  return interpolateField(interpolator);
674 }
675 
676 
678 {
679  os << "sampledIsoSurface: " << name() << " :"
680  << " field :" << isoField_
681  << " value :" << isoVal_;
682  //<< " faces:" << faces().size() // note: possibly no geom yet
683  //<< " points:" << points().size();
684 }
685 
686 
687 // ************************ vim: set sw=4 sts=4 et: ************************ //