writeFields.C
Go to the documentation of this file.
1#include "writeFields.H"
2#include "volFields.H"
3#include "surfaceFields.H"
4#include "polyMeshTools.H"
6#include "syncTools.H"
7#include "tetPointRef.H"
8#include "regionSplit.H"
9#include "wallDist.H"
10#include "cellAspectRatio.H"
11
12using namespace Foam;
13
14void maxFaceToCell
15(
16 const scalarField& faceData,
17 volScalarField& cellData
18)
19{
20 const cellList& cells = cellData.mesh().cells();
21
22 scalarField& cellFld = cellData.ref();
23
24 cellFld = -GREAT;
25 forAll(cells, cellI)
26 {
27 const cell& cFaces = cells[cellI];
28 forAll(cFaces, i)
29 {
30 cellFld[cellI] = max(cellFld[cellI], faceData[cFaces[i]]);
31 }
32 }
33
34 forAll(cellData.boundaryField(), patchI)
35 {
36 fvPatchScalarField& fvp = cellData.boundaryFieldRef()[patchI];
37
38 fvp = fvp.patch().patchSlice(faceData);
39 }
41}
42
43
44void minFaceToCell
45(
46 const scalarField& faceData,
47 volScalarField& cellData
48)
49{
50 const cellList& cells = cellData.mesh().cells();
51
52 scalarField& cellFld = cellData.ref();
53
54 cellFld = GREAT;
55 forAll(cells, cellI)
56 {
57 const cell& cFaces = cells[cellI];
58 forAll(cFaces, i)
59 {
60 cellFld[cellI] = min(cellFld[cellI], faceData[cFaces[i]]);
61 }
62 }
63
64 forAll(cellData.boundaryField(), patchI)
65 {
66 fvPatchScalarField& fvp = cellData.boundaryFieldRef()[patchI];
67
68 fvp = fvp.patch().patchSlice(faceData);
69 }
71}
72
73
74void minFaceToCell
75(
76 const surfaceScalarField& faceData,
77 volScalarField& cellData,
79)
80{
81 scalarField& cellFld = cellData.ref();
82
83 cellFld = GREAT;
84
85 const labelUList& own = cellData.mesh().owner();
86 const labelUList& nei = cellData.mesh().neighbour();
87
88 // Internal faces
89 forAll(own, facei)
90 {
91 cellFld[own[facei]] = min(cellFld[own[facei]], faceData[facei]);
92 cellFld[nei[facei]] = min(cellFld[nei[facei]], faceData[facei]);
93 }
94
95 // Patch faces
96 forAll(faceData.boundaryField(), patchi)
97 {
98 const fvsPatchScalarField& fvp = faceData.boundaryField()[patchi];
99 const labelUList& fc = fvp.patch().faceCells();
100
101 forAll(fc, i)
102 {
103 cellFld[fc[i]] = min(cellFld[fc[i]], fvp[i]);
104 }
105 }
106
108
109 forAll(bfld, patchi)
110 {
111 bfld[patchi] = faceData.boundaryField()[patchi];
112 }
114 {
115 cellData.correctBoundaryConditions();
116 }
117}
118
119
120void writeSurfaceField
121(
122 const fvMesh& mesh,
123 const fileName& fName,
124 const scalarField& faceData
125)
126{
127 // Write single surfaceScalarField
128
130 (
132 (
133 fName,
134 mesh.time().timeName(),
135 mesh,
136 IOobject::NO_READ,
137 IOobject::AUTO_WRITE,
138 false
139 ),
140 mesh,
141 dimensionedScalar(dimless, Zero),
142 calculatedFvsPatchScalarField::typeName
143 );
144 fld.primitiveFieldRef() = faceData;
145 //fld.correctBoundaryConditions();
146 Info<< " Writing face data to " << fName << endl;
147 fld.write();
148}
149
150
152(
153 const fvMesh& mesh,
154 const wordHashSet& selectedFields,
155 const bool writeFaceFields
156)
157{
158 if (selectedFields.empty())
159 {
160 return;
161 }
162
163 Info<< "Writing fields with mesh quality parameters" << endl;
164
165 if (selectedFields.found("nonOrthoAngle"))
166 {
167 //- Face based orthogonality
168 const scalarField faceOrthogonality
169 (
170 polyMeshTools::faceOrthogonality
171 (
172 mesh,
173 mesh.faceAreas(),
174 mesh.cellCentres()
175 )
176 );
177
178 //- Face based angle
179 const scalarField nonOrthoAngle
180 (
181 radToDeg
182 (
183 Foam::acos(min(scalar(1), max(scalar(-1), faceOrthogonality)))
184 )
185 );
186
187 //- Cell field - max of either face
188 volScalarField cellNonOrthoAngle
189 (
191 (
192 "nonOrthoAngle",
193 mesh.time().timeName(),
194 mesh,
195 IOobject::NO_READ,
196 IOobject::AUTO_WRITE,
197 false
198 ),
199 mesh,
200 dimensionedScalar(dimless, Zero),
201 calculatedFvPatchScalarField::typeName
202 );
203 //- Take max
204 maxFaceToCell(nonOrthoAngle, cellNonOrthoAngle);
205 Info<< " Writing non-orthogonality (angle) to "
206 << cellNonOrthoAngle.name() << endl;
207 cellNonOrthoAngle.write();
208
209 if (writeFaceFields)
210 {
211 writeSurfaceField
212 (
213 mesh,
214 "face_nonOrthoAngle",
215 SubField<scalar>(nonOrthoAngle, mesh.nInternalFaces())
216 );
217 }
218 }
219
220 if (selectedFields.found("faceWeight"))
221 {
222 volScalarField cellWeights
223 (
225 (
226 "faceWeight",
227 mesh.time().timeName(),
228 mesh,
229 IOobject::NO_READ,
230 IOobject::AUTO_WRITE,
231 false
232 ),
233 mesh,
234 dimensionedScalar(dimless, Zero),
235 wordList // wanted bc types
236 (
237 mesh.boundary().size(),
238 calculatedFvPatchScalarField::typeName
239 ),
240 mesh.weights().boundaryField().types() // current bc types
241 );
242 //- Take min
243 minFaceToCell(mesh.weights(), cellWeights, false);
244 Info<< " Writing face interpolation weights (0..0.5) to "
245 << cellWeights.name() << endl;
246 cellWeights.write();
247
248 if (writeFaceFields)
249 {
250 writeSurfaceField(mesh, "face_faceWeight", mesh.weights());
251 }
252 }
253
254
255 // Skewness
256 // ~~~~~~~~
257
258 if (selectedFields.found("skewness"))
259 {
260 //- Face based skewness
261 const scalarField faceSkewness
262 (
263 polyMeshTools::faceSkewness
264 (
265 mesh,
266 mesh.points(),
267 mesh.faceCentres(),
268 mesh.faceAreas(),
269 mesh.cellCentres()
270 )
271 );
272
273 //- Cell field - max of either face
274 volScalarField cellSkewness
275 (
277 (
278 "skewness",
279 mesh.time().timeName(),
280 mesh,
281 IOobject::NO_READ,
282 IOobject::AUTO_WRITE,
283 false
284 ),
285 mesh,
286 dimensionedScalar(dimless, Zero),
287 calculatedFvPatchScalarField::typeName
288 );
289 //- Take max
290 maxFaceToCell(faceSkewness, cellSkewness);
291 Info<< " Writing face skewness to " << cellSkewness.name() << endl;
292 cellSkewness.write();
293
294 if (writeFaceFields)
295 {
296 writeSurfaceField
297 (
298 mesh,
299 "face_skewness",
300 SubField<scalar>(faceSkewness, mesh.nInternalFaces())
301 );
302 }
303 }
304
305
306 // cellDeterminant
307 // ~~~~~~~~~~~~~~~
308
309 if (selectedFields.found("cellDeterminant"))
310 {
311 volScalarField cellDeterminant
312 (
314 (
315 "cellDeterminant",
316 mesh.time().timeName(),
317 mesh,
318 IOobject::NO_READ,
319 IOobject::AUTO_WRITE,
320 false
321 ),
322 mesh,
323 dimensionedScalar(dimless, Zero),
325 );
326 cellDeterminant.primitiveFieldRef() =
327 primitiveMeshTools::cellDeterminant
328 (
329 mesh,
330 mesh.geometricD(),
331 mesh.faceAreas(),
332 syncTools::getInternalOrCoupledFaces(mesh)
333 );
334 cellDeterminant.correctBoundaryConditions();
335 Info<< " Writing cell determinant to "
336 << cellDeterminant.name() << endl;
337 cellDeterminant.write();
338 }
339
340
341 // Aspect ratio
342 // ~~~~~~~~~~~~
343
344 if (selectedFields.found("aspectRatio"))
345 {
346 volScalarField aspectRatio
347 (
349 (
350 "aspectRatio",
351 mesh.time().timeName(),
352 mesh,
353 IOobject::NO_READ,
354 IOobject::AUTO_WRITE,
355 false
356 ),
357 mesh,
358 dimensionedScalar(dimless, Zero),
360 );
361
362
363 scalarField cellOpenness;
364 polyMeshTools::cellClosedness
365 (
366 mesh,
367 mesh.geometricD(),
368 mesh.faceAreas(),
369 mesh.cellVolumes(),
370 cellOpenness,
371 aspectRatio.ref()
372 );
373
374 aspectRatio.correctBoundaryConditions();
375 Info<< " Writing aspect ratio to " << aspectRatio.name() << endl;
376 aspectRatio.write();
377 }
378
379 if (selectedFields.found("cellAspectRatio"))
380 {
381 volScalarField aspectRatio
382 (
384 (
385 "cellAspectRatio",
386 mesh.time().timeName(),
387 mesh,
388 IOobject::NO_READ,
389 IOobject::AUTO_WRITE,
390 false
391 ),
392 mesh,
393 dimensionedScalar(dimless, Zero),
395 );
396
397 aspectRatio.ref().field() = cellAspectRatio(mesh);
398
399 aspectRatio.correctBoundaryConditions();
400 Info<< " Writing approximate aspect ratio to "
401 << aspectRatio.name() << endl;
402 aspectRatio.write();
403 }
404
405
406 // cell type
407 // ~~~~~~~~~
408
409 if (selectedFields.found("cellShapes"))
410 {
411 volScalarField shape
412 (
414 (
415 "cellShapes",
416 mesh.time().timeName(),
417 mesh,
418 IOobject::NO_READ,
419 IOobject::AUTO_WRITE,
420 false
421 ),
422 mesh,
423 dimensionedScalar(dimless, Zero),
425 );
426 const cellShapeList& cellShapes = mesh.cellShapes();
427 forAll(cellShapes, cellI)
428 {
429 const cellModel& model = cellShapes[cellI].model();
430 shape[cellI] = model.index();
431 }
432 shape.correctBoundaryConditions();
433 Info<< " Writing cell shape (hex, tet etc.) to " << shape.name()
434 << endl;
435 shape.write();
436 }
437
438 if (selectedFields.found("cellVolume"))
439 {
441 (
443 (
444 "cellVolume",
445 mesh.time().timeName(),
446 mesh,
447 IOobject::NO_READ,
448 IOobject::AUTO_WRITE,
449 false
450 ),
451 mesh,
452 dimensionedScalar(dimVolume, Zero),
453 calculatedFvPatchScalarField::typeName
454 );
455 V.ref() = mesh.V();
456 Info<< " Writing cell volume to " << V.name() << endl;
457 V.write();
458 }
459
460 if (selectedFields.found("cellVolumeRatio"))
461 {
462 const scalarField faceVolumeRatio
463 (
464 polyMeshTools::volRatio
465 (
466 mesh,
467 mesh.V()
468 )
469 );
470
471 volScalarField cellVolumeRatio
472 (
474 (
475 "cellVolumeRatio",
476 mesh.time().timeName(),
477 mesh,
478 IOobject::NO_READ,
479 IOobject::AUTO_WRITE,
480 false
481 ),
482 mesh,
483 dimensionedScalar(dimless, Zero),
484 calculatedFvPatchScalarField::typeName
485 );
486 //- Take min
487 minFaceToCell(faceVolumeRatio, cellVolumeRatio);
488 Info<< " Writing cell volume ratio to "
489 << cellVolumeRatio.name() << endl;
490 cellVolumeRatio.write();
491
492 if (writeFaceFields)
493 {
494 writeSurfaceField
495 (
496 mesh,
497 "face_cellVolumeRatio",
498 SubField<scalar>(faceVolumeRatio, mesh.nInternalFaces())
499 );
500 }
501 }
502
503 // minTetVolume
504 if (selectedFields.found("minTetVolume"))
505 {
506 volScalarField minTetVolume
507 (
509 (
510 "minTetVolume",
511 mesh.time().timeName(),
512 mesh,
513 IOobject::NO_READ,
514 IOobject::AUTO_WRITE,
515 false
516 ),
517 mesh,
518 dimensionedScalar("minTetVolume", dimless, GREAT),
520 );
521
522
523 const labelList& own = mesh.faceOwner();
524 const labelList& nei = mesh.faceNeighbour();
525 const pointField& p = mesh.points();
526 forAll(own, facei)
527 {
528 const face& f = mesh.faces()[facei];
529 const point& fc = mesh.faceCentres()[facei];
530
531 {
532 const point& ownCc = mesh.cellCentres()[own[facei]];
533 scalar& ownVol = minTetVolume[own[facei]];
534 forAll(f, fp)
535 {
536 scalar tetQual = tetPointRef
537 (
538 p[f[fp]],
539 p[f.nextLabel(fp)],
540 ownCc,
541 fc
542 ).quality();
543 ownVol = min(ownVol, tetQual);
544 }
545 }
546 if (mesh.isInternalFace(facei))
547 {
548 const point& neiCc = mesh.cellCentres()[nei[facei]];
549 scalar& neiVol = minTetVolume[nei[facei]];
550 forAll(f, fp)
551 {
552 scalar tetQual = tetPointRef
553 (
554 p[f[fp]],
555 p[f.nextLabel(fp)],
556 fc,
557 neiCc
558 ).quality();
559 neiVol = min(neiVol, tetQual);
560 }
561 }
562 }
563
564 minTetVolume.correctBoundaryConditions();
565 Info<< " Writing minTetVolume to " << minTetVolume.name() << endl;
566 minTetVolume.write();
567 }
568
569 // minPyrVolume
570 if (selectedFields.found("minPyrVolume"))
571 {
572 volScalarField minPyrVolume
573 (
575 (
576 "minPyrVolume",
577 mesh.time().timeName(),
578 mesh,
579 IOobject::NO_READ,
580 IOobject::AUTO_WRITE,
581 false
582 ),
583 mesh,
584 dimensionedScalar("minPyrVolume", dimless, GREAT),
586 );
587
588 // Get owner and neighbour pyr volumes
589 scalarField ownPyrVol(mesh.nFaces());
590 scalarField neiPyrVol(mesh.nInternalFaces());
591 primitiveMeshTools::facePyramidVolume
592 (
593 mesh,
594 mesh.points(),
595 mesh.cellCentres(),
596
597 ownPyrVol,
598 neiPyrVol
599 );
600
601 // Get min pyr vol per cell
602 scalarField& cellFld = minPyrVolume.ref();
603 cellFld = GREAT;
604
605 const labelUList& own = mesh.owner();
606 const labelUList& nei = mesh.neighbour();
607
608 // Internal faces
609 forAll(own, facei)
610 {
611 cellFld[own[facei]] = min(cellFld[own[facei]], ownPyrVol[facei]);
612 cellFld[nei[facei]] = min(cellFld[nei[facei]], neiPyrVol[facei]);
613 }
614
615 // Patch faces
616 for (const auto& fvp : minPyrVolume.boundaryField())
617 {
618 const labelUList& fc = fvp.patch().faceCells();
619
620 forAll(fc, i)
621 {
622 const label meshFacei = fvp.patch().start();
623 cellFld[fc[i]] = min(cellFld[fc[i]], ownPyrVol[meshFacei]);
624 }
625 }
626
627 minPyrVolume.correctBoundaryConditions();
628 Info<< " Writing minPyrVolume to " << minPyrVolume.name() << endl;
629 minPyrVolume.write();
630
631 if (writeFaceFields)
632 {
633 scalarField minFacePyrVol(neiPyrVol);
634 minFacePyrVol = min
635 (
636 minFacePyrVol,
637 SubField<scalar>(ownPyrVol, mesh.nInternalFaces())
638 );
639 writeSurfaceField(mesh, "face_minPyrVolume", minFacePyrVol);
640 }
641 }
642
643 if (selectedFields.found("cellRegion"))
644 {
645 volScalarField cellRegion
646 (
648 (
649 "cellRegion",
650 mesh.time().timeName(),
651 mesh,
652 IOobject::NO_READ,
653 IOobject::AUTO_WRITE,
654 false
655 ),
656 mesh,
657 dimensionedScalar(dimless, Zero),
658 calculatedFvPatchScalarField::typeName
659 );
660
661 regionSplit rs(mesh);
662 forAll(rs, celli)
663 {
664 cellRegion[celli] = rs[celli];
665 }
666 cellRegion.correctBoundaryConditions();
667 Info<< " Writing cell region to " << cellRegion.name() << endl;
668 cellRegion.write();
669 }
670
671 if (selectedFields.found("wallDistance"))
672 {
673 // See if wallDist.method entry in fvSchemes before calling factory
674 // method of wallDist. Have 'failing' version of wallDist::New instead?
675 const dictionary& schemesDict =
676 static_cast<const fvSchemes&>(mesh).schemesDict();
677 if (schemesDict.found("wallDist"))
678 {
679 if (schemesDict.subDict("wallDist").found("method"))
680 {
681 // Wall distance
682 volScalarField y("wallDistance", wallDist::New(mesh).y());
683 Info<< " Writing wall distance to " << y.name() << endl;
684 y.write();
685
686 // Wall-reflection vectors
687 //const volVectorField& n = wallDist::New(mesh).n();
688 //Info<< " Writing wall normal to " << n.name() << endl;
689 //n.write();
690 }
691 }
692 }
693
694 if (selectedFields.found("cellZone"))
695 {
697 (
699 (
700 "cellZone",
701 mesh.time().timeName(),
702 mesh,
703 IOobject::NO_READ,
704 IOobject::AUTO_WRITE,
705 false
706 ),
707 mesh,
708 dimensionedScalar(scalar(-1)),
709 calculatedFvPatchScalarField::typeName
710 );
711
712 const cellZoneMesh& czs = mesh.cellZones();
713 for (const auto& zone : czs)
714 {
716 }
717
718 cellZone.correctBoundaryConditions();
719 Info<< " Writing cell zoning to " << cellZone.name() << endl;
720 cellZone.write();
721 }
722 if (selectedFields.found("faceZone"))
723 {
724 // Determine for each face the zone index (scalar for ease of
725 // manipulation)
726 scalarField zoneID(mesh.nFaces(), -1);
727 const faceZoneMesh& czs = mesh.faceZones();
728 for (const auto& zone : czs)
729 {
731 }
732
733
734 // Split into internal and boundary values
736 (
738 (
739 "faceZone",
740 mesh.time().timeName(),
741 mesh,
742 IOobject::NO_READ,
743 IOobject::AUTO_WRITE,
744 false
745 ),
746 mesh,
747 dimensionedScalar(scalar(-1)),
748 calculatedFvsPatchScalarField::typeName
749 );
750
751 faceZone.primitiveFieldRef() =
752 SubField<scalar>(zoneID, mesh.nInternalFaces());
753 surfaceScalarField::Boundary& bfld = faceZone.boundaryFieldRef();
754 for (auto& pfld : bfld)
755 {
756 const fvPatch& fvp = pfld.patch();
757 pfld == SubField<scalar>(zoneID, fvp.size(), fvp.start());
758 }
759
760 //faceZone.correctBoundaryConditions();
761 Info<< " Writing face zoning to " << faceZone.name() << endl;
762 faceZone.write();
763 }
764
765 Info<< endl;
766}
scalar y
Y[inertIndex] max(0.0)
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
const Mesh & mesh() const
Return mesh.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
void correctBoundaryConditions()
Correct boundary field.
bool empty() const noexcept
True if the hash table is empty.
Definition: HashTableI.H:59
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
SubField is a Field obtained as a section of another Field, without its own allocation....
Definition: SubField.H:62
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
(Rough approximation of) cell aspect ratio
Maps a geometry to a set of cell primitives.
Definition: cellModel.H:73
label index() const noexcept
Return index of model in the model list.
Definition: cellModelI.H:37
A subset of mesh cells.
Definition: cellZone.H:65
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
virtual void write(Ostream &os) const
Write.
Definition: faceZone.C:616
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const fvPatch & patch() const
Return patch.
Definition: fvPatchField.H:362
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:71
virtual label size() const
Return size.
Definition: fvPatch.H:185
virtual label start() const
Return start label of this patch in the polyMesh face list.
Definition: fvPatch.H:179
const List< T >::subList patchSlice(const List< T > &l) const
Slice list to patch.
Definition: fvPatch.H:216
virtual const labelUList & faceCells() const
Return faceCells.
Definition: fvPatch.C:113
Selector class for finite volume differencing schemes. fvMesh is derived from fvSchemes so that all f...
Definition: fvSchemes.H:60
const fvPatch & patch() const
Return patch.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:144
A tetrahedron primitive.
Definition: tetrahedron.H:87
scalar quality() const
Return quality: Ratio of tetrahedron and circum-sphere.
Definition: tetrahedronI.H:232
label index() const noexcept
The index of this zone in the zone list.
const word & name() const noexcept
The zone name.
Base class for mesh zones.
Definition: zone.H:67
virtual void write(Ostream &os) const
Write.
Definition: zone.C:210
volScalarField & p
cellShapeList cellShapes
dynamicFvMesh & mesh
const cellShapeList & cells
const labelIOList & zoneID
Namespace for OpenFOAM.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void writeFields(const fvMesh &mesh, const wordHashSet &selectedFields, const bool writeFaceFields)
dimensionedScalar acos(const dimensionedScalar &ds)
labelList f(nPoints)
cellMask correctBoundaryConditions()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.
Foam::surfaceFields.