dynamicRefineFvMesh.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2018-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "dynamicRefineFvMesh.H"
31#include "surfaceInterpolate.H"
32#include "volFields.H"
33#include "polyTopoChange.H"
34#include "surfaceFields.H"
35#include "syncTools.H"
36#include "pointFields.H"
37#include "sigFpe.H"
38#include "cellSet.H"
39#include "HashOps.H"
40
41// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42
43namespace Foam
44{
48}
49
50// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
51
53(
54 bitSet& unrefineableCell
55) const
56{
58 {
59 unrefineableCell.clear();
60 return;
61 }
62
63 const labelList& cellLevel = meshCutter_.cellLevel();
64
65 unrefineableCell = protectedCell_;
66
67 // Get neighbouring cell level
68 labelList neiLevel(nBoundaryFaces());
69
70 for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
71 {
72 neiLevel[facei-nInternalFaces()] = cellLevel[faceOwner()[facei]];
73 }
74 syncTools::swapBoundaryFaceList(*this, neiLevel);
75
76
77 bitSet seedFace;
78
79 while (true)
80 {
81 // Pick up faces on border of protected cells
82 seedFace.reset();
83 seedFace.resize(nFaces());
84
85 for (label facei = 0; facei < nInternalFaces(); ++facei)
86 {
87 const label own = faceOwner()[facei];
88 const label nei = faceNeighbour()[facei];
89
90 if
91 (
92 // Protected owner
93 (
94 unrefineableCell.test(own)
95 && (cellLevel[nei] > cellLevel[own])
96 )
97 ||
98 // Protected neighbour
99 (
100 unrefineableCell.test(nei)
101 && (cellLevel[own] > cellLevel[nei])
102 )
103 )
104 {
105 seedFace.set(facei);
106 }
107 }
108 for (label facei = nInternalFaces(); facei < nFaces(); facei++)
109 {
110 const label own = faceOwner()[facei];
111
112 if
113 (
114 // Protected owner
115 (
116 unrefineableCell.test(own)
117 && (neiLevel[facei-nInternalFaces()] > cellLevel[own])
118 )
119 )
120 {
121 seedFace.set(facei);
122 }
123 }
124
126
127
128 // Extend unrefineableCell
129 bool hasExtended = false;
130
131 for (label facei = 0; facei < nInternalFaces(); ++facei)
132 {
133 if (seedFace.test(facei))
134 {
135 if (unrefineableCell.set(faceOwner()[facei]))
136 {
137 hasExtended = true;
138 }
139 if (unrefineableCell.set(faceNeighbour()[facei]))
140 {
141 hasExtended = true;
142 }
143 }
144 }
145 for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
146 {
147 if (seedFace.test(facei))
148 {
149 const label own = faceOwner()[facei];
150
151 if (unrefineableCell.set(own))
152 {
153 hasExtended = true;
154 }
155 }
156 }
157
158 if (!returnReduce(hasExtended, orOp<bool>()))
159 {
160 break;
161 }
162 }
163}
164
165
167{
168 const dictionary refineDict
169 (
171 (
173 (
174 "dynamicMeshDict",
175 time().constant(),
176 *this,
179 false
180 )
181 ).optionalSubDict(typeName + "Coeffs")
182 );
183
184 auto fluxVelocities = refineDict.get<List<Pair<word>>>("correctFluxes");
185
186 // Rework into hashtable.
187 correctFluxes_.resize(fluxVelocities.size());
188 for (const auto& pr : fluxVelocities)
189 {
190 correctFluxes_.insert(pr.first(), pr.second());
191 }
192
193 refineDict.readEntry("dumpLevel", dumpLevel_);
194}
195
196
198{
199 //dynamicFvMesh::mapFields(mpm);
201
202
203 // Correct old-time volumes for refined/unrefined cells. We know at this
204 // point that the points have not moved and the cells have only been split
205 // or merged. We hope that dynamicMotionSolverListFvMesh::mapFields
206 // does not use old-time volumes ...
207 {
208 const labelList& cellMap = mpm.cellMap();
209 const labelList& reverseCellMap = mpm.reverseCellMap();
210
211 // Each split cell becomes original + 7 additional
212 labelList nSubCells(mpm.nOldCells(), 0);
213
214 forAll(cellMap, celli)
215 {
216 const label oldCelli = cellMap[celli];
217 if (oldCelli >= 0 && reverseCellMap[oldCelli] >= 0)
218 {
219 // Found master cell.
220 nSubCells[oldCelli]++;
221 }
222 }
223
224 // Start off from mapped old volumes
225 scalarField correctedV0(this->V0());
226
227 // Correct any split cells
228 const auto& V = this->V();
229 forAll(cellMap, celli)
230 {
231 const label oldCelli = cellMap[celli];
232 if (oldCelli >= 0 && nSubCells[oldCelli] == 8)
233 {
234 // Found split cell. Use current volume instead of mapped
235 // old one
236 correctedV0[celli] = V[celli];
237 }
238 }
239
240 const auto& cellsFromCells = mpm.cellsFromCellsMap();
241 for (const auto& s : cellsFromCells)
242 {
243 // Reset unsplit cell
244 const label celli = s.index();
245 correctedV0[celli] = V[celli];
246 //? Or sum up old volumes?
247 //const labelList& oldCells = s.masterObjects();
248 //for (const label oldCelli : oldCells)
249 //{
250 //}
251 }
252
253 setV0().field() = correctedV0;
254 }
255
256
257 // Correct the flux for modified/added faces. All the faces which only
258 // have been renumbered will already have been handled by the mapping.
259 {
260 const labelList& faceMap = mpm.faceMap();
261 const labelList& reverseFaceMap = mpm.reverseFaceMap();
262
263 // Storage for any master faces. These will be the original faces
264 // on the coarse cell that get split into four (or rather the
265 // master face gets modified and three faces get added from the master)
266 // Estimate number of faces created
267
268 bitSet masterFaces(nFaces());
269
270 forAll(faceMap, facei)
271 {
272 const label oldFacei = faceMap[facei];
273
274 if (oldFacei >= 0)
275 {
276 const label masterFacei = reverseFaceMap[oldFacei];
277
278 if (masterFacei < 0)
279 {
281 << "Problem: should not have removed faces"
282 << " when refining."
283 << nl << "face:" << facei << endl
284 << abort(FatalError);
285 }
286 else if (masterFacei != facei)
287 {
288 masterFaces.set(masterFacei);
289 }
290 }
291 }
292
293 if (debug)
294 {
295 Pout<< "Found " << masterFaces.count() << " split faces " << endl;
296 }
297
299 (
300 lookupClass<surfaceScalarField>()
301 );
302
303 // Remove surfaceInterpolation to allow re-calculation on demand
304 // This could be done in fvMesh::updateMesh but some dynamicFvMesh
305 // might need the old interpolation fields (weights, etc).
307
308 forAllIters(fluxes, iter)
309 {
310 if (!correctFluxes_.found(iter.key()))
311 {
313 << "Cannot find surfaceScalarField " << iter.key()
314 << " in user-provided flux mapping table "
315 << correctFluxes_ << endl
316 << " The flux mapping table is used to recreate the"
317 << " flux on newly created faces." << endl
318 << " Either add the entry if it is a flux or use ("
319 << iter.key() << " none) to suppress this warning."
320 << endl;
321 continue;
322 }
323
324 const word& UName = correctFluxes_[iter.key()];
325
326 if (UName == "none")
327 {
328 continue;
329 }
330
331 surfaceScalarField& phi = *iter();
332
333 if (UName == "NaN")
334 {
335 Pout<< "Setting surfaceScalarField " << iter.key()
336 << " to NaN" << endl;
337
339
340 continue;
341 }
342
343 if (debug)
344 {
345 Pout<< "Mapping flux " << iter.key()
346 << " using interpolated flux " << UName
347 << endl;
348 }
349
350 const surfaceScalarField phiU
351 (
353 (
354 lookupObject<volVectorField>(UName)
355 )
356 & Sf()
357 );
358
359 // Recalculate new internal faces.
360 for (label facei = 0; facei < nInternalFaces(); ++facei)
361 {
362 const label oldFacei = faceMap[facei];
363
364 if (oldFacei == -1)
365 {
366 // Inflated/appended
367 phi[facei] = phiU[facei];
368 }
369 else if (reverseFaceMap[oldFacei] != facei)
370 {
371 // face-from-masterface
372 phi[facei] = phiU[facei];
373 }
374 }
375
376 // Recalculate new boundary faces.
378
379 forAll(phiBf, patchi)
380 {
381 fvsPatchScalarField& patchPhi = phiBf[patchi];
382 const fvsPatchScalarField& patchPhiU =
383 phiU.boundaryField()[patchi];
384
385 label facei = patchPhi.patch().start();
386
387 forAll(patchPhi, i)
388 {
389 const label oldFacei = faceMap[facei];
390
391 if (oldFacei == -1)
392 {
393 // Inflated/appended
394 patchPhi[i] = patchPhiU[i];
395 }
396 else if (reverseFaceMap[oldFacei] != facei)
397 {
398 // face-from-masterface
399 patchPhi[i] = patchPhiU[i];
400 }
401
402 ++facei;
403 }
404 }
405
406 // Update master faces
407 for (const label facei : masterFaces)
408 {
409 if (isInternalFace(facei))
410 {
411 phi[facei] = phiU[facei];
412 }
413 else
414 {
415 const label patchi = boundaryMesh().whichPatch(facei);
416 const label i = facei - boundaryMesh()[patchi].start();
417
418 const fvsPatchScalarField& patchPhiU =
419 phiU.boundaryField()[patchi];
420
421 fvsPatchScalarField& patchPhi = phiBf[patchi];
422
423 patchPhi[i] = patchPhiU[i];
424 }
425 }
426 }
427 }
428
429 // Correct the flux for injected faces - these are the faces which have
430 // no correspondence to the old mesh (i.e. added without a masterFace, edge
431 // or point). An example is the internal faces from hexRef8.
432 {
433 const labelList& faceMap = mpm.faceMap();
434
435 mapNewInternalFaces<scalar>(this->Sf(), this->magSf(), faceMap);
436 mapNewInternalFaces<vector>(this->Sf(), this->magSf(), faceMap);
437
438 // No oriented fields of more complex type
439 mapNewInternalFaces<sphericalTensor>(faceMap);
440 mapNewInternalFaces<symmTensor>(faceMap);
441 mapNewInternalFaces<tensor>(faceMap);
442 }
443}
444
445
446// Refines cells, maps fields and recalculates (an approximate) flux
449(
450 const labelList& cellsToRefine
451)
452{
453 // Mesh changing engine.
454 polyTopoChange meshMod(*this);
455
456 // Play refinement commands into mesh changer.
457 meshCutter_.setRefinement(cellsToRefine, meshMod);
458
459 // Create mesh (with inflation), return map from old to new mesh.
460 //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
461
462 // Clear moving flag. This is currently required since geometry calculation
463 // might get triggered when doing processor patches.
464 // (TBD: should be in changeMesh if no inflation?)
465 moving(false);
466 // Create mesh (no inflation), return map from old to new mesh.
467 autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
468
469 Info<< "Refined from "
470 << returnReduce(map().nOldCells(), sumOp<label>())
471 << " to " << globalData().nTotalCells() << " cells." << endl;
472
473 if (debug)
474 {
475 // Check map.
476 for (label facei = 0; facei < nInternalFaces(); ++facei)
477 {
478 const label oldFacei = map().faceMap()[facei];
479
480 if (oldFacei >= nInternalFaces())
481 {
483 << "New internal face:" << facei
484 << " fc:" << faceCentres()[facei]
485 << " originates from boundary oldFace:" << oldFacei
486 << abort(FatalError);
487 }
488 }
489 }
490
491 // // Remove the stored tet base points
492 // tetBasePtIsPtr_.clear();
493 // // Remove the cell tree
494 // cellTreePtr_.clear();
495
496 // Update fields
497 updateMesh(*map);
498
499
500 // Move mesh
501 /*
502 pointField newPoints;
503 if (map().hasMotionPoints())
504 {
505 newPoints = map().preMotionPoints();
506 }
507 else
508 {
509 newPoints = points();
510 }
511 movePoints(newPoints);
512 */
513
514
515
516 // Update numbering of cells/vertices.
517 meshCutter_.updateMesh(*map);
518
519 // Update numbering of protectedCell_
520 if (protectedCell_.size())
521 {
522 bitSet newProtectedCell(nCells());
523
524 forAll(newProtectedCell, celli)
525 {
526 const label oldCelli = map().cellMap()[celli];
527 if (protectedCell_.test(oldCelli))
528 {
529 newProtectedCell.set(celli);
530 }
531 }
532 protectedCell_.transfer(newProtectedCell);
533 }
534
535 // Debug: Check refinement levels (across faces only)
536 meshCutter_.checkRefinementLevels(-1, labelList());
537
538 return map;
539}
540
541
544(
545 const labelList& splitPoints
546)
547{
548 polyTopoChange meshMod(*this);
549
550 // Play refinement commands into mesh changer.
551 meshCutter_.setUnrefinement(splitPoints, meshMod);
552
553
554 // Save information on faces that will be combined
555 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
556
557 // Find the faceMidPoints on cells to be combined.
558 // for each face resulting of split of face into four store the
559 // midpoint
560 Map<label> faceToSplitPoint(3*splitPoints.size());
561
562 {
563 for (const label pointi : splitPoints)
564 {
565 const labelList& pEdges = pointEdges()[pointi];
566
567 for (const label edgei : pEdges)
568 {
569 const label otherPointi = edges()[edgei].otherVertex(pointi);
570
571 const labelList& pFaces = pointFaces()[otherPointi];
572
573 for (const label facei : pFaces)
574 {
575 faceToSplitPoint.insert(facei, otherPointi);
576 }
577 }
578 }
579 }
580
581
582 // Change mesh and generate map.
583 //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
584
585 // Clear moving flag. This is currently required since geometry calculation
586 // might get triggered when doing processor patches.
587 // (TBD: should be in changeMesh if no inflation?)
588 moving(false);
589 // Create mesh (no inflation), return map from old to new mesh.
590 autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
591
592 Info<< "Unrefined from "
593 << returnReduce(map().nOldCells(), sumOp<label>())
594 << " to " << globalData().nTotalCells() << " cells."
595 << endl;
596
597 // Update fields
598 updateMesh(*map);
599
600
601 // Move mesh
602 /*
603 pointField newPoints;
604 if (map().hasMotionPoints())
605 {
606 newPoints = map().preMotionPoints();
607 }
608 else
609 {
610 newPoints = points();
611 }
612 movePoints(newPoints);
613 */
614
615 // Correct the flux for modified faces.
616 {
617 const labelList& reversePointMap = map().reversePointMap();
618 const labelList& reverseFaceMap = map().reverseFaceMap();
619
621 (
622 lookupClass<surfaceScalarField>()
623 );
624 forAllIters(fluxes, iter)
625 {
626 if (!correctFluxes_.found(iter.key()))
627 {
629 << "Cannot find surfaceScalarField " << iter.key()
630 << " in user-provided flux mapping table "
631 << correctFluxes_ << endl
632 << " The flux mapping table is used to recreate the"
633 << " flux on newly created faces." << endl
634 << " Either add the entry if it is a flux or use ("
635 << iter.key() << " none) to suppress this warning."
636 << endl;
637 continue;
638 }
639
640 const word& UName = correctFluxes_[iter.key()];
641
642 if (UName == "none")
643 {
644 continue;
645 }
646
648 << "Mapping flux " << iter.key()
649 << " using interpolated flux " << UName
650 << endl;
651
652
653 surfaceScalarField& phi = *iter();
656
657 const surfaceScalarField phiU
658 (
660 (
661 lookupObject<volVectorField>(UName)
662 )
663 & Sf()
664 );
665
666
667 forAllConstIters(faceToSplitPoint, iter)
668 {
669 const label oldFacei = iter.key();
670 const label oldPointi = iter.val();
671
672 if (reversePointMap[oldPointi] < 0)
673 {
674 // midpoint was removed. See if face still exists.
675 const label facei = reverseFaceMap[oldFacei];
676
677 if (facei >= 0)
678 {
679 if (isInternalFace(facei))
680 {
681 phi[facei] = phiU[facei];
682 }
683 else
684 {
685 label patchi = boundaryMesh().whichPatch(facei);
686 label i = facei - boundaryMesh()[patchi].start();
687
688 const fvsPatchScalarField& patchPhiU =
689 phiU.boundaryField()[patchi];
690 fvsPatchScalarField& patchPhi = phiBf[patchi];
691 patchPhi[i] = patchPhiU[i];
692 }
693 }
694 }
695 }
696 }
697 }
698
699
700 // Update numbering of cells/vertices.
701 meshCutter_.updateMesh(*map);
702
703 // Update numbering of protectedCell_
704 if (protectedCell_.size())
705 {
706 bitSet newProtectedCell(nCells());
707
708 forAll(newProtectedCell, celli)
709 {
710 const label oldCelli = map().cellMap()[celli];
711 if (protectedCell_.test(oldCelli))
712 {
713 newProtectedCell.set(celli);
714 }
715 }
716 protectedCell_.transfer(newProtectedCell);
717 }
718
719 // Debug: Check refinement levels (across faces only)
720 meshCutter_.checkRefinementLevels(-1, labelList());
721
722 return map;
723}
724
725
728{
729 scalarField vFld(nCells(), -GREAT);
730
731 forAll(pointCells(), pointi)
732 {
733 const labelList& pCells = pointCells()[pointi];
734
735 for (const label celli : pCells)
736 {
737 vFld[celli] = max(vFld[celli], pFld[pointi]);
738 }
739 }
740 return vFld;
741}
742
743
746{
747 scalarField pFld(nPoints(), -GREAT);
748
749 forAll(pointCells(), pointi)
750 {
751 const labelList& pCells = pointCells()[pointi];
752
753 for (const label celli : pCells)
754 {
755 pFld[pointi] = max(pFld[pointi], vFld[celli]);
756 }
757 }
758 return pFld;
759}
760
761
764{
765 scalarField pFld(nPoints());
766
767 forAll(pointCells(), pointi)
768 {
769 const labelList& pCells = pointCells()[pointi];
770
771 scalar sum = 0.0;
772 for (const label celli : pCells)
773 {
774 sum += vFld[celli];
775 }
776 pFld[pointi] = sum/pCells.size();
777 }
778 return pFld;
779}
780
781
783(
784 const scalarField& fld,
785 const scalar minLevel,
786 const scalar maxLevel
787) const
788{
789 scalarField c(fld.size(), scalar(-1));
790
791 forAll(fld, i)
792 {
793 scalar err = min(fld[i]-minLevel, maxLevel-fld[i]);
794
795 if (err >= 0)
796 {
797 c[i] = err;
798 }
799 }
800 return c;
801}
802
803
805(
806 const scalar lowerRefineLevel,
807 const scalar upperRefineLevel,
808 const scalarField& vFld,
809 bitSet& candidateCell
810) const
811{
812 // Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
813 // higher more desirable to be refined).
814 scalarField cellError
815 (
816 maxPointField
817 (
818 error
819 (
820 cellToPoint(vFld),
821 lowerRefineLevel,
822 upperRefineLevel
823 )
824 )
825 );
826
827 // Mark cells that are candidates for refinement.
828 forAll(cellError, celli)
829 {
830 if (cellError[celli] > 0)
831 {
832 candidateCell.set(celli);
833 }
834 }
835}
836
837
839(
840 const label maxCells,
841 const label maxRefinement,
842 const bitSet& candidateCell
843) const
844{
845 // Every refined cell causes 7 extra cells
846 label nTotToRefine = (maxCells - globalData().nTotalCells()) / 7;
847
848 const labelList& cellLevel = meshCutter_.cellLevel();
849
850 // Mark cells that cannot be refined since they would trigger refinement
851 // of protected cells (since 2:1 cascade)
852 bitSet unrefineableCell;
853 calculateProtectedCells(unrefineableCell);
854
855 // Count current selection
856 label nLocalCandidates = candidateCell.count();
857 label nCandidates = returnReduce(nLocalCandidates, sumOp<label>());
858
859 // Collect all cells
860 DynamicList<label> candidates(nLocalCandidates);
861
862 if (nCandidates < nTotToRefine)
863 {
864 for (const label celli : candidateCell)
865 {
866 if
867 (
868 (!unrefineableCell.test(celli))
869 && cellLevel[celli] < maxRefinement
870 )
871 {
872 candidates.append(celli);
873 }
874 }
875 }
876 else
877 {
878 // Sort by error? For now just truncate.
879 for (label level = 0; level < maxRefinement; ++level)
880 {
881 for (const label celli : candidateCell)
882 {
883 if
884 (
885 (!unrefineableCell.test(celli))
886 && cellLevel[celli] == level
887 )
888 {
889 candidates.append(celli);
890 }
891 }
892
893 if (returnReduce(candidates.size(), sumOp<label>()) > nTotToRefine)
894 {
895 break;
896 }
897 }
898 }
899
900 // Guarantee 2:1 refinement after refinement
901 labelList consistentSet
902 (
903 meshCutter_.consistentRefinement
904 (
905 candidates.shrink(),
906 true // Add to set to guarantee 2:1
907 )
908 );
909
910 Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
911 << " cells for refinement out of " << globalData().nTotalCells()
912 << "." << endl;
913
914 return consistentSet;
915}
916
917
919(
920 const scalar unrefineLevel,
921 const bitSet& markedCell,
922 const scalarField& pFld
923) const
924{
925 // All points that can be unrefined
926 const labelList splitPoints(meshCutter_.getSplitPoints());
927
928
929 const labelListList& pointCells = this->pointCells();
930
931 // If we have any protected cells make sure they also are not being
932 // unrefined
933
934 bitSet protectedPoint(nPoints());
935
936 if (protectedCell_.size())
937 {
938 // Get all points on a protected cell
939 forAll(pointCells, pointi)
940 {
941 for (const label celli : pointCells[pointi])
942 {
943 if (protectedCell_.test(celli))
944 {
945 protectedPoint.set(pointi);
946 break;
947 }
948 }
949 }
950
952 (
953 *this,
954 protectedPoint,
956 0u
957 );
958
959 DebugInfo<< "From "
960 << returnReduce(protectedCell_.count(), sumOp<label>())
961 << " protected cells found "
962 << returnReduce(protectedPoint.count(), sumOp<label>())
963 << " protected points." << endl;
964 }
965
966
967 DynamicList<label> newSplitPoints(splitPoints.size());
968
969 for (const label pointi : splitPoints)
970 {
971 if (!protectedPoint[pointi] && pFld[pointi] < unrefineLevel)
972 {
973 // Check that all cells are not marked
974 bool hasMarked = false;
975
976 for (const label celli : pointCells[pointi])
977 {
978 if (markedCell.test(celli))
979 {
980 hasMarked = true;
981 break;
982 }
983 }
984
985 if (!hasMarked)
986 {
987 newSplitPoints.append(pointi);
988 }
989 }
990 }
991
992
993 newSplitPoints.shrink();
994
995 // Guarantee 2:1 refinement after unrefinement
996 labelList consistentSet
997 (
998 meshCutter_.consistentUnrefinement
999 (
1000 newSplitPoints,
1001 false
1002 )
1003 );
1004 Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
1005 << " split points out of a possible "
1006 << returnReduce(splitPoints.size(), sumOp<label>())
1007 << "." << endl;
1008
1009 return consistentSet;
1010}
1011
1012
1014(
1015 bitSet& markedCell
1016) const
1017{
1018 // Mark faces using any marked cell
1019 bitSet markedFace(nFaces());
1020
1021 for (const label celli : markedCell)
1022 {
1023 markedFace.set(cells()[celli]); // set multiple faces
1024 }
1025
1026 syncTools::syncFaceList(*this, markedFace, orEqOp<unsigned int>());
1027
1028 // Update cells using any markedFace
1029 for (label facei = 0; facei < nInternalFaces(); ++facei)
1030 {
1031 if (markedFace.test(facei))
1032 {
1033 markedCell.set(faceOwner()[facei]);
1034 markedCell.set(faceNeighbour()[facei]);
1035 }
1036 }
1037 for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
1038 {
1039 if (markedFace.test(facei))
1040 {
1041 markedCell.set(faceOwner()[facei]);
1042 }
1043 }
1044}
1045
1046
1048(
1049 bitSet& protectedCell
1050) const
1051{
1052 const labelList& cellLevel = meshCutter_.cellLevel();
1053 const labelList& pointLevel = meshCutter_.pointLevel();
1054
1055 labelList nAnchorPoints(nCells(), Zero);
1056
1057 forAll(pointLevel, pointi)
1058 {
1059 const labelList& pCells = pointCells(pointi);
1060
1061 for (const label celli : pCells)
1062 {
1063 if (pointLevel[pointi] <= cellLevel[celli])
1064 {
1065 // Check if cell has already 8 anchor points -> protect cell
1066 if (nAnchorPoints[celli] == 8)
1067 {
1068 protectedCell.set(celli);
1069 }
1070
1071 if (!protectedCell.test(celli))
1072 {
1073 ++nAnchorPoints[celli];
1074 }
1075 }
1076 }
1077 }
1078
1079
1080 forAll(protectedCell, celli)
1081 {
1082 if (nAnchorPoints[celli] != 8)
1083 {
1084 protectedCell.set(celli);
1085 }
1086 }
1087}
1088
1089
1090// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1091
1093(
1094 const IOobject& io,
1095 const bool doInit
1096)
1097:
1098 //dynamicFvMesh(io, doInit),
1100 meshCutter_(*this)
1101{
1102 if (doInit)
1103 {
1104 init(false); // do not initialise lower levels
1105 }
1106}
1107
1108
1109bool Foam::dynamicRefineFvMesh::init(const bool doInit)
1110{
1111 if (doInit)
1112 {
1113 //dynamicFvMesh::init(doInit);
1114 // Note: allow zero motion solvers
1116 }
1117
1118 protectedCell_.setSize(nCells());
1119 nRefinementIterations_ = 0;
1120 dumpLevel_ = false;
1121
1122 // Read static part of dictionary
1123 readDict();
1124
1125
1126 const labelList& cellLevel = meshCutter_.cellLevel();
1127 const labelList& pointLevel = meshCutter_.pointLevel();
1128
1129 // Set cells that should not be refined.
1130 // This is currently any cell which does not have 8 anchor points or
1131 // uses any face which does not have 4 anchor points.
1132 // Note: do not use cellPoint addressing
1133
1134 // Count number of points <= cellLevel
1135 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1136
1137 labelList nAnchors(nCells(), Zero);
1138
1139 forAll(pointCells(), pointi)
1140 {
1141 const labelList& pCells = pointCells()[pointi];
1142
1143 for (const label celli : pCells)
1144 {
1145 if (!protectedCell_.test(celli))
1146 {
1147 if (pointLevel[pointi] <= cellLevel[celli])
1148 {
1149 ++nAnchors[celli];
1150
1151 if (nAnchors[celli] > 8)
1152 {
1153 protectedCell_.set(celli);
1154 }
1155 }
1156 }
1157 }
1158 }
1159
1160
1161 // Count number of points <= faceLevel
1162 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1163 // Bit tricky since proc face might be one more refined than the owner since
1164 // the coupled one is refined.
1165
1166 {
1167 labelList neiLevel(nFaces());
1168
1169 for (label facei = 0; facei < nInternalFaces(); ++facei)
1170 {
1171 neiLevel[facei] = cellLevel[faceNeighbour()[facei]];
1172 }
1173 for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
1174 {
1175 neiLevel[facei] = cellLevel[faceOwner()[facei]];
1176 }
1177 syncTools::swapFaceList(*this, neiLevel);
1178
1179
1180 bitSet protectedFace(nFaces());
1181
1182 forAll(faceOwner(), facei)
1183 {
1184 const label faceLevel = max
1185 (
1186 cellLevel[faceOwner()[facei]],
1187 neiLevel[facei]
1188 );
1189
1190 const face& f = faces()[facei];
1191
1192 label nAnchors = 0;
1193
1194 for (const label pointi : f)
1195 {
1196 if (pointLevel[pointi] <= faceLevel)
1197 {
1198 ++nAnchors;
1199
1200 if (nAnchors > 4)
1201 {
1202 protectedFace.set(facei);
1203 break;
1204 }
1205 }
1206 }
1207 }
1208
1209 syncTools::syncFaceList(*this, protectedFace, orEqOp<unsigned int>());
1210
1211 for (label facei = 0; facei < nInternalFaces(); ++facei)
1212 {
1213 if (protectedFace.test(facei))
1214 {
1215 protectedCell_.set(faceOwner()[facei]);
1216 protectedCell_.set(faceNeighbour()[facei]);
1217 }
1218 }
1219 for (label facei = nInternalFaces(); facei < nFaces(); ++facei)
1220 {
1221 if (protectedFace.test(facei))
1222 {
1223 protectedCell_.set(faceOwner()[facei]);
1224 }
1225 }
1226
1227 // Also protect any cells that are less than hex
1228 forAll(cells(), celli)
1229 {
1230 const cell& cFaces = cells()[celli];
1231
1232 if (cFaces.size() < 6)
1233 {
1234 protectedCell_.set(celli);
1235 }
1236 else
1237 {
1238 for (const label cfacei : cFaces)
1239 {
1240 if (faces()[cfacei].size() < 4)
1241 {
1242 protectedCell_.set(celli);
1243 break;
1244 }
1245 }
1246 }
1247 }
1248
1249 // Check cells for 8 corner points
1250 checkEightAnchorPoints(protectedCell_);
1251 }
1252
1253 if (!returnReduce(protectedCell_.any(), orOp<bool>()))
1254 {
1255 protectedCell_.clear();
1256 }
1257 else
1258 {
1259 cellSet protectedCells
1260 (
1261 *this,
1262 "protectedCells",
1263 HashSetOps::used(protectedCell_)
1264 );
1265
1266 Info<< "Detected "
1267 << returnReduce(protectedCells.size(), sumOp<label>())
1268 << " cells that are protected from refinement."
1269 << " Writing these to cellSet "
1270 << protectedCells.name()
1271 << "." << endl;
1272
1273 protectedCells.write();
1274 }
1275
1276 return true;
1277}
1278
1279
1280// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1281
1283{
1284 // Re-read dictionary. Chosen since usually -small so trivial amount
1285 // of time compared to actual refinement. Also very useful to be able
1286 // to modify on-the-fly.
1287 dictionary refineDict
1288 (
1290 (
1291 IOobject
1292 (
1293 "dynamicMeshDict",
1294 time().constant(),
1295 *this,
1298 false
1299 )
1300 ).optionalSubDict(typeName + "Coeffs")
1301 );
1302
1303 const label refineInterval = refineDict.get<label>("refineInterval");
1304
1305 bool hasChanged = false;
1306
1307 if (refineInterval == 0)
1308 {
1309 topoChanging(hasChanged);
1310
1311 return false;
1312 }
1313 else if (refineInterval < 0)
1314 {
1316 << "Illegal refineInterval " << refineInterval << nl
1317 << "The refineInterval setting in the dynamicMeshDict should"
1318 << " be >= 1." << nl
1319 << exit(FatalError);
1320 }
1321
1322
1323 // Note: cannot refine at time 0 since no V0 present since mesh not
1324 // moved yet.
1325
1326 if (time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0)
1327 {
1328 const label maxCells = refineDict.get<label>("maxCells");
1329
1330 if (maxCells <= 0)
1331 {
1333 << "Illegal maximum number of cells " << maxCells << nl
1334 << "The maxCells setting in the dynamicMeshDict should"
1335 << " be > 0." << nl
1336 << exit(FatalError);
1337 }
1338
1339 const label maxRefinement = refineDict.get<label>("maxRefinement");
1340
1341 if (maxRefinement <= 0)
1342 {
1344 << "Illegal maximum refinement level " << maxRefinement << nl
1345 << "The maxCells setting in the dynamicMeshDict should"
1346 << " be > 0." << nl
1347 << exit(FatalError);
1348 }
1349
1350 const word fieldName(refineDict.get<word>("field"));
1351
1352 const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
1353
1354 const scalar lowerRefineLevel =
1355 refineDict.get<scalar>("lowerRefineLevel");
1356 const scalar upperRefineLevel =
1357 refineDict.get<scalar>("upperRefineLevel");
1358 const scalar unrefineLevel = refineDict.getOrDefault<scalar>
1359 (
1360 "unrefineLevel",
1361 GREAT
1362 );
1363 const label nBufferLayers = refineDict.get<label>("nBufferLayers");
1364
1365 // Cells marked for refinement or otherwise protected from unrefinement.
1366 bitSet refineCell(nCells());
1367
1368 // Determine candidates for refinement (looking at field only)
1369 selectRefineCandidates
1370 (
1371 lowerRefineLevel,
1372 upperRefineLevel,
1373 vFld,
1375 );
1376
1377 if (globalData().nTotalCells() < maxCells)
1378 {
1379 // Select subset of candidates. Take into account max allowable
1380 // cells, refinement level, protected cells.
1381 labelList cellsToRefine
1382 (
1383 selectRefineCells
1384 (
1385 maxCells,
1386 maxRefinement,
1388 )
1389 );
1390
1391 const label nCellsToRefine = returnReduce
1392 (
1393 cellsToRefine.size(), sumOp<label>()
1394 );
1395
1396 if (nCellsToRefine > 0)
1397 {
1398 // Refine/update mesh and map fields
1399 autoPtr<mapPolyMesh> map = refine(cellsToRefine);
1400
1401 // Update refineCell. Note that some of the marked ones have
1402 // not been refined due to constraints.
1403 {
1404 const labelList& cellMap = map().cellMap();
1405 const labelList& reverseCellMap = map().reverseCellMap();
1406
1407 bitSet newRefineCell(cellMap.size());
1408
1409 forAll(cellMap, celli)
1410 {
1411 const label oldCelli = cellMap[celli];
1412
1413 if
1414 (
1415 (oldCelli < 0)
1416 || (reverseCellMap[oldCelli] != celli)
1417 || (refineCell.test(oldCelli))
1418 )
1419 {
1420 newRefineCell.set(celli);
1421 }
1422 }
1423 refineCell.transfer(newRefineCell);
1424 }
1425
1426 // Extend with a buffer layer to prevent neighbouring points
1427 // being unrefined.
1428 for (label i = 0; i < nBufferLayers; ++i)
1429 {
1430 extendMarkedCells(refineCell);
1431 }
1432
1433 hasChanged = true;
1434 }
1435 }
1436
1437
1438 {
1439 // Select unrefineable points that are not marked in refineCell
1440 labelList pointsToUnrefine
1441 (
1442 selectUnrefinePoints
1443 (
1444 unrefineLevel,
1445 refineCell,
1446 maxCellField(vFld)
1447 )
1448 );
1449
1450 const label nSplitPoints = returnReduce
1451 (
1452 pointsToUnrefine.size(),
1453 sumOp<label>()
1454 );
1455
1456 if (nSplitPoints > 0)
1457 {
1458 // Refine/update mesh
1459 unrefine(pointsToUnrefine);
1460
1461 hasChanged = true;
1462 }
1463 }
1464
1465
1466 if ((nRefinementIterations_ % 10) == 0)
1467 {
1468 // Compact refinement history occasionally (how often?).
1469 // Unrefinement causes holes in the refinementHistory.
1470 const_cast<refinementHistory&>(meshCutter().history()).compact();
1471 }
1472 nRefinementIterations_++;
1473 }
1474
1475 topoChanging(hasChanged);
1476 if (hasChanged)
1477 {
1478 // Reset moving flag (if any). If not using inflation we'll not move,
1479 // if are using inflation any follow on movePoints will set it.
1480 moving(false);
1481 }
1482
1483 return hasChanged;
1484}
1485
1486
1488{
1489 bool hasChanged = updateTopology();
1490 // Do any mesh motion (resets mesh.moving() if it does any mesh motion)
1491 hasChanged = dynamicMotionSolverListFvMesh::update() && hasChanged;
1492
1493 return hasChanged;
1494}
1495
1496
1498(
1499 IOstreamOption streamOpt,
1500 const bool valid
1501) const
1502{
1503 // Force refinement data to go to the current time directory.
1504 const_cast<hexRef8&>(meshCutter_).setInstance(time().timeName());
1505
1506 bool writeOk =
1507 (
1508 //dynamicFvMesh::writeObject(streamOpt, valid)
1510 && meshCutter_.write(valid)
1511 );
1512
1513 if (dumpLevel_)
1514 {
1515 volScalarField scalarCellLevel
1516 (
1517 IOobject
1518 (
1519 "cellLevel",
1520 time().timeName(),
1521 *this,
1524 false
1525 ),
1526 *this,
1528 );
1529
1530 const labelList& cellLevel = meshCutter_.cellLevel();
1531
1532 forAll(cellLevel, celli)
1533 {
1534 scalarCellLevel[celli] = cellLevel[celli];
1535 }
1536
1537 writeOk = writeOk && scalarCellLevel.write();
1538 }
1539
1540 return writeOk;
1541}
1542
1543
1544// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
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))
surfaceScalarField & phi
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Boundary & boundaryFieldRef(const bool updateAccessTime=true)
Return a reference to the boundary field.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
@ MUST_READ_IF_MODIFIED
Definition: IOobject.H:180
The IOstreamOption is a simple container for options an IOstream can normally have.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:330
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
bool empty() const noexcept
True if the list is empty (ie, size() is zero).
Definition: PackedListI.H:384
void resize(const label numElem, const unsigned int val=0u)
Reset addressable list size, does not shrink the allocated size.
Definition: PackedListI.H:409
void clear()
Clear the list, i.e. set addressable size to zero.
Definition: PackedListI.H:512
void reset()
Clear all bits but do not adjust the addressable size.
Definition: PackedListI.H:505
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
unsigned int count(const bool on=true) const
Count number of bits set.
Definition: bitSetI.H:500
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A collection of cell labels.
Definition: cellSet.H:54
A topoSetPointSource to select all the points from given cellSet(s).
Definition: cellToPoint.H:176
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
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Abstract base class for geometry and/or topology changing fvMesh.
Definition: dynamicFvMesh.H:81
Dynamic mesh able to handle multiple motion solvers. NOTE: If the word entry "solvers" is not found i...
virtual bool update()
Update the mesh for both mesh motion and topology change.
A fvMesh with built-in refinement.
void readDict()
Read the projection parameters from dictionary.
bool updateTopology()
Update topology (refinement, unrefinement)
hexRef8 meshCutter_
Mesh cutting engine.
virtual bool init(const bool doInit)
Initialise all non-demand-driven data.
virtual autoPtr< mapPolyMesh > refine(const labelList &)
Refine cells. Update mesh and fields.
scalarField maxPointField(const scalarField &) const
Get per cell max of connected point.
void calculateProtectedCells(bitSet &unrefineableCell) const
Calculate cells that cannot be refined since would trigger.
virtual labelList selectRefineCells(const label maxCells, const label maxRefinement, const bitSet &candidateCell) const
Subset candidate cells for refinement.
bitSet protectedCell_
Protected cells (usually since not hexes)
virtual bool writeObject(IOstreamOption streamOpt, const bool valid) const
Write using stream options.
virtual bool update()
Update the mesh for both mesh motion and topology change.
virtual void selectRefineCandidates(const scalar lowerRefineLevel, const scalar upperRefineLevel, const scalarField &vFld, bitSet &candidateCell) const
Select candidate cells for refinement.
scalarField maxCellField(const volScalarField &) const
Get point max of connected cell.
virtual labelList selectUnrefinePoints(const scalar unrefineLevel, const bitSet &markedCell, const scalarField &pFld) const
Select points that can be unrefined.
void extendMarkedCells(bitSet &markedCell) const
Extend markedCell with cell-face-cell.
virtual autoPtr< mapPolyMesh > unrefine(const labelList &)
Unrefine cells. Gets passed in centre points of cells to combine.
void checkEightAnchorPoints(bitSet &protectedCell) const
Check all cells have 8 anchor points.
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:77
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Maps input fields from local mesh to secondary mesh at runtime.
Definition: mapFields.H:225
virtual label start() const
Return start label of this patch in the polyMesh face list.
Definition: fvPatch.H:179
const fvPatch & patch() const
Return patch.
Refinement of (split) hexes using polyTopoChange.
Definition: hexRef8.H:68
const labelIOList & cellLevel() const
Definition: hexRef8.H:397
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
label nOldCells() const
Number of old cells.
Definition: mapPolyMesh.H:387
const labelList & faceMap() const
Old face map.
Definition: mapPolyMesh.H:410
const labelList & reverseCellMap() const
Reverse cell map.
Definition: mapPolyMesh.H:532
const labelList & cellMap() const
Old cell map.
Definition: mapPolyMesh.H:435
const List< objectMap > & cellsFromCellsMap() const
Cells originating from cells.
Definition: mapPolyMesh.H:459
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:501
Cuts (splits) cells.
Definition: meshCutter.H:141
constant condensation/saturation model.
Smooth ATC in cells having a point to a set of patches supplied by type.
Definition: pointCells.H:59
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
Direct mesh changes based on v1.3 polyTopoChange syntax.
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
label nInternalFaces() const noexcept
Number of internal faces.
label nFaces() const noexcept
Number of mesh faces.
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:57
All refinement history. Used in unrefinement.
virtual bool write(const bool valid=true) const
Write using setting from DB.
static void fillNan(UList< scalar > &list)
Fill data block with NaN values.
Definition: sigFpe.C:225
void clearOut()
Clear all geometry and addressing.
static void swapFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled face values. Uses eqOp.
Definition: syncTools.H:478
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
label nPoints
const cellShapeList & cells
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
word timeName
Definition: getTimeIndex.H:3
#define DebugInfo
Report an information message using Foam::Info.
#define WarningInFunction
Report a warning using Foam::Warning.
labelHashSet used(const bitSet &select)
Convert a bitset to a labelHashSet of the indices used.
Definition: HashOps.C:35
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
Definition: List.H:66
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
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
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label timeIndex
Definition: getTimeIndex.H:30
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:260
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
Foam::surfaceFields.
const word UName(propsDict.getOrDefault< word >("U", "U"))