diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..273b4dad2085f94fcdb458d476d3ca83ac62eee9
--- /dev/null
+++ b/README.txt
@@ -0,0 +1,47 @@
+Notes from merging Integration-TUD-addOns and plus.
+
+0)
+- Changed
+dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+to created internal faces out-of-nothing.
+- Added the mapNewInternalFaces to dynamicRefineFvMesh
+- unset FOAM_SETNAN, FOAM_SIGFPE and run
+
+ testCases/testAMRandLoadBalancing/damBreakWithObstacle
+
+- However when writing V0 we noticed a nan in the V0 field.
+ (this doesn't get written anymore but it might indicate a bug)
+
+
+1) Tried moving additional mapping of surface fields to
+
+ virtual dynamicRefineFvMesh::mapFields(const mapPolyMesh&);
+
+so it would get called from the fvMesh::updateMesh. However
+this mapping (explicitly) gets done using unadapted addressing
+and it causes the cell addressing to be wrong:
+
+
+[1] --> FOAM FATAL ERROR:
+[1] index 826015320 out of range 0 ... 71079
+[1]
+[1] From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = Foam::cell; Foam::label = int]
+[1] in file /home/preston2/mattijs/OpenFOAM/work/OpenFOAM-plus.integration-TUD/src/OpenFOAM/lnInclude/UListI.H at line 106.
+
+
+From dynamicRefineFvMesh::mapNewInternalFaces : (I think)
+
+ cell faceOwner = this->cells()[owner[facei]];
+
+So this adaptation should be done in a separate pass
+
+
+2) Mapping new internal faces:
+- currently done by averaging the values of 'properly' mapped
+ faces from owner and neighbour.
+- this would not work if a cell gets split into 3x3x3
+- instead this should be done by some geometric interpolation
+ from point values?
+- have selection mechanism instead or use oriented flag on surfaceFields
+
+3) Current averaging is valid only for internal faces & only scalar/vector
diff --git a/src/dynamicFvMesh/Make/files b/src/dynamicFvMesh/Make/files
index b8432c07c905108747f9b94eb4cffefd26d1e0bd..bb04e8ac64705ac6dc159ca970a4bc5dc7320ca9 100644
--- a/src/dynamicFvMesh/Make/files
+++ b/src/dynamicFvMesh/Make/files
@@ -5,6 +5,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
dynamicMultiMotionSolverFvMesh/dynamicMultiMotionSolverFvMesh.C
dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
dynamicRefineFvMesh/dynamicRefineFvMesh.C
+dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C
diff --git a/src/dynamicFvMesh/Make/options b/src/dynamicFvMesh/Make/options
index 3540c13bfc8d77a85b5642874b6806561cfdacd8..718fde7d5c5de0e42db7e1a715730dd1bcc68f47 100644
--- a/src/dynamicFvMesh/Make/options
+++ b/src/dynamicFvMesh/Make/options
@@ -1,10 +1,13 @@
EXE_INC = \
+ -g -DFULLDEBUG -O0 -g \
+ -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
+ -ldecompositionMethods \
-lmeshTools \
-ldynamicMesh \
-lfiniteVolume
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict
new file mode 100644
index 0000000000000000000000000000000000000000..08d1ef65388aa50771de1f3b78493a2d4bbdad0b
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 2.1.1 |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object balanceParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method ptscotch;
+
+constraints
+{
+ refinementHistory
+ {
+ type refinementHistory;
+ }
+}
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..a8df1cb7efa27c0931063127ac6dd3d65dd2cefc
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict
@@ -0,0 +1,125 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "constant";
+ object dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh dynamicRefineBalancedFvMesh;
+
+refinementControls
+{
+ enableRefinementControl true;
+
+ fields // must be scalarFields
+ (
+ //alpha (min max refineLevel)
+ alpha (0.01 0.99 2) // refine cells where alpha in [0.01:0.99] with maximal 2 refinement layers
+ );
+
+ interface // must be a scalarField (only one dictionary!)
+ (
+ alpha // refine interface (found based on snGrad of alpha > 0.1)
+ {
+ innerRefLayers 2; // describes how many cell layers inside phase alpha are to be refined
+ outerRefLayers 5; // describes how many cell layers outside phase alpha are to be refined
+
+ // optional settings:
+ maxRefineLevel 4; // max refinement layers; Default: maxRefinement from dynamicRefineFvMeshCoeffs is used
+ // to get slower than 2:1 refinement; add #nAddLayers between each refinement level at the interface
+ nAddLayers 1; //Default: 0
+ }
+ );
+
+ gradients // must be scalars
+ (
+ // arguments as in 'fields'
+ // min/max values are based on mag(fvc::grad(volScalarField)) * cellVolume
+ T (0.01 10 1)
+ );
+
+ curls // must be vectors
+ (
+ // arguments as in 'fields'
+ // min/max values are based on mag(fvc::curl(volVectorField))
+ U (0.5 1 2)
+ );
+
+ regions
+ (
+ boxToCell
+ {
+ minLevel 1;
+
+ box (-1 0.001 0.002)(1 0.005 0.003);
+ }
+
+ );
+}
+
+dynamicRefineFvMeshCoeffs
+{
+ // Extra entries for balancing
+ enableBalancing true;
+ allowableImbalance 0.15;
+
+ // How often to refine
+ refineInterval 10;
+
+ // Field to be refinement on (set it to 'internalRefinementField' to use the
+ // refinementControls dictionary entries above)
+ field internalRefinementField;
+
+ // Refine field inbetween lower..upper
+ lowerRefineLevel 0.5; // do not change
+ upperRefineLevel 3.5; // maxRefinement+0.5
+
+ // If value < unrefineLevel unrefine
+ unrefineLevel -0.5; // do not change
+
+ // Have slower than 2:1 refinement
+ nBufferLayers 4;
+
+ // Refine cells only up to maxRefinement levels
+ maxRefinement 3;
+
+ // Stop refinement if maxCells reached
+ maxCells 200000;
+
+ // Flux field and corresponding velocity field. Fluxes on changed
+ // faces get recalculated by interpolating the velocity. Use 'none'
+ // on surfaceScalarFields that do not need to be reinterpolated.
+ correctFluxes
+ (
+ (phi Urel)
+ (phiAbs U)
+ (phiAbs_0 U_0)
+ (nHatf none)
+ (rho*phi none)
+ (ghf none)
+ );
+
+ // List of non-flux surface<Type>Fields to be mapped
+ // only for new internal faces (AMR refine)
+ mapSurfaceFields
+ (
+ Uf
+ Uf_0
+ );
+
+ // Write the refinement level as a volScalarField
+ dumpLevel true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..d5098ffe6805cbe0c7b58dd3c47c91423036ca70
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
@@ -0,0 +1,843 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 Tyler Voskuilen
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is a derivative work of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "dynamicRefineBalancedFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceInterpolate.H"
+#include "volFields.H"
+#include "polyTopoChange.H"
+#include "surfaceFields.H"
+#include "syncTools.H"
+#include "pointFields.H"
+#include "fvCFD.H"
+#include "volPointInterpolation.H"
+#include "pointMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(dynamicRefineBalancedFvMesh, 0);
+ addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineBalancedFvMesh, IOobject);
+}
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+Foam::label Foam::dynamicRefineBalancedFvMesh::topParentID(label p)
+{
+ label nextP = meshCutter().history().splitCells()[p].parent_;
+ if( nextP < 0 )
+ {
+ return p;
+ }
+ else
+ {
+ return topParentID(nextP);
+ }
+}
+
+Foam::List<Foam::scalar> Foam::dynamicRefineBalancedFvMesh::readRefinementPoints()
+{
+ dictionary refineDict
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "dynamicMeshDict",
+ time().constant(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE,
+ false
+ )
+ ).subDict("dynamicRefineFvMeshCoeffs")
+ );
+
+ List<scalar> refData(4, scalar(0));
+
+ refData[0] = readScalar(refineDict.lookup("unrefineLevel"));
+ refData[1] = readScalar(refineDict.lookup("lowerRefineLevel"));
+ refData[2] = readScalar(refineDict.lookup("refineInterval"));
+ refData[3] = readScalar(refineDict.lookup("maxRefinement"));
+
+ return refData;
+}
+
+void Foam::dynamicRefineBalancedFvMesh::updateRefinementField()
+{
+ Info<< "Calculating internal refinement field" << endl;
+
+ volScalarField& intRefFld = *internalRefinementFieldPtr_;
+ volScalarField& targetFld = *targetLevelPtr_;
+ volScalarField& currFld = *isLevelPtr_;
+
+ // Set the internal refinement field to zero to start with
+ intRefFld = dimensionedScalar("zero",dimless,0.0);
+
+ // Get the cell level field from dynamicRefineFvMesh
+ const labelList& cellLevel = meshCutter().cellLevel();
+
+ // Read the points at which refinement and unrefinement occur from the
+ // dynamicMeshDict entries
+ List<scalar> refinePoints = readRefinementPoints();
+
+ // init list for refine/unrefine field (-1=unrefine, 1=refine, 0=do nothing)
+ labelList markRefineCells (this->nCells(), 0);
+
+ // init list for target refinement level per cell
+ labelList targetLevel (this->nCells(), 0);
+
+ // First fields
+ List<word> fieldNames = fields_.toc();
+ Field<scalar> refFld(nCells(),0.0);
+
+ forAll(fieldNames, i)
+ {
+ word fldName = fieldNames[i];
+ scalar minValue = fields_[fldName][0];
+ scalar maxValue = fields_[fldName][1];
+ label refineLevel = static_cast<label>(fields_[fldName][2]);
+
+ const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+
+ // Limit the value of refFld based on its max level
+ forAll(fld, cellI)
+ {
+ if ((fld[cellI] >= minValue) && (fld[cellI] <= maxValue))
+ {
+ // increase targetLevel up to refineLevel
+ // BUT: do not decrease if cell already marked for higher refinement level by previous criterion
+ targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+ }
+ }
+ }
+
+ // Then gradients
+ List<word> gradFieldNames = gradFields_.toc();
+
+ Field<scalar> cubeRtV = Foam::pow(this->V(),1.0/3.0);
+
+ forAll(gradFieldNames, i)
+ {
+ word fldName = gradFieldNames[i];
+ scalar minValue = gradFields_[fldName][0];
+ scalar maxValue = gradFields_[fldName][1];
+ label refineLevel = static_cast<label>(gradFields_[fldName][2]);
+
+ const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+
+ refFld = mag(fvc::grad(fld)) * cubeRtV;
+
+ // Limit the value of refFld based on its max level
+ forAll(refFld, cellI)
+ {
+ if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
+ {
+ targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+ }
+ }
+ }
+
+ // Then curls
+ List<word> curlFieldNames = curlFields_.toc();
+
+
+ forAll(curlFieldNames, i)
+ {
+ word fldName = curlFieldNames[i];
+ scalar minValue = curlFields_[fldName][0];
+ scalar maxValue = curlFields_[fldName][1];
+ label refineLevel = static_cast<label>(curlFields_[fldName][2]);
+
+ const volVectorField& fld = this->lookupObject<volVectorField>(fldName);
+
+ refFld = mag(fvc::curl(fld));
+
+ // Limit the value of refFld based on its max level
+ forAll(refFld, cellI)
+ {
+ if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
+ {
+ targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+ }
+ }
+ }
+
+ // At the interface, assumed to be always the maximum refinement level
+ List<word> interfaceRefineField = interface_.toc();
+ forAll(interfaceRefineField, i)
+ {
+ word fldName = interfaceRefineField[i];
+ // read region of maximum refinement levels inside and outside of interface indicator field
+ // (does not need to be alpha, can also be a concentration field)
+
+ scalar innerRefLayers = readScalar(interface_[fldName].lookup("innerRefLayers"));
+ scalar outerRefLayers = readScalar(interface_[fldName].lookup("outerRefLayers"));
+
+ scalar nAddLayers(0);
+ if (interface_[fldName].found("nAddLayers"))
+ {
+ nAddLayers = readScalar(interface_[fldName].lookup("nAddLayers"));
+ }
+
+ label refineLevel = refinePoints[3];
+
+ if (interface_[fldName].found("maxRefineLevel"))
+ {
+ refineLevel = readScalar(interface_[fldName].lookup("maxRefineLevel"));
+
+ // to avoid user input mistakes, limit the value with the maximum allowed
+ refineLevel = min(refinePoints[3], refineLevel);
+ }
+
+ const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
+
+ volScalarField isInterface(intRefFld * 0.0);
+ isInterface = dimensionedScalar("isInterface",dimless,0.0);
+
+ surfaceScalarField deltaAlpha = mag(fvc::snGrad(fld) / this->deltaCoeffs());
+
+ const labelUList& owner = this->owner();
+ const labelUList& neighbour = this->neighbour();
+
+ forAll(deltaAlpha, faceI)
+ {
+ //TODO: add a min max value of the specified field to be marked for refinement
+ // NOW: hard-coded method snGradAlpha: checks mag(alpha_N - alpha_P) > 0.1 ? 1:0
+ if (deltaAlpha[faceI] > 0.1) // currently a fixed prescribed value; should be read from dÃctionary
+ {
+ label own = owner[faceI];
+ label nei = neighbour[faceI];
+
+ // set isInterface field to one
+ isInterface[own] = 1.0;
+ isInterface[nei] = 1.0;
+ }
+ }
+
+ // assumed fld=0.5*(fldMax+fldMin) defines the interface
+ dimensionedScalar fldInterfaceValue(0.5*(max(fld)+min(fld)));
+
+ //-DD: old implementation based on face interpolation
+ // which results in slower transport in diagonal direction
+ // add inner refinement layers
+ //for(label i=0; i < innerRefLayers; i++)
+ //{
+ // isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fld - fldInterfaceValue));
+ // isInterface = neg(- isInterface);
+ //}
+ //
+ // add outer refinement layers
+ //for(label i=0; i < outerRefLayers; i++)
+ //{
+ // isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fldInterfaceValue - fld));
+ // isInterface = neg(- isInterface);
+ //}
+ //
+ //forAll(isInterface, cellI)
+ //{
+ // if (isInterface[cellI] > 0.5)
+ // {
+ // targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+ // }
+ //}
+
+ //-DD: new version using volPointInterpolation (direction independent buffer layer)
+ const volPointInterpolation& pInterp = volPointInterpolation::New(*this);
+ const fvMesh& mesh = fld.mesh();
+
+ // add inner refinement layers
+ for(label i=0; i < innerRefLayers; i++)
+ {
+ volScalarField markInner(isInterface*pos(fld - fldInterfaceValue));
+ pointScalarField markLayerP(pInterp.interpolate(markInner));
+
+ forAll(mesh.C(), cellI)
+ {
+ scalar sum = 0.;
+ label nPoints = 0;
+
+ forAll(mesh.cellPoints()[cellI], pointI)
+ {
+ sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
+ nPoints++;
+ }
+ if (nPoints > 0)
+ {
+ sum /= nPoints;
+ }
+ isInterface[cellI] += sum;
+ }
+ }
+ isInterface = pos(isInterface - SMALL);
+
+ // add outer refinement layers
+ for(label i=0; i < outerRefLayers; i++)
+ {
+ volScalarField markOuter(isInterface*pos(fldInterfaceValue - fld));
+ pointScalarField markLayerP(pInterp.interpolate(markOuter));
+
+ forAll(mesh.C(), cellI)
+ {
+ scalar sum = 0.;
+ label nPoints = 0;
+
+ forAll(mesh.cellPoints()[cellI], pointI)
+ {
+ sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
+ nPoints++;
+ }
+ if (nPoints > 0)
+ {
+ sum /= nPoints;
+ }
+ isInterface[cellI] += sum;
+ }
+ }
+ isInterface = pos(isInterface - SMALL);
+
+ forAll(isInterface, cellI)
+ {
+ if (isInterface[cellI] > 0.5)
+ {
+ targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
+ }
+ }
+
+ //-DD: old implementation based on face interpolation
+ // which results in slower transport in diagonal direction
+ // expand additional layers if specified:
+ if (nAddLayers > 0)
+ {
+ // loop over additional layers
+ for(label i=1; i < refineLevel; i++)
+ {
+ // #nAddLayers cells per additional level
+ for(label j=0; j < nAddLayers; j++)
+ {
+ isInterface += neg(- fvc::average(fvc::interpolate(isInterface)));
+ isInterface = neg(- isInterface);
+ }
+
+ forAll(isInterface, cellI)
+ {
+ if (isInterface[cellI] == 1.0)
+ {
+ targetLevel[cellI] = max(targetLevel[cellI], (refineLevel - i));
+ }
+ }
+ }
+ }
+ }
+
+ // The set refinement physical regions (force the mesh to stay refined
+ // near key features)
+ forAll(refinedRegions_, regionI)
+ {
+ const entry& region = refinedRegions_[regionI];
+
+ autoPtr<topoSetSource> source =
+ topoSetSource::New(region.keyword(), *this, region.dict());
+
+ cellSet selectedCellSet
+ (
+ *this,
+ "cellSet",
+ nCells()/10+1 //estimate
+ );
+
+ source->applyToSet
+ (
+ topoSetSource::NEW,
+ selectedCellSet
+ );
+
+ const labelList cells = selectedCellSet.toc();
+
+ label minLevel = readLabel(region.dict().lookup("minLevel"));
+
+ forAll(cells, i)
+ {
+ const label& cellI = cells[i];
+
+ targetLevel[cellI] = max(targetLevel[cellI], minLevel);
+ }
+ }
+
+ //-DD: add buffer layer based on targetLevel field to prevent 2-to-1 refinement
+
+ volScalarField blockedLevel = targetFld * 0.;
+
+ for (label currLayer=refinePoints[3]; currLayer>=1; currLayer--)
+ {
+ forAll (targetLevel, cellI)
+ {
+ if (targetLevel[cellI] >= currLayer)
+ {
+ blockedLevel[cellI] = targetLevel[cellI];
+ }
+ }
+
+ for (label i=0; i<nBufferLayers_; i++)
+ {
+ blockedLevel = max(blockedLevel, pos(fvc::average(fvc::interpolate(blockedLevel)) - SMALL)*(currLayer-1));
+ }
+
+ labelList blockLev(blockedLevel.internalField().size(), 0);
+ forAll (blockLev, cellI)
+ {
+ blockLev[cellI] = blockedLevel[cellI];
+ }
+ targetLevel = max(targetLevel, blockLev);
+ }
+
+ // mark cells to be refined/unrefined based on above criteria:
+ // simply check, if targetLevel lower or higher cellLevel
+ forAll (intRefFld.internalField(), cellI)
+ {
+ intRefFld.primitiveFieldRef()[cellI] = targetLevel[cellI] - cellLevel[cellI];
+ targetFld.primitiveFieldRef()[cellI] = targetLevel[cellI];
+ currFld.primitiveFieldRef()[cellI] = cellLevel[cellI];
+ }
+
+ intRefFld.correctBoundaryConditions();
+
+ Info<<"Min,max refinement field = " << Foam::min(intRefFld).value() << ", "
+ << Foam::max(intRefFld).value() << endl;
+}
+
+void Foam::dynamicRefineBalancedFvMesh::readRefinementDict()
+{
+ dictionary dynamicMeshDict
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "dynamicMeshDict",
+ time().constant(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE,
+ false
+ )
+ )
+ );
+
+ const dictionary refineDict
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "dynamicMeshDict",
+ time().constant(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE,
+ false
+ )
+ ).subDict("dynamicRefineFvMeshCoeffs")
+ );
+
+ {
+ wordList surfFlds;
+ if (refineDict.readIfPresent("mapSurfaceFields", surfFlds))
+ {
+ // Rework into hashtable.
+ mapSurfaceFields_.resize(surfFlds.size());
+ forAll(surfFlds, i)
+ {
+ mapSurfaceFields_.insert(surfFlds[i], surfFlds[i]);
+ }
+ }
+ }
+
+ nBufferLayers_ = readLabel(refineDict.lookup("nBufferLayers"));
+
+ if( dynamicMeshDict.isDict("refinementControls") )
+ {
+ dictionary refineControlDict =
+ dynamicMeshDict.subDict("refinementControls");
+
+ enableRefinementControl_ =
+ Switch(refineControlDict.lookup("enableRefinementControl"));
+
+ if( enableRefinementControl_ )
+ {
+ // Overwrite field name entry in dynamicRefineFvMeshCoeffs?
+ // For now you just have to be smart and enter
+ // 'internalRefinementField' for the field name manually
+
+ // Read HashTable of field-refinement scalars
+ if( refineControlDict.found("fields") )
+ {
+ fields_ = HashTable< List<scalar> >
+ (
+ refineControlDict.lookup("fields")
+ );
+ }
+
+ // Read HashTable of gradient-refinement scalars
+ if( refineControlDict.found("gradients") )
+ {
+ gradFields_ = HashTable< List<scalar> >
+ (
+ refineControlDict.lookup("gradients")
+ );
+ }
+
+ // Read HashTable of curl-refinement vectors
+ if( refineControlDict.found("curls") )
+ {
+ curlFields_ = HashTable< List<scalar> >
+ (
+ refineControlDict.lookup("curls")
+ );
+ }
+
+ // Read interface refinement data
+ if( refineControlDict.found("interface") )
+ {
+ interface_ = HashTable< dictionary >
+ (
+ refineControlDict.lookup("interface")
+ );
+ }
+
+ // Read refinement regions
+ if( refineControlDict.found("regions") )
+ {
+ refinedRegions_ = PtrList<entry>
+ (
+ refineControlDict.lookup("regions")
+ );
+ }
+ }
+ }
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::dynamicRefineBalancedFvMesh::dynamicRefineBalancedFvMesh
+(
+ const IOobject& io
+)
+:
+ dynamicRefineFvMesh(io),
+ internalRefinementFieldPtr_(NULL),
+ targetLevelPtr_(NULL),
+ isLevelPtr_(NULL),
+ fields_(),
+ gradFields_(),
+ curlFields_(),
+ interface_(),
+ refinedRegions_(),
+ enableRefinementControl_(false),
+ nBufferLayers_(0),
+ rebalance_(false)
+{
+ readRefinementDict();
+
+ if( enableRefinementControl_ )
+ {
+ internalRefinementFieldPtr_ = new volScalarField
+ (
+ IOobject
+ (
+ "internalRefinementField",
+ time().timeName(),
+ *this,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ *this,
+ dimensionedScalar("zero", dimless, 0.0)
+ );
+
+ targetLevelPtr_ = new volScalarField
+ (
+ IOobject
+ (
+ "targetLevel",
+ time().timeName(),
+ *this,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ *this,
+ dimensionedScalar("zero", dimless, 0.0)
+ );
+
+ isLevelPtr_ = new volScalarField
+ (
+ IOobject
+ (
+ "isLevel",
+ time().timeName(),
+ *this,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ *this,
+ dimensionedScalar("zero", dimless, 0.0)
+ );
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::dynamicRefineBalancedFvMesh::~dynamicRefineBalancedFvMesh()
+{
+
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::dynamicRefineBalancedFvMesh::mapFields(const mapPolyMesh& mpm)
+{
+DebugVar(mpm.nOldCells());
+ dynamicRefineFvMesh::mapFields(mpm);
+
+ // Correct surface fields on introduced internal faces. These get
+ // created out-of-nothing so get an interpolated value.
+ mapNewInternalFaces<scalar>(mpm.faceMap());
+ mapNewInternalFaces<vector>(mpm.faceMap());
+ mapNewInternalFaces<sphericalTensor>(mpm.faceMap());
+ mapNewInternalFaces<symmTensor>(mpm.faceMap());
+ mapNewInternalFaces<tensor>(mpm.faceMap());
+}
+
+
+bool Foam::dynamicRefineBalancedFvMesh::update()
+{
+ //Part 0 - Update internally calculated refinement field
+ readRefinementDict();
+
+ List<scalar> refinePoints = readRefinementPoints();
+ label refineInterval = refinePoints[2];
+
+ if( enableRefinementControl_ && time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0 )
+ {
+ updateRefinementField();
+ }
+
+ //Part 1 - Call normal update from dynamicRefineFvMesh
+ bool hasChanged = dynamicRefineFvMesh::update();
+
+ if( Pstream::parRun() && hasChanged)
+ {
+ //Correct values on all coupled patches
+ correctBoundaries<scalar>();
+ correctBoundaries<vector>();
+ correctBoundaries<sphericalTensor>();
+ correctBoundaries<symmTensor>();
+ correctBoundaries<tensor>();
+ }
+
+ dictionary decomposeParDict
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "decomposeParDict",
+ time().system(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE,
+ false
+ )
+ )
+ );
+
+ rebalance_ = false;
+
+ // Part 2 - Load Balancing
+ {
+ dictionary refineDict
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "dynamicMeshDict",
+ time().constant(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE,
+ false
+ )
+ ).subDict("dynamicRefineFvMeshCoeffs")
+ );
+
+ Switch enableBalancing = refineDict.lookup("enableBalancing");
+
+ if ( Pstream::parRun() && hasChanged )
+ {
+ const scalar allowableImbalance =
+ readScalar(refineDict.lookup("allowableImbalance"));
+
+ //First determine current level of imbalance - do this for all
+ // parallel runs with a changing mesh, even if balancing is disabled
+ label nGlobalCells = globalData().nTotalCells();
+ scalar idealNCells = scalar(nGlobalCells)/scalar(Pstream::nProcs());
+ scalar localImbalance = mag(scalar(nCells()) - idealNCells);
+ Foam::reduce(localImbalance, maxOp<scalar>());
+ scalar maxImbalance = localImbalance/idealNCells;
+
+ Info<<"Maximum imbalance = " << 100*maxImbalance << " %" << endl;
+
+ //If imbalanced, construct weighted coarse graph (level 0) with node
+ // weights equal to their number of subcells. This partitioning works
+ // as long as the number of level 0 cells is several times greater than
+ // the number of processors.
+ if( maxImbalance > allowableImbalance && enableBalancing)
+ {
+ Info << "\n**Solver hold for redistribution at time = " << time().timeName() << " s" << endl;
+
+ rebalance_ = true;
+
+ const labelIOList& cellLevel = meshCutter().cellLevel();
+ Map<label> coarseIDmap(100);
+ labelList uniqueIndex(nCells(),0);
+
+ label nCoarse = 0;
+
+ forAll(cells(), cellI)
+ {
+ if( cellLevel[cellI] > 0 )
+ {
+ uniqueIndex[cellI] = nCells() + topParentID
+ (
+ meshCutter().history().parentIndex(cellI)
+ );
+ }
+ else
+ {
+ uniqueIndex[cellI] = cellI;
+ }
+
+ if( coarseIDmap.insert(uniqueIndex[cellI], nCoarse) )
+ {
+ ++nCoarse;
+ }
+ }
+
+ // Convert to local sequential indexing and calculate coarse
+ // points and weights
+ labelList localIndex(nCells(),0);
+ pointField coarsePoints(nCoarse,vector::zero);
+ scalarField coarseWeights(nCoarse,0.0);
+
+ forAll(uniqueIndex, cellI)
+ {
+ localIndex[cellI] = coarseIDmap[uniqueIndex[cellI]];
+
+ // If 2D refinement (quadtree) is ever implemented, this '3'
+ // should be set in general as the number of refinement
+ // dimensions.
+ label w = (1 << (3*cellLevel[cellI]));
+
+ coarseWeights[localIndex[cellI]] += 1.0;
+ coarsePoints[localIndex[cellI]] += C()[cellI]/w;
+ }
+
+ // Set up decomposer - a separate dictionary is used here so
+ // you can use a simple partitioning for decomposePar and
+ // ptscotch for the rebalancing (or any chosen algorithms)
+ autoPtr<decompositionMethod> decomposer
+ (
+ decompositionMethod::New
+ (
+ IOdictionary
+ (
+ IOobject
+ (
+ "balanceParDict",
+ time().system(),
+ *this,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+ )
+ )
+ );
+
+ labelList finalDecomp = decomposer().decompose
+ (
+ *this,
+ localIndex,
+ coarsePoints,
+ coarseWeights
+ );
+
+ scalar tolDim = globalMeshData::matchTol_ * bounds().mag();
+
+ Info<< "Distributing the mesh ..." << endl;
+ fvMeshDistribute distributor(*this, tolDim);
+
+ Info<< "Mapping the fields ..." << endl;
+ autoPtr<mapDistributePolyMesh> map =
+ distributor.distribute(finalDecomp);
+
+ Info<< "Distribute the map ..." << endl;
+ meshCutter_.distribute(map);
+
+
+ Info << "Successfully distributed mesh" << endl;
+
+ scalarList procLoadNew (Pstream::nProcs(), 0.0);
+ procLoadNew[Pstream::myProcNo()] = this->nCells();
+
+ reduce(procLoadNew, sumOp<List<scalar> >());
+
+ scalar overallLoadNew = sum(procLoadNew);
+ scalar averageLoadNew = overallLoadNew/double(Pstream::nProcs());
+
+ Info << "New distribution: " << procLoadNew << endl;
+ Info << "Max deviation: " << max(Foam::mag(procLoadNew-averageLoadNew)/averageLoadNew)*100.0 << " %" << endl;
+ }
+ }
+ }
+
+ if( Pstream::parRun() && rebalance_)
+ {
+ //Correct values on all coupled patches
+ correctBoundaries<scalar>();
+ correctBoundaries<vector>();
+ correctBoundaries<sphericalTensor>();
+ correctBoundaries<symmTensor>();
+ correctBoundaries<tensor>();
+ }
+
+ return hasChanged;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
new file mode 100644
index 0000000000000000000000000000000000000000..f4944d5c850a8c0df73c88e02eeb1654bf4d7aa9
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 Tyler Voskuilen
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is a derivative work of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+
+Class
+ Foam::dynamicRefineBalancedFvMesh
+
+SourceFiles
+ dynamicRefineBalancedFvMesh.C
+ dynamicRefineBalancedFvMeshTemplates.C
+
+Authors
+ T.G. Voskuilen ( https://github.com/tgvoskuilen/meshBalancing )
+ Daniel Deising <deising@mma.tu-darmstadt.de>
+ Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+ All rights reserved.
+
+Description
+ A fvMesh with built-in multi-criterion refinement
+ and run-time load balancing.
+
+ You may refer to this software as :
+ //- full bibliographic data to be provided
+
+ This code has been developed by :
+ Daniel Rettenmaier (main developer).
+
+ Method Development and Intellectual Property :
+ T.G. Voskuilen (Purdue University)
+ Timothée Pourpoint <timothee@purdue.edu> (Purdue University
+ Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+ Daniel Deising <deising@mma.tu-darmstadt.de>
+ Holger Marschall <marschall@csi.tu-darmstadt.de>
+ Dieter Bothe <bothe@csi.tu-darmstadt.de>
+ Cameron Tropea <ctropea@sla.tu-darmstadt.de>
+
+ School of Aeronautics and Astronautics
+ Purdue University
+ Mathematical Modeling and Analysis
+ Institute for Fluid Mechanics and Aerodynamics
+ Center of Smart Interfaces
+ Technische Universitaet Darmstadt
+
+ If you use this software for your scientific work or your publications,
+ please don't forget to acknowledge explicitly the use of it.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicRefineBalancedFvMesh_H
+#define dynamicRefineBalancedFvMesh_H
+
+
+#include "dynamicFvMesh.H"
+#include "dynamicRefineFvMesh.H"
+#include "hexRef8.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
+#include "decompositionMethod.H"
+#include "fvMeshDistribute.H"
+#include "mapDistributePolyMesh.H"
+#include "HashTable.H"
+#include "topoSetSource.H"
+#include "cellSet.H"
+#include "PtrDictionary.H"
+#include "dictionaryEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class dynamicRefineBalancedFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicRefineBalancedFvMesh
+:
+ public dynamicRefineFvMesh
+{
+
+private:
+
+ //- (non-flux)SurfaceFields to map
+ HashTable<word> mapSurfaceFields_;
+
+ //-
+ volScalarField* internalRefinementFieldPtr_;
+ volScalarField* targetLevelPtr_;
+ volScalarField* isLevelPtr_;
+
+ //-
+ HashTable< List<scalar> > fields_;
+
+ //-
+ HashTable< List<scalar> > gradFields_;
+
+ //-
+ HashTable< List<scalar> > curlFields_;
+
+ //-
+ HashTable< dictionary > interface_;
+
+ //-
+ PtrList<entry> refinedRegions_;
+
+ //-
+ Switch enableRefinementControl_;
+
+ //-
+ void updateRefinementField();
+
+ //-
+ void readRefinementDict();
+
+ //-
+ List<scalar> readRefinementPoints();
+
+ //- Disallow default bitwise copy construct
+ dynamicRefineBalancedFvMesh(const dynamicRefineBalancedFvMesh&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const dynamicRefineBalancedFvMesh&);
+
+ //-
+ label topParentID(label p);
+
+ //-
+ label nBufferLayers_;
+
+ //-
+ bool rebalance_;
+
+ //- Map non-flux surface<Type>Fields for new internal faces
+ // (from cell splitting)
+ template <class T>
+ void mapNewInternalFaces(const labelList& faceMap);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("dynamicRefineBalancedFvMesh");
+
+
+ // Constructors
+
+ //- Construct from IOobject
+ explicit dynamicRefineBalancedFvMesh(const IOobject& io);
+
+
+ //- Destructor
+ virtual ~dynamicRefineBalancedFvMesh();
+
+
+ // Member Functions
+
+ //- Update the mesh for both mesh motion and topology change
+ virtual bool update();
+
+ //- Map all fields in time using given map
+ virtual void mapFields(const mapPolyMesh& mpm);
+
+ //- Template to update all volField boundaries
+ template<class Type> void correctBoundaries();
+
+ //-
+ bool rebalance()
+ {
+ return rebalance_;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "dynamicRefineBalancedFvMeshTemplates.C"
+#endif
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..200e3f831d3356fb5316e8785ae55715f6c4add3
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
@@ -0,0 +1,217 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2014 Tyler Voskuilen
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volMesh.H"
+#include "fvPatchField.H"
+#include "surfaceFields.H"
+
+template <class T>
+void Foam::dynamicRefineBalancedFvMesh::mapNewInternalFaces
+(
+ const labelList& faceMap
+)
+{
+ typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+ HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+ const labelUList& owner = this->faceOwner();
+ const labelUList& neighbour = this->faceNeighbour();
+ const dimensionedScalar deltaN = 1e-8 / pow(average(this->V()), 1.0 / 3.0);
+
+ forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+ {
+ GeoField& sFld = *iter();
+ if (mapSurfaceFields_.found(iter.key()))
+ {
+ if (debug)
+ {
+ InfoInFunction << iter.key()<< endl;
+ }
+
+
+ // Create flat version of field
+ Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
+ {
+ forAll(sFld.internalField(), iFace)
+ {
+ tsFld[iFace] = sFld.internalField()[iFace];
+ }
+ forAll(sFld.boundaryField(), iPatch)
+ {
+ const fvsPatchField<T>& pfld = sFld.boundaryField()[iPatch];
+
+ label start = pfld.patch().start();
+ forAll(pfld, faceI)
+ {
+ tsFld[faceI+start] = pfld[faceI];
+ }
+ }
+ }
+
+
+ // Loop over all faces
+ for (label facei = 0; facei < nInternalFaces(); facei++)
+ {
+ label oldFacei = faceMap[facei];
+
+ //- map surface field on newly generated faces
+ if (oldFacei == -1)
+ {
+ //- find owner and neighbour cell
+ const cell& ownFaces = this->cells()[owner[facei]];
+ const cell& neiFaces = this->cells()[neighbour[facei]];
+
+ //- loop over all owner/neighbour cell faces
+ //- and find already mapped ones (master-faces):
+ T tmpValue = pTraits<T>::zero;
+ scalar magFld = 0;
+ label counter = 0;
+
+ //- simple averaging of all neighbour master-faces
+ forAll(ownFaces, iFace)
+ {
+ if (faceMap[ownFaces[iFace]] != -1)
+ {
+ tmpValue += tsFld[ownFaces[iFace]];
+ magFld += mag(tsFld[ownFaces[iFace]]);
+ counter++;
+ }
+ }
+
+ forAll(neiFaces, iFace)
+ {
+ if (faceMap[neiFaces[iFace]] != -1)
+ {
+ tmpValue = tmpValue + tsFld[neiFaces[iFace]];
+ magFld += mag(tsFld[neiFaces[iFace]]);
+ counter++;
+ }
+ }
+
+ if(counter > 0)
+ {
+ if
+ (
+ GeometricField<T, fvsPatchField, surfaceMesh>::typeName
+ == "surfaceScalarField"
+ )
+ {
+ tmpValue /= counter;
+ }
+ else if
+ (
+ GeometricField<T, fvsPatchField, surfaceMesh>::typeName
+ == "surfaceVectorField"
+ )
+ {
+ magFld /= counter;
+ tmpValue *= magFld/(mag(tmpValue)+deltaN.value());
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "mapping implementation only valid for"
+ << " scalar and vector fields! \n Field "
+ << sFld.name() << " is of type: "
+ << GeometricField<T, fvsPatchField, surfaceMesh>::typeName
+ << abort(FatalError);
+ }
+ }
+
+ sFld[facei] = tmpValue;
+ }
+ }
+ }
+ }
+}
+
+
+template<class Type>
+void Foam::dynamicRefineBalancedFvMesh::correctBoundaries()
+{
+ typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+ HashTable<GeoField*> flds(this->objectRegistry::lookupClass<GeoField>());
+
+ forAllIter(typename HashTable<GeoField*>, flds, iter)
+ {
+ GeoField& fld = *iter();
+
+ //mimic "evaluate" but only for coupled patches (processor or cyclic)
+ // and only for blocking or nonBlocking comms (no scheduled comms)
+ if
+ (
+ Pstream::defaultCommsType == Pstream::commsTypes::blocking
+ || Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+ )
+ {
+ label nReq = Pstream::nRequests();
+
+ forAll(fld.boundaryField(), patchi)
+ {
+ if(fld.boundaryField()[patchi].coupled())
+ {
+ fld.boundaryFieldRef()[patchi].initEvaluate
+ (
+ Pstream::defaultCommsType
+ );
+ }
+ }
+
+ // Block for any outstanding requests
+ if
+ (
+ Pstream::parRun()
+ && Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+ )
+ {
+ Pstream::waitRequests(nReq);
+ }
+
+ forAll(fld.boundaryField(), patchi)
+ {
+ if(fld.boundaryField()[patchi].coupled())
+ {
+ fld.boundaryFieldRef()[patchi].evaluate
+ (
+ Pstream::defaultCommsType
+ );
+ }
+ }
+ }
+ else
+ {
+ //Scheduled patch updates not supported
+ FatalErrorIn
+ (
+ "dynamicRefineBalancedFvMeshTemplates::correctBoundaries"
+ ) << "Unsuported communications type "
+ << Pstream::commsTypeNames[Pstream::defaultCommsType]
+ << exit(FatalError);
+ }
+ }
+}
+
+// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 466b7d4e5ef2a904e4ce8b5ea35623f0f002e495..91aad27d1eea8d8f94e4ce4d92038d50c6d1f1d5 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -46,7 +46,7 @@ namespace Foam
Foam::label Foam::dynamicRefineFvMesh::count
(
- const bitSet& l,
+ const PackedBoolList& l,
const unsigned int val
)
{
@@ -58,9 +58,9 @@ Foam::label Foam::dynamicRefineFvMesh::count
n++;
}
- // debug also serves to get-around Clang compiler trying to optimise
+ // debug also serves to get-around Clang compiler trying to optimsie
// out this forAll loop under O3 optimisation
- if (debug)
+ if (debug & 128)
{
Info<< "n=" << n << endl;
}
@@ -72,7 +72,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
void Foam::dynamicRefineFvMesh::calculateProtectedCells
(
- bitSet& unrefineableCell
+ PackedBoolList& unrefineableCell
) const
{
if (protectedCell_.empty())
@@ -86,7 +86,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
unrefineableCell = protectedCell_;
// Get neighbouring cell level
- labelList neiLevel(nBoundaryFaces());
+ labelList neiLevel(nFaces()-nInternalFaces());
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
@@ -103,10 +103,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
forAll(faceNeighbour(), facei)
{
label own = faceOwner()[facei];
+ bool ownProtected = unrefineableCell.get(own);
label nei = faceNeighbour()[facei];
-
- bool ownProtected = unrefineableCell.test(own);
- bool neiProtected = unrefineableCell.test(nei);
+ bool neiProtected = unrefineableCell.get(nei);
if (ownProtected && (cellLevel[nei] > cellLevel[own]))
{
@@ -120,7 +119,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
label own = faceOwner()[facei];
- bool ownProtected = unrefineableCell.test(own);
+ bool ownProtected = unrefineableCell.get(own);
if
(
ownProtected
@@ -142,14 +141,16 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
if (seedFace[facei])
{
label own = faceOwner()[facei];
- if (unrefineableCell.set(own))
+ if (unrefineableCell.get(own) == 0)
{
+ unrefineableCell.set(own, 1);
hasExtended = true;
}
label nei = faceNeighbour()[facei];
- if (unrefineableCell.set(nei))
+ if (unrefineableCell.get(nei) == 0)
{
+ unrefineableCell.set(nei, 1);
hasExtended = true;
}
}
@@ -159,8 +160,9 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
if (seedFace[facei])
{
label own = faceOwner()[facei];
- if (unrefineableCell.set(own))
+ if (unrefineableCell.get(own) == 0)
{
+ unrefineableCell.set(own, 1);
hasExtended = true;
}
}
@@ -205,63 +207,27 @@ void Foam::dynamicRefineFvMesh::readDict()
}
-// Refines cells, maps fields and recalculates (an approximate) flux
-Foam::autoPtr<Foam::mapPolyMesh>
-Foam::dynamicRefineFvMesh::refine
-(
- const labelList& cellsToRefine
-)
+void Foam::dynamicRefineFvMesh::mapFields(const mapPolyMesh& mpm)
{
- // Mesh changing engine.
- polyTopoChange meshMod(*this);
-
- // Play refinement commands into mesh changer.
- meshCutter_.setRefinement(cellsToRefine, meshMod);
-
- // Create mesh (with inflation), return map from old to new mesh.
- autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
- mapPolyMesh& map = *mapPtr;
-
- Info<< "Refined from "
- << returnReduce(map.nOldCells(), sumOp<label>())
- << " to " << globalData().nTotalCells() << " cells." << endl;
-
- if (debug)
- {
- // Check map.
- for (label facei = 0; facei < nInternalFaces(); facei++)
- {
- label oldFacei = map.faceMap()[facei];
-
- if (oldFacei >= nInternalFaces())
- {
- FatalErrorInFunction
- << "New internal face:" << facei
- << " fc:" << faceCentres()[facei]
- << " originates from boundary oldFace:" << oldFacei
- << abort(FatalError);
- }
- }
- }
-
- // // Remove the stored tet base points
- // tetBasePtIsPtr_.clear();
- // // Remove the cell tree
- // cellTreePtr_.clear();
-
- // Update fields
- updateMesh(map);
-
+ dynamicFvMesh::mapFields(mpm);
// Correct the flux for modified/added faces. All the faces which only
// have been renumbered will already have been handled by the mapping.
{
- const labelList& faceMap = map.faceMap();
- const labelList& reverseFaceMap = map.reverseFaceMap();
+ const labelList& faceMap = mpm.faceMap();
+ const labelList& reverseFaceMap = mpm.reverseFaceMap();
// Storage for any master faces. These will be the original faces
// on the coarse cell that get split into four (or rather the
// master face gets modified and three faces get added from the master)
- labelHashSet masterFaces(4*cellsToRefine.size());
+ // Estimate number of faces created
+ labelHashSet masterFaces
+ (
+ max
+ (
+ mag(nFaces()-mpm.nOldFaces())/4,
+ nFaces()/100
+ )
+ );
forAll(faceMap, facei)
{
@@ -393,8 +359,10 @@ Foam::dynamicRefineFvMesh::refine
}
// Update master faces
- for (const label facei : masterFaces)
+ forAllConstIter(labelHashSet, masterFaces, iter)
{
+ label facei = iter.key();
+
if (isInternalFace(facei))
{
phi[facei] = phiU[facei];
@@ -415,6 +383,86 @@ Foam::dynamicRefineFvMesh::refine
}
}
+ // Correct the flux for injected faces - these are the faces which have
+ // no correspondence to the old mesh (i.e. added without a masterFace, edge
+ // or point). An example is the internal faces from hexRef8.
+ {
+ const labelList& faceMap = mpm.faceMap();
+
+ mapNewInternalFaces<scalar>(this->Sf(), this->magSf(), faceMap);
+ mapNewInternalFaces<vector>(this->Sf(), this->magSf(), faceMap);
+
+ // No oriented fields of more complex type
+ mapNewInternalFaces<sphericalTensor>(faceMap);
+ mapNewInternalFaces<symmTensor>(faceMap);
+ mapNewInternalFaces<tensor>(faceMap);
+ }
+}
+
+
+// Refines cells, maps fields and recalculates (an approximate) flux
+Foam::autoPtr<Foam::mapPolyMesh>
+Foam::dynamicRefineFvMesh::refine
+(
+ const labelList& cellsToRefine
+)
+{
+ // Mesh changing engine.
+ polyTopoChange meshMod(*this);
+
+ // Play refinement commands into mesh changer.
+ meshCutter_.setRefinement(cellsToRefine, meshMod);
+
+ // Create mesh (with inflation), return map from old to new mesh.
+ //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+ autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
+
+ Info<< "Refined from "
+ << returnReduce(map().nOldCells(), sumOp<label>())
+ << " to " << globalData().nTotalCells() << " cells." << endl;
+
+ if (debug)
+ {
+ // Check map.
+ for (label facei = 0; facei < nInternalFaces(); facei++)
+ {
+ label oldFacei = map().faceMap()[facei];
+
+ if (oldFacei >= nInternalFaces())
+ {
+ FatalErrorInFunction
+ << "New internal face:" << facei
+ << " fc:" << faceCentres()[facei]
+ << " originates from boundary oldFace:" << oldFacei
+ << abort(FatalError);
+ }
+ }
+ }
+
+ // // Remove the stored tet base points
+ // tetBasePtIsPtr_.clear();
+ // // Remove the cell tree
+ // cellTreePtr_.clear();
+
+ // Update fields
+ updateMesh(map);
+
+
+ // Move mesh
+ /*
+ pointField newPoints;
+ if (map().hasMotionPoints())
+ {
+ newPoints = map().preMotionPoints();
+ }
+ else
+ {
+ newPoints = points();
+ }
+ movePoints(newPoints);
+ */
+
+
// Update numbering of cells/vertices.
meshCutter_.updateMesh(map);
@@ -422,12 +470,12 @@ Foam::dynamicRefineFvMesh::refine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
- bitSet newProtectedCell(nCells());
+ PackedBoolList newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
- const label oldCelli = map.cellMap()[celli];
- newProtectedCell.set(celli, protectedCell_.test(oldCelli));
+ label oldCelli = map().cellMap()[celli];
+ newProtectedCell.set(celli, protectedCell_.get(oldCelli));
}
protectedCell_.transfer(newProtectedCell);
}
@@ -435,7 +483,7 @@ Foam::dynamicRefineFvMesh::refine
// Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0));
- return mapPtr;
+ return map;
}
@@ -482,21 +530,36 @@ Foam::dynamicRefineFvMesh::unrefine
// Change mesh and generate map.
- autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(*this, false);
- mapPolyMesh& map = *mapPtr;
+ //autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, true);
+ autoPtr<mapPolyMesh> map = meshMod.changeMesh(*this, false);
Info<< "Unrefined from "
- << returnReduce(map.nOldCells(), sumOp<label>())
+ << returnReduce(map().nOldCells(), sumOp<label>())
<< " to " << globalData().nTotalCells() << " cells."
<< endl;
// Update fields
updateMesh(map);
+
+ // Move mesh
+ /*
+ pointField newPoints;
+ if (map().hasMotionPoints())
+ {
+ newPoints = map().preMotionPoints();
+ }
+ else
+ {
+ newPoints = points();
+ }
+ movePoints(newPoints);
+ */
+
// Correct the flux for modified faces.
{
- const labelList& reversePointMap = map.reversePointMap();
- const labelList& reverseFaceMap = map.reverseFaceMap();
+ const labelList& reversePointMap = map().reversePointMap();
+ const labelList& reverseFaceMap = map().reverseFaceMap();
HashTable<surfaceScalarField*> fluxes
(
@@ -585,14 +648,14 @@ Foam::dynamicRefineFvMesh::unrefine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
- bitSet newProtectedCell(nCells());
+ PackedBoolList newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
- label oldCelli = map.cellMap()[celli];
+ label oldCelli = map().cellMap()[celli];
if (oldCelli >= 0)
{
- newProtectedCell.set(celli, protectedCell_.test(oldCelli));
+ newProtectedCell.set(celli, protectedCell_.get(oldCelli));
}
}
protectedCell_.transfer(newProtectedCell);
@@ -601,7 +664,7 @@ Foam::dynamicRefineFvMesh::unrefine
// Debug: Check refinement levels (across faces only)
meshCutter_.checkRefinementLevels(-1, labelList(0));
- return mapPtr;
+ return map;
}
@@ -688,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
const scalar lowerRefineLevel,
const scalar upperRefineLevel,
const scalarField& vFld,
- bitSet& candidateCell
+ PackedBoolList& candidateCell
) const
{
// Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
@@ -711,7 +774,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
{
if (cellError[celli] > 0)
{
- candidateCell.set(celli);
+ candidateCell.set(celli, 1);
}
}
}
@@ -721,7 +784,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
(
const label maxCells,
const label maxRefinement,
- const bitSet& candidateCell
+ const PackedBoolList& candidateCell
) const
{
// Every refined cell causes 7 extra cells
@@ -731,7 +794,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
// Mark cells that cannot be refined since they would trigger refinement
// of protected cells (since 2:1 cascade)
- bitSet unrefineableCell;
+ PackedBoolList unrefineableCell;
calculateProtectedCells(unrefineableCell);
// Count current selection
@@ -748,8 +811,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
if
(
cellLevel[celli] < maxRefinement
- && candidateCell.test(celli)
- && !unrefineableCell.test(celli)
+ && candidateCell.get(celli)
+ && (
+ unrefineableCell.empty()
+ || !unrefineableCell.get(celli)
+ )
)
{
candidates.append(celli);
@@ -766,8 +832,11 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
if
(
cellLevel[celli] == level
- && candidateCell.test(celli)
- && !unrefineableCell.test(celli)
+ && candidateCell.get(celli)
+ && (
+ unrefineableCell.empty()
+ || !unrefineableCell.get(celli)
+ )
)
{
candidates.append(celli);
@@ -802,7 +871,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
(
const scalar unrefineLevel,
- const bitSet& markedCell,
+ const PackedBoolList& markedCell,
const scalarField& pFld
) const
{
@@ -815,7 +884,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
// If we have any protected cells make sure they also are not being
// unrefined
- bitSet protectedPoint(nPoints());
+ PackedBoolList protectedPoint(nPoints());
if (protectedCell_.size())
{
@@ -828,9 +897,9 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
{
label cellI = pCells[pCellI];
- if (protectedCell_.test(cellI))
+ if (protectedCell_[cellI])
{
- protectedPoint.set(pointI);
+ protectedPoint[pointI] = true;
break;
}
}
@@ -870,7 +939,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
forAll(pCells, pCelli)
{
- if (markedCell.test(pCells[pCelli]))
+ if (markedCell.get(pCells[pCelli]))
{
hasMarked = true;
break;
@@ -907,7 +976,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
void Foam::dynamicRefineFvMesh::extendMarkedCells
(
- bitSet& markedCell
+ PackedBoolList& markedCell
) const
{
// Mark faces using any marked cell
@@ -915,7 +984,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
forAll(markedCell, celli)
{
- if (markedCell.test(celli))
+ if (markedCell.get(celli))
{
const cell& cFaces = cells()[celli];
@@ -933,15 +1002,15 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
{
if (markedFace[facei])
{
- markedCell.set(faceOwner()[facei]);
- markedCell.set(faceNeighbour()[facei]);
+ markedCell.set(faceOwner()[facei], 1);
+ markedCell.set(faceNeighbour()[facei], 1);
}
}
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
if (markedFace[facei])
{
- markedCell.set(faceOwner()[facei]);
+ markedCell.set(faceOwner()[facei], 1);
}
}
}
@@ -949,7 +1018,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
(
- bitSet& protectedCell,
+ PackedBoolList& protectedCell,
label& nProtected
) const
{
@@ -971,7 +1040,7 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
// Check if cell has already 8 anchor points -> protect cell
if (nAnchorPoints[celli] == 8)
{
- if (protectedCell.set(celli))
+ if (protectedCell.set(celli, true))
{
nProtected++;
}
@@ -988,9 +1057,9 @@ void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
forAll(protectedCell, celli)
{
- if (!protectedCell.test(celli) && nAnchorPoints[celli] != 8)
+ if (!protectedCell[celli] && nAnchorPoints[celli] != 8)
{
- protectedCell.set(celli);
+ protectedCell.set(celli, true);
nProtected++;
}
}
@@ -1005,7 +1074,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
meshCutter_(*this),
dumpLevel_(false),
nRefinementIterations_(0),
- protectedCell_(nCells(), false)
+ protectedCell_(nCells(), 0)
{
// Read static part of dictionary
readDict();
@@ -1034,7 +1103,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
label celli = pCells[i];
- if (!protectedCell_.test(celli))
+ if (!protectedCell_.get(celli))
{
if (pointLevel[pointi] <= cellLevel[celli])
{
@@ -1042,7 +1111,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
if (nAnchors[celli] > 8)
{
- protectedCell_.set(celli);
+ protectedCell_.set(celli, 1);
nProtected++;
}
}
@@ -1105,9 +1174,9 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (protectedFace[facei])
{
- protectedCell_.set(faceOwner()[facei]);
+ protectedCell_.set(faceOwner()[facei], 1);
nProtected++;
- protectedCell_.set(faceNeighbour()[facei]);
+ protectedCell_.set(faceNeighbour()[facei], 1);
nProtected++;
}
}
@@ -1115,7 +1184,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (protectedFace[facei])
{
- protectedCell_.set(faceOwner()[facei]);
+ protectedCell_.set(faceOwner()[facei], 1);
nProtected++;
}
}
@@ -1127,7 +1196,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
if (cFaces.size() < 6)
{
- if (protectedCell_.set(celli))
+ if (protectedCell_.set(celli, 1))
{
nProtected++;
}
@@ -1138,7 +1207,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
{
if (faces()[cFaces[cFacei]].size() < 4)
{
- if (protectedCell_.set(celli))
+ if (protectedCell_.set(celli, 1))
{
nProtected++;
}
@@ -1158,6 +1227,7 @@ Foam::dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
}
else
{
+
cellSet protectedCells(*this, "protectedCells", nProtected);
forAll(protectedCell_, celli)
{
@@ -1188,9 +1258,9 @@ Foam::dynamicRefineFvMesh::~dynamicRefineFvMesh()
bool Foam::dynamicRefineFvMesh::update()
{
- // Re-read dictionary. Usually small so takes trivial amount of time
- // compared to actual refinement. Also very useful to be able to modify
- // on-the-fly.
+ // Re-read dictionary. Choosen since usually -small so trivial amount
+ // of time compared to actual refinement. Also very useful to be able
+ // to modify on-the-fly.
dictionary refineDict
(
IOdictionary
@@ -1264,18 +1334,16 @@ bool Foam::dynamicRefineFvMesh::update()
refineDict.get<scalar>("lowerRefineLevel");
const scalar upperRefineLevel =
refineDict.get<scalar>("upperRefineLevel");
-
const scalar unrefineLevel = refineDict.lookupOrDefault<scalar>
(
"unrefineLevel",
GREAT
);
-
const label nBufferLayers =
refineDict.get<label>("nBufferLayers");
// Cells marked for refinement or otherwise protected from unrefinement.
- bitSet refineCell(nCells());
+ PackedBoolList refineCell(nCells());
// Determine candidates for refinement (looking at field only)
selectRefineCandidates
@@ -1316,7 +1384,7 @@ bool Foam::dynamicRefineFvMesh::update()
const labelList& cellMap = map().cellMap();
const labelList& reverseCellMap = map().reverseCellMap();
- bitSet newRefineCell(cellMap.size());
+ PackedBoolList newRefineCell(cellMap.size());
forAll(cellMap, celli)
{
@@ -1324,15 +1392,15 @@ bool Foam::dynamicRefineFvMesh::update()
if (oldCelli < 0)
{
- newRefineCell.set(celli);
+ newRefineCell.set(celli, 1);
}
else if (reverseCellMap[oldCelli] != celli)
{
- newRefineCell.set(celli);
+ newRefineCell.set(celli, 1);
}
else
{
- newRefineCell.set(celli, refineCell.test(oldCelli));
+ newRefineCell.set(celli, refineCell.get(oldCelli));
}
}
refineCell.transfer(newRefineCell);
@@ -1380,7 +1448,7 @@ bool Foam::dynamicRefineFvMesh::update()
if ((nRefinementIterations_ % 10) == 0)
{
- // Compact refinement history occasionally (how often?).
+ // Compact refinement history occassionally (how often?).
// Unrefinement causes holes in the refinementHistory.
const_cast<refinementHistory&>(meshCutter().history()).compact();
}
@@ -1403,8 +1471,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
(
IOstream::streamFormat fmt,
IOstream::versionNumber ver,
- IOstream::compressionType cmp,
- const bool valid
+ IOstream::compressionType cmp
) const
{
// Force refinement data to go to the current time directory.
@@ -1412,8 +1479,8 @@ bool Foam::dynamicRefineFvMesh::writeObject
bool writeOk =
(
- dynamicFvMesh::writeObject(fmt, ver, cmp, valid)
- && meshCutter_.write(valid)
+ dynamicFvMesh::writeObject(fmt, ver, cmp)
+ && meshCutter_.write()
);
if (dumpLevel_)
@@ -1430,7 +1497,7 @@ bool Foam::dynamicRefineFvMesh::writeObject
false
),
*this,
- dimensionedScalar(dimless, Zero)
+ dimensionedScalar("level", dimless, 0)
);
const labelList& cellLevel = meshCutter_.cellLevel();
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
index babbceabea16533c071745a7768ccd93f4100865..2ccd3244866eb69ed6af56835f049857b6004650 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -56,6 +56,7 @@ Description
(rho*phi none) //none)
(ghf none) //NaN) //none)
);
+
// Write the refinement level as a volScalarField
dumpLevel true;
@@ -117,10 +118,10 @@ protected:
//- Refine cells. Update mesh and fields.
- autoPtr<mapPolyMesh> refine(const labelList&);
+ virtual autoPtr<mapPolyMesh> refine(const labelList&);
//- Unrefine cells. Gets passed in centre points of cells to combine.
- autoPtr<mapPolyMesh> unrefine(const labelList&);
+ virtual autoPtr<mapPolyMesh> unrefine(const labelList&);
// Selection of cells to un/refine
@@ -185,6 +186,32 @@ protected:
label& nProtected
) const;
+ //- Map single non-flux surface<Type>Field
+ // for new internal faces (e.g. AMR refine). This currently
+ // interpolates values from surrounding faces (faces on
+ // neighbouring cells) that do have a value.
+ template <class T>
+ void mapNewInternalFaces
+ (
+ const labelList& faceMap,
+ GeometricField<T, fvsPatchField, surfaceMesh>&
+ );
+
+ //- Correct surface fields for new faces
+ template <class T>
+ void mapNewInternalFaces(const labelList& faceMap);
+
+ //- Correct surface fields for new faces. Converts any oriented
+ // fields into non-oriented (i.e. phi -> Uf) before mapping
+ template <class T>
+ void mapNewInternalFaces
+ (
+ const surfaceVectorField& Sf,
+ const surfaceScalarField& magSf,
+ const labelList& faceMap
+ );
+
+
private:
//- No copy construct
@@ -232,6 +259,9 @@ public:
//- Update the mesh for both mesh motion and topology change
virtual bool update();
+ //- Map all fields in time using given map.
+ virtual void mapFields(const mapPolyMesh& mpm);
+
// Writing
@@ -243,16 +273,20 @@ public:
IOstream::compressionType cmp,
const bool valid
) const;
-
};
-
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+ #include "dynamicRefineFvMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
#endif
// ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..6bf655a7ec36c7dca235a428bf164c47a2b3beb6
--- /dev/null
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMeshTemplates.C
@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+ const labelList& faceMap,
+ GeometricField<T, fvsPatchField, surfaceMesh>& sFld
+)
+{
+ typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+
+ //- Make flat field for ease of looping
+ Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
+ SubField<T>(tsFld, this->nInternalFaces()) = sFld.internalField();
+
+ const typename GeoField::Boundary& bFld = sFld.boundaryField();
+ forAll(bFld, patchi)
+ {
+ label facei = this->boundaryMesh()[patchi].start();
+ for (const T& val : bFld[patchi])
+ {
+ tsFld[facei++] = val;
+ }
+ }
+
+ const labelUList& owner = this->faceOwner();
+ const labelUList& neighbour = this->faceNeighbour();
+ const cellList& cells = this->cells();
+
+ for (label facei = 0; facei < nInternalFaces(); facei++)
+ {
+ label oldFacei = faceMap[facei];
+
+ // Map surface field on newly generated faces by obtaining the
+ // hull of the outside faces
+ if (oldFacei == -1)
+ {
+ // Loop over all owner/neighbour cell faces
+ // and find already mapped ones (master-faces):
+ T tmpValue = pTraits<T>::zero;
+ label counter = 0;
+
+ const cell& ownFaces = cells[owner[facei]];
+ for (auto ownFacei : ownFaces)
+ {
+ if (faceMap[ownFacei] != -1)
+ {
+ tmpValue += tsFld[ownFacei];
+ counter++;
+ }
+ }
+
+ const cell& neiFaces = cells[neighbour[facei]];
+ for (auto neiFacei : neiFaces)
+ {
+ if (faceMap[neiFacei] != -1)
+ {
+ tmpValue += tsFld[neiFacei];
+ counter++;
+ }
+ }
+
+ if (counter > 0)
+ {
+ sFld[facei] = tmpValue/counter;
+ }
+ }
+ }
+}
+
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+ const labelList& faceMap
+)
+{
+ typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+ HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+ forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+ {
+ //if (mapSurfaceFields_.found(iter.key()))
+ {
+ if (debug)
+ {
+ Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
+ << " Mapping new internal faces by interpolation on "
+ << iter.key()<< endl;
+ }
+
+ GeoField& sFld = *iter();
+
+ if (sFld.oriented()())
+ {
+ WarningInFunction << "Ignoring mapping oriented field "
+ << sFld.name() << " since of type " << sFld.type()
+ << endl;
+ }
+ else
+ {
+ mapNewInternalFaces(faceMap, sFld);
+ }
+ }
+ }
+}
+
+template<class T>
+void Foam::dynamicRefineFvMesh::mapNewInternalFaces
+(
+ const surfaceVectorField& Sf,
+ const surfaceScalarField& magSf,
+ const labelList& faceMap
+)
+{
+ typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+ HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+ forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+ {
+ //if (mapSurfaceFields_.found(iter.key()))
+ {
+ if (debug)
+ {
+ Info<< "dynamicRefineFvMesh::mapNewInternalFaces():"
+ << " Mapping new internal faces by interpolation on "
+ << iter.key()<< endl;
+ }
+
+ GeoField& sFld = *iter();
+
+ if (sFld.oriented()())
+ {
+ if (debug)
+ {
+ Info<< "dynamicRefineFvMesh::mapNewInternalFaces(): "
+ << "Converting oriented field " << iter.key()
+ << " to intensive field and mapping" << endl;
+ }
+
+ // Assume any oriented field is face area weighted (i.e. a flux)
+ // Convert to intensive (& oriented) before mapping. Untested.
+
+ typedef GeometricField
+ <
+ typename outerProduct<vector, T>::type,
+ fvsPatchField,
+ surfaceMesh
+ > NormalGeoField;
+
+ // Convert to intensive and non oriented
+ NormalGeoField fFld(sFld*Sf/Foam::sqr(magSf));
+
+ // Interpolate
+ mapNewInternalFaces(faceMap, fFld);
+
+ // Convert back to extensive and oriented
+ sFld = (fFld & Sf);
+ }
+ else
+ {
+ mapNewInternalFaces(faceMap, sFld);
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index a66bd7c48636b3114e5c095d5173c7b87682b264..f535a71cd4a438ee53555e343d41190b11e9dec8 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -210,7 +210,7 @@ Foam::label Foam::hexRef8::addInternalFace
nei, // neighbour
-1, // master point
-1, // master edge
- meshFacei, // master face for addition
+ -1, // master face for addition
false, // flux flip
-1, // patch for face
-1, // zone for face
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index a0a032935b2884b08bd2e027389fae842084c908..41dc37f78422fc70228995cc609e23d2ae824268 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -571,6 +571,13 @@ const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
{
if (!lduPtr_)
{
+ if (debug)
+ {
+ InfoInFunction
+ << " calculating fvMeshLduAddressing from nFaces:"
+ << nFaces() << endl;
+ }
+
lduPtr_ = new fvMeshLduAddressing(*this);
}
@@ -821,6 +828,9 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
// Update polyMesh. This needs to keep volume existent!
polyMesh::updateMesh(mpm);
+ // Our slice of the addressing is no longer valid
+ deleteDemandDrivenData(lduPtr_);
+
if (VPtr_)
{
// Grab old time volumes if the time has been incremented
diff --git a/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict b/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
index 88613b395f8a76e75bc238b10ca6751a78d7c52f..88a8b82bb86c09881cfec7f563c7f927f306de48 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
+++ b/tutorials/multiphase/interFoam/laminar/damBreakWithObstacle/constant/dynamicMeshDict
@@ -49,6 +49,7 @@ correctFluxes
(rhoPhi none)
(alphaPhi0.water none)
(ghf none)
+ (alphaPhiUn none)
);
// Write the refinement level as a volScalarField