diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
index 5f1fc7451c2e41d21061dba1cdeb96f2cee151d0..ba9cbf752c7374aa9f25608c89df044f6fd1e90a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H
@@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
-//mesh.interpolate(HbyA);
+volVectorField HbyA("HbyA", U);
+HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
if (pimple.nCorrPISO() <= 1)
{
diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H
index 5e9861ccc383ec482ff1e04ff06dbc461e87ebc4..68b208f41b1d30938ac5bbcf2e571c0082dc474c 100644
--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H
@@ -16,6 +16,9 @@
fvOptions.constrain(UEqn);
- solve(UEqn == -cellMask*fvc::grad(p));
+ if (simple.momentumPredictor())
+ {
+ solve(UEqn == -cellMask*fvc::grad(p));
+ }
fvOptions.correct(U);
diff --git a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
index c64f92e33f798c6bd762797bfd6a9118e468d838..78555ec5ce8436c7a7e64e9765518295c63421a5 100644
--- a/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H
@@ -5,7 +5,7 @@
surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
- HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
+ HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
//mesh.interpolate(HbyA);
if (massFluxInterpolation)
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
index 266b339b67570e4b2db6134d9daec29dfa0a6059..c6bde4234f6f3452665a8b943deafe791d30d7b0 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H
@@ -138,3 +138,4 @@ volScalarField pDivU
mesh,
dimensionedScalar("pDivU", p.dimensions()/dimTime, 0)
);
+
diff --git a/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H b/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
index d8481c88704f42298d1b7f2c2fad8f843ce8c135..a2acf920692b33ac049e3a166a3c7c6981e5bacf 100644
--- a/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H
@@ -1,6 +1,6 @@
#include "readTimeControls.H"
-correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", false);
checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
diff --git a/applications/test/rigidBodyDynamics/pendulum/pendulum b/applications/test/rigidBodyDynamics/pendulum/pendulum
index d179af4438547ee888802d8074f94b04aeedbad3..3a6b4b18433f54fa2e9de49221406620864a5290 100644
--- a/applications/test/rigidBodyDynamics/pendulum/pendulum
+++ b/applications/test/rigidBodyDynamics/pendulum/pendulum
@@ -32,6 +32,7 @@ bodies
type sphere;
mass 1;
radius 0.05;
+ centreOfMass (0 10 0);
transform (1 0 0 0 1 0 0 0 1) (0 -1 0);
mergeWith hinge;
}
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring b/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring
index ec9c3517b7dd78e5e2f715591fb9e531003e74a9..4f4fdc44c1c8fdd47ee424df748f4f36b3bc1ac3 100644
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring
+++ b/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring
@@ -14,6 +14,7 @@ bodies
parent root;
mass 6e4;
radius 0.01;
+ centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
@@ -37,6 +38,7 @@ bodies
parent M;
mass 6e3;
radius 0.01;
+ centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 -5 0);
mergeWith M;
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index b0bfc8d58571d94fedfc4af53ca2ef102f13e2d3..b2783c5563e76c41e82a3f3a26ae742c77f9625c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -233,10 +233,10 @@ public:
void reorder(const labelUList&, const bool validBoundary);
//- writeData member function required by regIOobject
- bool writeData(Ostream&) const;
+ virtual bool writeData(Ostream&) const;
//- Write using given format, version and form uncompression
- bool writeObject
+ virtual bool writeObject
(
IOstream::streamFormat fmt,
IOstream::versionNumber ver,
diff --git a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
index 7dd34b9343e3f0d5f5f1183115df65b2c441c7ba..542dc624cf7bdff18b3057d2c30b8424788f9428 100644
--- a/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
+++ b/src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C
@@ -101,7 +101,18 @@ bool Foam::patchDistMethods::Poisson::correct
volVectorField gradyPsi(fvc::grad(yPsi));
volScalarField magGradyPsi(mag(gradyPsi));
- y = sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi;
+ // Need to stabilise the y for overset meshes since the holed cells
+ // keep the initial value (0.0) so the gradient of that will be
+ // zero as well. Turbulence models do not like zero wall distance.
+ y = max
+ (
+ sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi,
+ dimensionedScalar("smallY", dimLength, SMALL)
+ );
+
+ // For overset: enforce smooth y field (yPsi is smooth, magGradyPsi is
+ // not)
+ mesh_.interpolate(y);
// Need to stabilise the y for overset meshes since the holed cells
// keep the initial value (0.0) so the gradient of that will be
diff --git a/src/overset/Make/files b/src/overset/Make/files
index a285c4979cc925fc794b52b95129686c2bf1ccfa..14ce43d1f27894ad3e70dcc5ffe470a9d9806aaa 100644
--- a/src/overset/Make/files
+++ b/src/overset/Make/files
@@ -6,6 +6,7 @@ cellCellStencil/inverseDistance/waveMethod.C
cellCellStencil/inverseDistance/meshToMeshData.C
cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
+cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
dynamicOversetFvMesh/dynamicOversetFvMesh.C
@@ -21,6 +22,9 @@ oversetPolyPatch/oversetGAMGInterfaceField.C
oversetPolyPatch/oversetPointPatch.C
oversetPolyPatch/oversetPointPatchFields.C
+oversetPolyPatch/semiImplicitOversetFvPatchFields.C
+oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C
+
oversetAdjustPhi/oversetAdjustPhi.C
regionsToCell/regionsToCell.C
diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C
index 9a2851870d5b682fd916a6d179efb2986f2f1d6b..15b2c00a6ba335533bc319d28dd567e38a194d96 100644
--- a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C
+++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C
@@ -94,22 +94,22 @@ Foam::cellCellStencil::~cellCellStencil()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
+const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
{
- if (!mesh_.foundObject<labelIOList>("zoneID"))
+ if (!mesh.foundObject<labelIOList>("zoneID"))
{
labelIOList* zoneIDPtr = new labelIOList
(
IOobject
(
"zoneID",
- mesh_.facesInstance(),
+ mesh.facesInstance(),
polyMesh::meshSubDir,
- mesh_,
+ mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
- mesh_.nCells()
+ mesh.nCells()
);
labelIOList& zoneID = *zoneIDPtr;
@@ -118,13 +118,13 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
IOobject
(
"zoneID",
- mesh_.time().findInstance(mesh_.dbDir(), "zoneID"),
- mesh_,
+ mesh.time().findInstance(mesh.dbDir(), "zoneID"),
+ mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
- mesh_
+ mesh
);
forAll(volZoneID, cellI)
{
@@ -133,7 +133,7 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
zoneIDPtr->store();
}
- return mesh_.lookupObject<labelIOList>("zoneID");
+ return mesh.lookupObject<labelIOList>("zoneID");
}
diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H
index b46f5339932f504f7e94e431a75ea92cdec9110a..05a3d7cf590163ade3b6e00e012222808733222d 100644
--- a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H
+++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H
@@ -168,12 +168,27 @@ public:
// calculated, 1 = use interpolated)
virtual const scalarList& cellInterpolationWeight() const = 0;
+ //- Calculate weights for a single acceptor
+ virtual void stencilWeights
+ (
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+ ) const = 0;
+
//- Helper: is stencil fully local
bool localStencil(const labelUList&) const;
//- Helper: get reference to registered zoneID. Loads volScalarField
// if not registered.
- const labelIOList& zoneID() const;
+ static const labelIOList& zoneID(const fvMesh&);
+
+ //- Helper: get reference to registered zoneID. Loads volScalarField
+ // if not registered.
+ const labelIOList& zoneID() const
+ {
+ return zoneID(mesh_);
+ }
//- Helper: create cell-cell addressing in global numbering
static void globalCellCells
diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
index 2a0e3ca15df14c9560d36d9f46ab1025b4b965f7..c465b7c101dd5972a9e5277f52ce55487f0b12f9 100644
--- a/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
+++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H
@@ -154,6 +154,17 @@ public:
{
return stencilPtr_().cellInterpolationWeight();
}
+
+ //- Calculate weights for a single acceptor
+ virtual void stencilWeights
+ (
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+ ) const
+ {
+ stencilPtr_().stencilWeights(sample, donorCcs, weights);
+ }
};
diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
index 85c20494d7038fc188e459cee1d419eb943c0dd1..d030bc1623da30d8f09b38f64a8a2adb6f23c1c7 100644
--- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
+++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
@@ -1100,4 +1100,39 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
}
+void Foam::cellCellStencils::cellVolumeWeight::stencilWeights
+(
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+) const
+{
+ // Inverse-distance weighting
+
+ weights.setSize(donorCcs.size());
+ scalar sum = 0.0;
+ forAll(donorCcs, i)
+ {
+ scalar d = mag(sample-donorCcs[i]);
+
+ if (d > ROOTVSMALL)
+ {
+ weights[i] = 1.0/d;
+ sum += weights[i];
+ }
+ else
+ {
+ // Short circuit
+ weights = 0.0;
+ weights[i] = 1.0;
+ return;
+ }
+ }
+ forAll(weights, i)
+ {
+ weights[i] /= sum;
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H
index 4331b1421f5e4c0b82b6c10fac66492e15c319ef..6c1a8003be768d320d44c8a484d095a6e8d954af 100644
--- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H
+++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H
@@ -227,6 +227,15 @@ public:
{
return cellInterpolationWeight_;
}
+
+ //- Calculate inverse distance weights for a single acceptor. Revert
+ // to inverse distance (so not consistent with volume overlap!)
+ virtual void stencilWeights
+ (
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+ ) const;
};
diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
index eb0d12915947bb00da7d9cb95b247efa0b871bcc..ef5a35724a3f3c9f0a8edd4838224541666a71d0 100644
--- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C
@@ -155,6 +155,7 @@ void Foam::cellCellStencils::inverseDistance::fill
void Foam::cellCellStencils::inverseDistance::markBoundaries
(
const fvMesh& mesh,
+ const vector& smallVec,
const boundBox& bb,
const labelVector& nDivs,
@@ -188,6 +189,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
+ faceBb.min() -= smallVec;
+ faceBb.max() += smallVec;
+
if (bb.overlaps(faceBb))
{
fill(patchTypes, bb, nDivs, faceBb, patchCellType::PATCH);
@@ -213,6 +217,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
+ faceBb.min() -= smallVec;
+ faceBb.max() += smallVec;
+
if (bb.overlaps(faceBb))
{
fill(patchTypes, bb, nDivs, faceBb, patchCellType::OVERSET);
@@ -333,6 +340,10 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
+ treeBoundBox cBb(cellBb(mesh_, celli));
+ cBb.min() -= smallVec_;
+ cBb.max() += smallVec_;
+
if
(
overlaps
@@ -340,7 +351,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
srcPatchBb,
zoneDivs,
srcPatchTypes,
- cellBb(mesh_, celli),
+ cBb,
patchCellType::PATCH
)
)
@@ -399,6 +410,9 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
+ treeBoundBox cBb(cellBb(mesh_, celli));
+ cBb.min() -= smallVec_;
+ cBb.max() += smallVec_;
if
(
overlaps
@@ -406,7 +420,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
srcPatchBb,
zoneDivs,
srcPatchTypes,
- cellBb(mesh_, celli),
+ cBb,
patchCellType::PATCH
)
)
@@ -501,7 +515,9 @@ void Foam::cellCellStencils::inverseDistance::markDonors
label celli = tgtCellMap[tgtCelli];
if (allStencil[celli].empty())
{
- const treeBoundBox subBb(cellBb(mesh_, celli));
+ treeBoundBox subBb(cellBb(mesh_, celli));
+ subBb.min() -= smallVec_;
+ subBb.max() += smallVec_;
forAll(srcOverlapProcs, i)
{
@@ -1194,9 +1210,17 @@ void Foam::cellCellStencils::inverseDistance::walkFront
{
if (isA<oversetFvPatch>(fvm[patchI]))
{
- forAll(fvm[patchI], i)
+ const labelList& fc = fvm[patchI].faceCells();
+ forAll(fc, i)
{
- isFront[fvm[patchI].start()+i] = true;
+ label cellI = fc[i];
+ if (allCellTypes[cellI] == INTERPOLATED)
+ {
+ // Note that acceptors might have been marked hole if
+ // there are no donors in which case we do not want to
+ // walk this out. This is an extreme situation.
+ isFront[fvm[patchI].start()+i] = true;
+ }
}
}
}
@@ -1353,12 +1377,12 @@ void Foam::cellCellStencils::inverseDistance::walkFront
}
-void Foam::cellCellStencils::inverseDistance::calcStencilWeights
+void Foam::cellCellStencils::inverseDistance::stencilWeights
(
const point& sample,
const pointList& donorCcs,
scalarList& weights
-)
+) const
{
// Inverse-distance weighting
@@ -1500,7 +1524,7 @@ void Foam::cellCellStencils::inverseDistance::createStencil
{
label cellI = donorCells[i];
const pointList& donorCentres = donorCellCentres[cellI];
- calcStencilWeights
+ stencilWeights
(
samples[cellI],
donorCentres,
@@ -1565,6 +1589,7 @@ Foam::cellCellStencils::inverseDistance::inverseDistance
:
cellCellStencil(mesh),
dict_(dict),
+ smallVec_(vector::zero),
cellTypes_(labelList(mesh.nCells(), CALCULATED)),
interpolationCells_(0),
cellInterpolationMap_(),
@@ -1605,6 +1630,9 @@ bool Foam::cellCellStencils::inverseDistance::update()
{
scalar layerRelax(dict_.lookupOrDefault("layerRelax", 1.0));
+ scalar tol = dict_.lookupOrDefault("tolerance", 1e-10);
+ smallVec_ = mesh_.bounds().span()*tol;
+
const labelIOList& zoneID = this->zoneID();
label nZones = gMax(zoneID)+1;
@@ -1745,6 +1773,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
markBoundaries
(
meshParts[zoneI].subMesh(),
+ smallVec_,
patchBb[zoneI][Pstream::myProcNo()],
patchDivisions[zoneI],
@@ -1758,7 +1787,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Print a bit
{
- Info<< typeName << " : detected " << nZones
+ Info<< type() << " : detected " << nZones
<< " mesh regions" << endl;
Info<< incrIndent;
forAll(nCellsPerZone, zoneI)
@@ -1862,6 +1891,16 @@ bool Foam::cellCellStencils::inverseDistance::update()
}
else
{
+ // If there are no donors we can either set the
+ // acceptor cell to 'hole' i.e. nothing gets calculated
+ // if the acceptor cells go outside the donor meshes or
+ // we can set it to 'calculated' to have something
+ // like an acmi type behaviour where only covered
+ // acceptor cell use the interpolation and non-covered
+ // become calculated. Uncomment below line. Note: this
+ // should be switchable?
+ //allCellTypes[cellI] = CALCULATED;
+
allCellTypes[cellI] = HOLE;
}
}
@@ -1915,7 +1954,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
- Pout<< typeName << " : dumping injectionStencil to "
+ Pout<< type() << " : dumping injectionStencil to "
<< str.name() << endl;
pointField cc(mesh_.cellCentres());
cellInterpolationMap_.distribute(cc);
@@ -1973,7 +2012,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
// Dump stencil
mkDir(mesh_.time().timePath());
OBJstream str(mesh_.time().timePath()/"stencil.obj");
- Pout<< typeName << " : dumping to " << str.name() << endl;
+ Pout<< type() << " : dumping to " << str.name() << endl;
pointField cc(mesh_.cellCentres());
cellInterpolationMap_.distribute(cc);
diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H
index ffc9a291ef4b00cd6633f855bc57b38c62cb8b76..1a57c6034048986a1a4d0dc5be6a4d2c17cbacda 100644
--- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H
+++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H
@@ -74,6 +74,9 @@ protected:
//- Dictionary of motion control parameters
const dictionary dict_;
+ //- Small perturbation vector for geometric tests
+ vector smallVec_;
+
//- Per cell the cell type
labelList cellTypes_;
@@ -143,6 +146,8 @@ protected:
static void markBoundaries
(
const fvMesh& mesh,
+ const vector& smallVec,
+
const boundBox& bb,
const labelVector& nDivs,
PackedList<2>& patchTypes,
@@ -237,15 +242,6 @@ protected:
scalarField& allWeight
) const;
-
- //- Calculate inverse distance weights
- static void calcStencilWeights
- (
- const point& sample,
- const pointList& donorCcs,
- scalarList& weights
- );
-
//- Create stencil starting from the donor containing the acceptor
virtual void createStencil(const globalIndex&);
@@ -319,6 +315,14 @@ public:
{
return cellInterpolationWeight_;
}
+
+ //- Calculate inverse distance weights for a single acceptor
+ virtual void stencilWeights
+ (
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+ ) const;
};
diff --git a/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.C b/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
new file mode 100644
index 0000000000000000000000000000000000000000..8f2da7a0ba010e66b977be872121847797d8e5cd
--- /dev/null
+++ b/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify i
+ 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 "leastSquaresCellCellStencil.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SVD.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cellCellStencils
+{
+ defineTypeNameAndDebug(leastSquares, 0);
+ addToRunTimeSelectionTable(cellCellStencil, leastSquares, mesh);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::cellCellStencils::leastSquares::stencilWeights
+(
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+) const
+{
+ // Implicit least squares weighting
+ // Number of donors
+ label nD = donorCcs.size();
+
+ weights.setSize(nD);
+
+ // List for distance vectors and LSQ weights
+ List<vector> d(nD);
+ scalarList LSQw(nD);
+
+ // Sum of weights
+ scalar W = 0;
+
+ // Sum of weighted distance vectors
+ vector dw = vector::zero;
+
+ RectangularMatrix<scalar> A(nD, 3);
+
+ bool shortC = false;
+
+ // Compute distance vectors and fill rectangular matrix
+ forAll(donorCcs, j)
+ {
+ // Neighbour weights
+ d[j] = donorCcs[j] - sample;
+
+ // Check for short-circuiting if zero distance
+ // is detected with respect to any donor
+ if (mag(d[j]) < ROOTVSMALL)
+ {
+ shortC = true;
+ break;
+ }
+
+ LSQw[j] = 1.0/magSqr(d[j]);
+
+ // T matrix
+ vector wd = LSQw[j]*d[j];
+ A[j][0] = wd.x();
+ A[j][1] = wd.y();
+ A[j][2] = wd.z();
+
+ // Sum of weighted distance vectors
+ dw += wd;
+
+ // Sum of weights
+ W += LSQw[j];
+ }
+
+ if (!shortC)
+ {
+ // Use Singular Value Decomposition to avoid problems
+ // with 1D, 2D stencils
+ SVD svd(A.T()*A, SMALL);
+
+ // Least squares vectors
+ RectangularMatrix<scalar> ATAinvAT(svd.VSinvUt()*A.T());
+
+ scalar saveDiag(W);
+
+ // Diagonal coefficient
+ for (label i = 0; i < 3; i++)
+ {
+ // Get row
+ scalarList Row(UList<scalar>(ATAinvAT[i], nD));
+
+ forAll(donorCcs, j)
+ {
+ W -= Row[j]*LSQw[j]*dw[i];
+ }
+ }
+
+ if (mag(W) < SMALL*mag(saveDiag))
+ {
+ shortC = true;
+ }
+ else
+ {
+ // Compute final neighbour weights with additional scaling
+ forAll(donorCcs, j)
+ {
+ weights[j] =
+ (
+ LSQw[j]
+ - ATAinvAT[0][j]*LSQw[j]*dw[0]
+ - ATAinvAT[1][j]*LSQw[j]*dw[1]
+ - ATAinvAT[2][j]*LSQw[j]*dw[2]
+ )
+ /W;
+ }
+ }
+ }
+
+ if (shortC)
+ {
+ // Matrix ill conditioned. Use straight injection from central
+ // donor.
+ weights = 0.0;
+ weights[0] = 1.0;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::cellCellStencils::leastSquares::leastSquares
+(
+ const fvMesh& mesh,
+ const dictionary& dict,
+ const bool doUpdate
+)
+:
+ inverseDistance(mesh, dict, false)
+{
+ if (doUpdate)
+ {
+ update();
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::cellCellStencils::leastSquares::~leastSquares()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.H b/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.H
new file mode 100644
index 0000000000000000000000000000000000000000..4d9b88eb0fe31347c5dc844fe31d208b950f79d3
--- /dev/null
+++ b/src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.H
@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+Class
+ Foam::cellCellStencils::leastSquares
+
+Description
+ Least-squares-weighted interpolation stencil.
+
+ Base machinery is similar to inverse distance interpolation stencil
+ but weights minimize error in LSQ sense recovering exact solution
+ for linear solution problems. Gradient and values are found
+ simultaneously.
+
+SourceFiles
+ leastSquaresCellCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellCellStencils_leastSquares_H
+#define cellCellStencils_leastSquares_H
+
+#include "inverseDistanceCellCellStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace cellCellStencils
+{
+
+/*---------------------------------------------------------------------------*\
+ Class leastSquares Declaration
+\*---------------------------------------------------------------------------*/
+
+class leastSquares
+:
+ public inverseDistance
+{
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ leastSquares(const leastSquares&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const leastSquares&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("leastSquares");
+
+
+ // Constructors
+
+ //- Construct from fvMesh
+ leastSquares(const fvMesh&, const dictionary&, const bool);
+
+
+ //- Destructor
+ virtual ~leastSquares();
+
+
+ // Member Functions
+
+ //- Calculate lsq weights for single acceptor
+ virtual void stencilWeights
+ (
+ const point& sample,
+ const pointList& donorCcs,
+ scalarList& weights
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace cellCellStencils
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
index 7a9b2906f6bb08bcf37a48ef03a287182d6fd767..b2e3efccd50b3fd5129fae3e263fd3775c79c586 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
+++ b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
@@ -51,6 +51,7 @@ namespace cellCellStencils
void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
(
const fvMesh& mesh,
+ const vector& smallVec,
const boundBox& bb,
const labelVector& nDivs,
@@ -84,6 +85,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
+ faceBb.min() -= smallVec;
+ faceBb.max() += smallVec;
+
if (bb.overlaps(faceBb))
{
voxelMeshSearch::fill
@@ -116,6 +120,8 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
// Mark in voxel mesh
boundBox faceBb(pp.points(), pp[i], false);
+ faceBb.min() -= smallVec;
+ faceBb.max() += smallVec;
if (bb.overlaps(faceBb))
{
voxelMeshSearch::fill
@@ -165,6 +171,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
+ boundBox cBb(allPoints, allCellPoints[celli], false);
+ cBb.min() -= smallVec_;
+ cBb.max() += smallVec_;
if
(
@@ -172,7 +181,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
(
srcPatchBb,
srcDivs,
- boundBox(allPoints, allCellPoints[celli], false),
+ cBb,
srcPatchTypes,
static_cast<unsigned int>(patchCellType::PATCH)
)
@@ -231,13 +240,17 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
forAll(tgtCellMap, tgtCelli)
{
label celli = tgtCellMap[tgtCelli];
+ boundBox cBb(allPoints, allCellPoints[celli], false);
+ cBb.min() -= smallVec_;
+ cBb.max() += smallVec_;
+
if
(
voxelMeshSearch::overlaps
(
srcPatchBb,
srcDivs,
- boundBox(allPoints, allCellPoints[celli], false),
+ cBb,
srcPatchTypes,
static_cast<unsigned int>(patchCellType::PATCH)
)
@@ -628,6 +641,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
markBoundaries
(
meshParts_[zonei].subMesh(),
+ smallVec_,
patchBb[zonei][Pstream::myProcNo()],
patchDivisions[zonei],
diff --git a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H
index f0cafa3521c7fa1ad541d4ff2fe54ab50710407b..bcc49000c445110fcedb022929d6af609261abc5 100644
--- a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H
+++ b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H
@@ -77,6 +77,7 @@ protected:
static void markBoundaries
(
const fvMesh& mesh,
+ const vector& smallVec,
const boundBox& bb,
const labelVector& nDivs,
diff --git a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C
index deb346f0816a0fcf838e15ec09ecd91b9693e271..526e813b6c62bc8ba5c8695f2cc3bfa8dbad281a 100644
--- a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C
+++ b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C
@@ -228,7 +228,10 @@ Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io)
:
dynamicMotionSolverFvMesh(io),
active_(false)
-{}
+{
+ // Force loading zoneID field before time gets incremented
+ (void)cellCellStencil::zoneID(*this);
+}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
diff --git a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H
index 0b61ed42a70b2eb9bc51d0fe3e6f6b1e78d3fc96..3f8d9782f4846f6f855e417cf07a1ba5b37ef20e 100644
--- a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H
+++ b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H
@@ -38,6 +38,8 @@ SourceFiles
#include "dynamicMotionSolverFvMesh.H"
#include "labelIOList.H"
#include "fvMeshPrimitiveLduAddressing.H"
+#include "lduInterfaceFieldPtrsList.H"
+#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -100,6 +102,14 @@ protected:
template<class GeoField>
void interpolate(GeoField& psi) const;
+ //- Freeze values at holes
+ template<class Type>
+ void freezeHoles(fvMatrix<Type>&) const;
+
+ //- Get scalar interfaces of certain type
+ template<class GeoField, class PatchType>
+ lduInterfaceFieldPtrsList scalarInterfaces(const GeoField& psi) const;
+
//- Add interpolation to matrix (coefficients)
template<class Type>
void addInterpolation(fvMatrix<Type>&) const;
@@ -112,25 +122,6 @@ protected:
template<class GeoField>
static void correctCoupledBoundaryConditions(GeoField& fld);
- //- Use extended addressing
- void active(const bool f) const
- {
- active_ = f;
-
- if (active_)
- {
- DebugInfo<< "Switching to extended addressing with nFaces:"
- << primitiveLduAddr().lowerAddr().size()
- << endl;
- }
- else
- {
- DebugInfo<< "Switching to base addressing with nFaces:"
- << fvMesh::lduAddr().lowerAddr().size()
- << endl;
- }
- }
-
private:
@@ -164,12 +155,6 @@ public:
// Extended addressing
- //- Use extended addressing
- bool active() const
- {
- return active_;
- }
-
//- Return extended ldu addressing
const fvMeshPrimitiveLduAddressing& primitiveLduAddr() const;
@@ -177,6 +162,37 @@ public:
// primitiveLduAddr
virtual const lduAddressing& lduAddr() const;
+ //- Return old to new face addressing
+ const labelList& reverseFaceMap() const
+ {
+ return reverseFaceMap_;
+ }
+
+ //- Return true if using extended addressing
+ bool active() const
+ {
+ return active_;
+ }
+
+ //- Enable/disable extended addressing
+ void active(const bool f) const
+ {
+ active_ = f;
+
+ if (active_)
+ {
+ DebugInfo<< "Switching to extended addressing with nFaces:"
+ << primitiveLduAddr().lowerAddr().size()
+ << endl;
+ }
+ else
+ {
+ DebugInfo<< "Switching to base addressing with nFaces:"
+ << fvMesh::lduAddr().lowerAddr().size()
+ << endl;
+ }
+ }
+
// Overset
diff --git a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C
index dce8e9cd7a86c4d0848476f47db121147b6d4319..006c576d096686634fa77a71f733ef771d194771 100644
--- a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C
+++ b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C
@@ -26,6 +26,7 @@ License
#include "volFields.H"
#include "fvMatrix.H"
#include "cellCellStencilObject.H"
+#include "oversetFvPatchField.H"
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
@@ -77,41 +78,34 @@ void Foam::dynamicOversetFvMesh::interpolate(GeoField& psi) const
}
-//template<class T>
-//void Foam::dynamicOversetFvMesh::setReference(fvMatrix<T>& m) const
-//{
-// const mapDistribute& map = cellInterpolationMap();
-// const List<scalarList>& wghts = cellInterpolationWeights();
-// const labelListList& stencil = cellStencil();
-// const labelList& cellIDs = interpolationCells();
-// const scalarList& factor = cellInterpolationWeight();
-//
-// Field<T> work(m.psi());
-// map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
-//
-// forAll(cellIDs, i)
-// {
-// label celli = cellIDs[i];
-//
-// const scalarList& w = wghts[celli];
-// const labelList& nbrs = stencil[celli];
-// const scalar f = factor[celli];
-//
-// //Pout<< "Interpolating " << celli << " from values:" << endl;
-// T s(pTraits<T>::zero);
-// forAll(nbrs, nbrI)
-// {
-// //Pout<< " " << work[nbrs[nbrI]]
-// // << " from slot " << nbrs[nbrI] << endl;
-// s += w[nbrI]*work[nbrs[nbrI]];
-// }
-// //Pout<< "Interpolated value:" << s << endl;
-// const T oldPsi = m.psi()[celli];
-// m.setReference(celli, (1.0-f)*oldPsi + f*s, true);
-// Pout<< "psi was:" << oldPsi << " now forcing to:"
-// << (1.0-f)*oldPsi + f*s << endl;
-// }
-//}
+template<class GeoField, class PatchType>
+Foam::lduInterfaceFieldPtrsList
+Foam::dynamicOversetFvMesh::scalarInterfaces(const GeoField& psi) const
+{
+ const typename GeoField::Boundary& bpsi = psi.boundaryField();
+
+ lduInterfaceFieldPtrsList interfaces(bpsi.size());
+
+ forAll(interfaces, patchi)
+ {
+ if
+ (
+ isA<lduInterfaceField>(bpsi[patchi])
+ && isA<PatchType>(bpsi[patchi])
+ )
+ {
+ interfaces.set
+ (
+ patchi,
+ &refCast<const lduInterfaceField>
+ (
+ bpsi[patchi]
+ )
+ );
+ }
+ }
+ return interfaces;
+}
template<class Type>
@@ -185,6 +179,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
}
}
+
forAll(m.internalCoeffs(), patchI)
{
const labelUList& fc = addr.patchAddr(patchI);
@@ -209,6 +204,76 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
}
}
+ // NOTE: For a non-scalar matrix, in the function fvMatrix::solveSegregated
+ // it uses updateInterfaceMatrix to correct for remote source contributions.
+ // This unfortunately also triggers the overset interpolation which then
+ // produces a non-zero source for interpolated cells.
+ // The procedure below calculates this contribution 'correctionSource'
+ // and substracts it from the source later in order to compensate.
+ Field<Type> correctionSource(diag.size(), pTraits<Type>::zero);
+
+ if (pTraits<Type>::nComponents > 1)
+ {
+ typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+ typename pTraits<Type>::labelType validComponents
+ (
+ this->template validComponents<Type>()
+ );
+
+ // Get all the overset lduInterfaces
+ lduInterfaceFieldPtrsList interfaces
+ (
+ scalarInterfaces<GeoField, oversetFvPatchField<Type>>
+ (
+ m.psi()
+ )
+ );
+
+ const Field<Type>& psiInt = m.psi().primitiveField();
+
+ for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+ {
+ if (component(validComponents, cmpt) != -1)
+ {
+ const scalarField psiCmpt(psiInt.component(cmpt));
+ scalarField corrSourceCmpt(correctionSource.component(cmpt));
+
+ forAll(interfaces, patchi)
+ {
+ if (interfaces.set(patchi))
+ {
+ interfaces[patchi].initInterfaceMatrixUpdate
+ (
+ corrSourceCmpt,
+ true,
+ psiCmpt,
+ m.boundaryCoeffs()[patchi].component(cmpt),
+ cmpt,
+ Pstream::defaultCommsType
+ );
+ }
+ }
+ forAll(interfaces, patchi)
+ {
+ if (interfaces.set(patchi))
+ {
+ interfaces[patchi].updateInterfaceMatrix
+ (
+ corrSourceCmpt,
+ true,
+ psiCmpt,
+ m.boundaryCoeffs()[patchi].component(cmpt),
+ cmpt,
+ Pstream::defaultCommsType
+ );
+ }
+ }
+ correctionSource.replace(cmpt, corrSourceCmpt);
+ }
+ }
+ }
+
// Modify matrix
@@ -271,12 +336,14 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
{
// Create interpolation stencil. Leave any non-local contributions
// to be done by oversetFvPatchField updateMatrixInterface
+ // 'correctionSource' is only applied for vector and higher
+ // fvMatrix; it is zero for scalar matrices (see above).
diag[celli] *= (1.0-f);
diag[celli] += normalisation*f;
source[celli] *= (1.0-f);
- source[celli] += normalisation*f*pTraits<Type>::zero; // dummy
+ source[celli] -= correctionSource[celli];
//Pout<< "Interpolating " << celli << " with factor:" << f
// << " new diag:" << diag[celli]
@@ -328,6 +395,39 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
const dictionary& dict
) const
{
+ // Check we're running with bcs that can handle implicit matrix manipulation
+ const typename GeometricField<Type, fvPatchField, volMesh>::Boundary bpsi =
+ m.psi().boundaryField();
+
+ bool hasOverset = false;
+ forAll(bpsi, patchi)
+ {
+ if (isA<oversetFvPatchField<Type>>(bpsi[patchi]))
+ {
+ hasOverset = true;
+ break;
+ }
+ }
+
+ if (!hasOverset)
+ {
+ if (debug)
+ {
+ Pout<< "dynamicOversetFvMesh::solve() :"
+ << " bypassing matrix adjustment for field " << m.psi().name()
+ << endl;
+ }
+ return dynamicMotionSolverFvMesh::solve(m, dict);
+ }
+
+ if (debug)
+ {
+ Pout<< "dynamicOversetFvMesh::solve() :"
+ << " adjusting matrix for interpolation for field "
+ << m.psi().name() << endl;
+ }
+
+
// Switch to extended addressing (requires mesh::update() having been
// called)
active(true);
@@ -346,8 +446,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
// Use lower level solver
SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
- // Reset matrix (shuffle faces into original order)
- //resetMatrix(m);
+ // Restore matrix
m.upper().transfer(oldUpper);
m.lower().transfer(oldLower);
m.internalCoeffs().transfer(oldInt);
diff --git a/src/overset/oversetAdjustPhi/oversetAdjustPhi.C b/src/overset/oversetAdjustPhi/oversetAdjustPhi.C
index b6eb644b5502d8841ea896559eb4a116abff4933..28c533c761e9f0c22f5b4c9a38f0bd1f7f350145 100644
--- a/src/overset/oversetAdjustPhi/oversetAdjustPhi.C
+++ b/src/overset/oversetAdjustPhi/oversetAdjustPhi.C
@@ -92,7 +92,7 @@ bool Foam::oversetAdjustPhi
else
{
adjustableMassOut[zonei] += flux;
- }
+ }
}
}
}
diff --git a/src/overset/oversetPolyPatch/oversetFvPatchField.C b/src/overset/oversetPolyPatch/oversetFvPatchField.C
index de61b8c118e6c335a50c73a4ff3ffd7f690041b0..da288827af0ab0ac8beeacb4f8ad92f0b559a3fc 100644
--- a/src/overset/oversetPolyPatch/oversetFvPatchField.C
+++ b/src/overset/oversetPolyPatch/oversetFvPatchField.C
@@ -24,6 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "volFields.H"
+#include "cellCellStencil.H"
+#include "dynamicOversetFvMesh.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -34,9 +36,7 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
- LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
- zeroGradientFvPatchField<Type>(p, iF),
- oversetPatch_(refCast<const oversetFvPatch>(p))
+ semiImplicitOversetFvPatchField<Type>(p, iF)
{}
@@ -49,21 +49,8 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const fvPatchFieldMapper& mapper
)
:
- LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
- zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
- oversetPatch_(refCast<const oversetFvPatch>(p))
-{
- if (!isA<oversetFvPatch>(this->patch()))
- {
- FatalErrorInFunction
- << " patch type '" << p.type()
- << "' not constraint type '" << typeName << "'"
- << "\n for patch " << p.name()
- << " of field " << this->internalField().name()
- << " in file " << this->internalField().objectPath()
- << exit(FatalIOError);
- }
-}
+ semiImplicitOversetFvPatchField<Type>(ptf, p, iF, mapper)
+{}
template<class Type>
@@ -74,26 +61,8 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const dictionary& dict
)
:
- LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
- zeroGradientFvPatchField<Type>(p, iF, dict),
- oversetPatch_(refCast<const oversetFvPatch>(p))
-{
- if (!isA<oversetFvPatch>(p))
- {
- FatalIOErrorInFunction(dict)
- << " patch type '" << p.type()
- << "' not constraint type '" << typeName << "'"
- << "\n for patch " << p.name()
- << " of field " << this->internalField().name()
- << " in file " << this->internalField().objectPath()
- << exit(FatalIOError);
- }
-
- if (!dict.found("value") && this->coupled())
- {
- this->evaluate(Pstream::commsTypes::blocking);
- }
-}
+ semiImplicitOversetFvPatchField<Type>(p, iF, dict)
+{}
template<class Type>
@@ -102,9 +71,7 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const oversetFvPatchField<Type>& ptf
)
:
- LduInterfaceField<Type>(ptf.oversetPatch_),
- zeroGradientFvPatchField<Type>(ptf),
- oversetPatch_(ptf.oversetPatch_)
+ semiImplicitOversetFvPatchField<Type>(ptf)
{}
@@ -115,22 +82,19 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
const DimensionedField<Type, volMesh>& iF
)
:
- LduInterfaceField<Type>(ptf.oversetPatch_),
- zeroGradientFvPatchField<Type>(ptf, iF),
- oversetPatch_(ptf.oversetPatch_)
+ semiImplicitOversetFvPatchField<Type>(ptf, iF)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
template<class Type>
void Foam::oversetFvPatchField<Type>::initEvaluate
(
const Pstream::commsTypes commsType
)
{
- if (oversetPatch_.master())
+ if (this->oversetPatch_.master())
{
// Trigger interpolation
const fvMesh& mesh = this->internalField().mesh();
@@ -139,24 +103,14 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
if (&mesh.lduAddr() != &mesh.fvMesh::lduAddr())
{
- // Running extended addressing. Interpolate always.
+ // Running extended addressing. Flux correction already done
+ // in the linear solver (through the initUpdateInterfaceMatrix
+ // method below)
if (debug)
{
- Info<< "Interpolating solved-for field " << fldName << endl;
+ Info<< "Skipping overset interpolation for solved-for field "
+ << fldName << endl;
}
-
- // Interpolate without bc update (since would trigger infinite
- // recursion back to oversetFvPatchField<Type>::evaluate)
- // The only problem is bcs that use the new cell values
- // (e.g. zeroGradient, processor). These need to appear -after-
- // the 'overset' bc.
- mesh.interpolate
- (
- const_cast<Field<Type>&>
- (
- this->primitiveField()
- )
- );
}
else if
(
@@ -206,31 +160,6 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
}
-template<class Type>
-void Foam::oversetFvPatchField<Type>::evaluate(const Pstream::commsTypes)
-{
- if (!this->updated())
- {
- this->updateCoeffs();
- }
-
- zeroGradientFvPatchField<Type>::evaluate();
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
-(
- scalarField& result,
- const bool add,
- const scalarField& psiInternal,
- const scalarField& coeffs,
- const direction cmpt,
- const Pstream::commsTypes commsType
-) const
-{}
-
-
template<class Type>
void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
(
@@ -243,7 +172,10 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
) const
{
// Add remote values
- if (oversetPatch_.master())
+
+ const oversetFvPatch& ovp = this->oversetPatch_;
+
+ if (ovp.master())
{
const fvMesh& mesh = this->patch().boundaryMesh().mesh();
@@ -251,11 +183,10 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
// TBD. This should be cleaner.
if (&mesh.lduAddr() == &mesh.fvMesh::lduAddr())
{
- //Pout<< "** not running extended addressing..." << endl;
return;
}
- const labelListList& stencil = oversetPatch_.stencil();
+ const labelListList& stencil = ovp.stencil();
if (stencil.size() != psiInternal.size())
{
@@ -263,12 +194,11 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
<< " stencil:" << stencil.size() << exit(FatalError);
}
- const mapDistribute& map = oversetPatch_.cellInterpolationMap();
- const List<scalarList>& wghts =
- oversetPatch_.cellInterpolationWeights();
- const labelList& cellIDs = oversetPatch_.interpolationCells();
- const scalarList& factor = oversetPatch_.cellInterpolationWeight();
- const scalarField& normalisation = oversetPatch_.normalisation();
+ const mapDistribute& map = ovp.cellInterpolationMap();
+ const List<scalarList>& wghts = ovp.cellInterpolationWeights();
+ const labelList& cellIDs = ovp.interpolationCells();
+ const scalarList& factor = ovp.cellInterpolationWeight();
+ const scalarField& normalisation = ovp.normalisation();
// Since we're inside initEvaluate/evaluate there might be processor
@@ -306,41 +236,4 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
}
-template<class Type>
-void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
-(
- Field<Type>&,
- const bool add,
- const Field<Type>&,
- const scalarField&,
- const Pstream::commsTypes commsType
-) const
-{
- NotImplemented;
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
-(
- Field<Type>& result,
- const bool add,
- const Field<Type>& psiInternal,
- const scalarField&,
- const Pstream::commsTypes
-) const
-{
- NotImplemented;
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::write(Ostream& os) const
-{
- zeroGradientFvPatchField<Type>::write(os);
- // Make sure to write the value for ease of postprocessing.
- this->writeEntry("value", os);
-}
-
-
// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/oversetFvPatchField.H b/src/overset/oversetPolyPatch/oversetFvPatchField.H
index b7c923fcd5f7875f74a0e091eb1a30fbd116c4e8..9cae63dc43caff5a9882fcb9114da016b45dc458 100644
--- a/src/overset/oversetPolyPatch/oversetFvPatchField.H
+++ b/src/overset/oversetPolyPatch/oversetFvPatchField.H
@@ -28,8 +28,9 @@ Group
grpCoupledBoundaryConditions
Description
-
-SeeAlso
+ Boundary condition for use on overset patches. To be run in combination
+ with special dynamicFvMesh type that inserts local contributions into
+ the matrix. This boundary condition adds the remote contributions.
SourceFiles
oversetFvPatchField.C
@@ -39,9 +40,7 @@ SourceFiles
#ifndef oversetFvPatchField_H
#define oversetFvPatchField_H
-#include "oversetFvPatch.H"
-#include "zeroGradientFvPatchField.H"
-#include "LduInterfaceField.H"
+#include "semiImplicitOversetFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -55,18 +54,8 @@ namespace Foam
template<class Type>
class oversetFvPatchField
:
- public LduInterfaceField<Type>,
- public zeroGradientFvPatchField<Type>
+ public semiImplicitOversetFvPatchField<Type>
{
- // Private data
-
- //- Local reference cast into the overset patch
- const oversetFvPatch& oversetPatch_;
-
- //- Master patch ID
- mutable label masterPatchID_;
-
-
public:
//- Runtime type information
@@ -133,34 +122,11 @@ public:
// Member functions
- // Access
-
- //- Return local reference cast into the cyclic AMI patch
- const oversetFvPatch& oversetPatch() const
- {
- return oversetPatch_;
- }
-
-
// Evaluation functions
//- Initialise the evaluation of the patch field
virtual void initEvaluate(const Pstream::commsTypes commsType);
- //- Evaluate the patch field
- virtual void evaluate(const Pstream::commsTypes commsType);
-
- //- Initialise neighbour matrix update
- virtual void initInterfaceMatrixUpdate
- (
- scalarField&,
- const bool add,
- const scalarField&,
- const scalarField&,
- const direction,
- const Pstream::commsTypes commsType
- ) const;
-
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
@@ -171,32 +137,6 @@ public:
const direction cmpt,
const Pstream::commsTypes commsType
) const;
-
- //- Initialise neighbour matrix update
- virtual void initInterfaceMatrixUpdate
- (
- Field<Type>&,
- const bool add,
- const Field<Type>&,
- const scalarField&,
- const Pstream::commsTypes commsType
- ) const;
-
- //- Update result field based on interface functionality
- virtual void updateInterfaceMatrix
- (
- Field<Type>&,
- const bool add,
- const Field<Type>&,
- const scalarField&,
- const Pstream::commsTypes commsType
- ) const;
-
-
- // I-O
-
- //- Write
- virtual void write(Ostream& os) const;
};
diff --git a/src/overset/oversetPolyPatch/oversetGAMGInterfaceField.H b/src/overset/oversetPolyPatch/oversetGAMGInterfaceField.H
index 3cd44fa3dba456cfa1505ef2f15d73b924a19020..e706be9c14ebba645885a1a125abeef519f15a18 100644
--- a/src/overset/oversetPolyPatch/oversetGAMGInterfaceField.H
+++ b/src/overset/oversetPolyPatch/oversetGAMGInterfaceField.H
@@ -52,12 +52,16 @@ class oversetGAMGInterfaceField
:
public GAMGInterfaceField
{
- // Private data
+protected:
+
+ // Protected data
//- Local reference cast into the interface
const oversetGAMGInterface& oversetInterface_;
+private:
+
// Private Member Functions
//- Disallow default bitwise copy construct
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.C b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.C
new file mode 100644
index 0000000000000000000000000000000000000000..cd361f3ad8571c7abaa4a273a50ef610d214b670
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.C
@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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 "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF
+)
+:
+ LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+ zeroGradientFvPatchField<Type>(p, iF),
+ oversetPatch_(refCast<const oversetFvPatch>(p))
+{}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+ const semiImplicitOversetFvPatchField<Type>& ptf,
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF,
+ const fvPatchFieldMapper& mapper
+)
+:
+ LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+ zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
+ oversetPatch_(refCast<const oversetFvPatch>(p))
+{
+ if (!isA<oversetFvPatch>(this->patch()))
+ {
+ FatalErrorInFunction
+ << " patch type '" << p.type()
+ << "' not constraint type '" << typeName << "'"
+ << "\n for patch " << p.name()
+ << " of field " << this->internalField().name()
+ << " in file " << this->internalField().objectPath()
+ << exit(FatalIOError);
+ }
+}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+ const fvPatch& p,
+ const DimensionedField<Type, volMesh>& iF,
+ const dictionary& dict
+)
+:
+ LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+ zeroGradientFvPatchField<Type>(p, iF, dict),
+ oversetPatch_(refCast<const oversetFvPatch>(p))
+{
+ if (!isA<oversetFvPatch>(p))
+ {
+ FatalIOErrorInFunction(dict)
+ << " patch type '" << p.type()
+ << "' not constraint type '" << typeName << "'"
+ << "\n for patch " << p.name()
+ << " of field " << this->internalField().name()
+ << " in file " << this->internalField().objectPath()
+ << exit(FatalIOError);
+ }
+
+ if (!dict.found("value") && this->coupled())
+ {
+ this->evaluate(Pstream::commsTypes::blocking);
+ }
+}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+ const semiImplicitOversetFvPatchField<Type>& ptf
+)
+:
+ LduInterfaceField<Type>(ptf.oversetPatch_),
+ zeroGradientFvPatchField<Type>(ptf),
+ oversetPatch_(ptf.oversetPatch_)
+{}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+ const semiImplicitOversetFvPatchField<Type>& ptf,
+ const DimensionedField<Type, volMesh>& iF
+)
+:
+ LduInterfaceField<Type>(ptf.oversetPatch_),
+ zeroGradientFvPatchField<Type>(ptf, iF),
+ oversetPatch_(ptf.oversetPatch_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::evaluate
+(
+ const Pstream::commsTypes
+)
+{
+ if (debug)
+ {
+ Pout<< FUNCTION_NAME << " field " << this->internalField().name()
+ << " patch " << this->patch().name() << endl;
+ }
+
+
+ if (!this->updated())
+ {
+ this->updateCoeffs();
+ }
+
+ zeroGradientFvPatchField<Type>::evaluate();
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::initEvaluate
+(
+ const Pstream::commsTypes commsType
+)
+{
+ if (debug)
+ {
+ Pout<< FUNCTION_NAME << " field " << this->internalField().name()
+ << " patch " << this->patch().name() << endl;
+ }
+
+ if (this->oversetPatch().master())
+ {
+ // Trigger interpolation
+ const fvMesh& mesh = this->internalField().mesh();
+ //const dictionary& fvSchemes = mesh.schemesDict();
+ const word& fldName = this->internalField().name();
+
+ if (&mesh.lduAddr() != &mesh.fvMesh::lduAddr())
+ {
+ // Running extended addressing so called from within fvMatrix::solve
+ FatalErrorInFunction
+ << "Running extended addressing is not allowed when solving "
+ << fldName
+ << " Please choose a dynamicFvMesh without matrix adaptation"
+ << exit(FatalError);
+ }
+ else
+ {
+ if (debug)
+ {
+ Info<< "Interpolating field " << fldName << endl;
+ }
+
+ // Interpolate without bc update (since would trigger infinite
+ // recursion back to
+ // semiImplicitOversetFvPatchField<Type>::evaluate)
+ // The only problem is bcs that use the new cell values
+ // (e.g. zeroGradient, processor). These need to appear -after-
+ // the 'overset' bc.
+ mesh.interpolate
+ (
+ const_cast<Field<Type>&>
+ (
+ this->primitiveField()
+ )
+ );
+ }
+ }
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::updateInterfaceMatrix
+(
+ scalarField& result,
+ const bool add,
+ const scalarField& psiInternal,
+ const scalarField& coeffs,
+ const direction cmpt,
+ const Pstream::commsTypes
+) const
+{
+ if (debug)
+ {
+ Pout<< FUNCTION_NAME << " field " << this->internalField().name()
+ << " patch " << this->patch().name() << endl;
+ }
+
+ const oversetFvPatch& ovp = this->oversetPatch();
+
+ // Set all interpolated cells
+ if (ovp.master())
+ {
+ const labelListList& stencil = ovp.stencil();
+
+ if (stencil.size() != psiInternal.size())
+ {
+ FatalErrorInFunction << "psiInternal:" << psiInternal.size()
+ << " stencil:" << stencil.size() << exit(FatalError);
+ }
+
+ const mapDistribute& map = ovp.cellInterpolationMap();
+ const List<scalarList>& wghts = ovp.cellInterpolationWeights();
+ const labelList& cellIDs = ovp.interpolationCells();
+ //const scalarList& factor = ovp.cellInterpolationWeight();
+
+ // Since we're inside initEvaluate/evaluate there might be processor
+ // comms underway. Change the tag we use.
+ scalarField work(psiInternal);
+ map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+
+ forAll(cellIDs, i)
+ {
+ label celli = cellIDs[i];
+
+ const scalarList& w = wghts[celli];
+ const labelList& nbrs = stencil[celli];
+
+ //scalar f = factor[celli];
+
+ scalar s(0.0);
+ forAll(nbrs, nbrI)
+ {
+ s += w[nbrI]*work[nbrs[nbrI]];
+ }
+
+ //Pout<< "cell:" << celli << " interpolated value:" << s << endl;
+ result[celli] = s; //(1.0-f)*result[celli] + f*s;
+ }
+ }
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::write(Ostream& os) const
+{
+ zeroGradientFvPatchField<Type>::write(os);
+ // Make sure to write the value for ease of postprocessing.
+ this->writeEntry("value", os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.H b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.H
new file mode 100644
index 0000000000000000000000000000000000000000..f01a412ec1fc62c020dfe43318bfbd7fcd6d527b
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.H
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+Class
+ Foam::semiImplicitOversetFvPatchField
+
+Group
+ grpCoupledBoundaryConditions
+
+Description
+ Boundary condition for use on overset patches. Implements update of
+ interpolated cells inside the linear solvers.
+
+SeeAlso
+
+SourceFiles
+ semiImplicitOversetFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetFvPatchField_H
+#define semiImplicitOversetFvPatchField_H
+
+#include "oversetFvPatch.H"
+#include "zeroGradientFvPatchField.H"
+#include "LduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class semiImplicitOversetFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class semiImplicitOversetFvPatchField
+:
+ public LduInterfaceField<Type>,
+ public zeroGradientFvPatchField<Type>
+{
+protected:
+
+ // Protected data
+
+ //- Local reference cast into the overset patch
+ const oversetFvPatch& oversetPatch_;
+
+ //- Master patch ID
+ mutable label masterPatchID_;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("semiImplicitOverset");
+
+
+ // Constructors
+
+ //- Construct from patch and internal field
+ semiImplicitOversetFvPatchField
+ (
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&
+ );
+
+ //- Construct from patch, internal field and dictionary
+ semiImplicitOversetFvPatchField
+ (
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&,
+ const dictionary&
+ );
+
+ //- Construct by mapping given semiImplicitOversetFvPatchField onto a
+ // new patch
+ semiImplicitOversetFvPatchField
+ (
+ const semiImplicitOversetFvPatchField<Type>&,
+ const fvPatch&,
+ const DimensionedField<Type, volMesh>&,
+ const fvPatchFieldMapper&
+ );
+
+ //- Construct as copy
+ semiImplicitOversetFvPatchField
+ (
+ const semiImplicitOversetFvPatchField<Type>&
+ );
+
+ //- Construct and return a clone
+ virtual tmp<fvPatchField<Type> > clone() const
+ {
+ return tmp<fvPatchField<Type> >
+ (
+ new semiImplicitOversetFvPatchField<Type>(*this)
+ );
+ }
+
+ //- Construct as copy setting internal field reference
+ semiImplicitOversetFvPatchField
+ (
+ const semiImplicitOversetFvPatchField<Type>&,
+ const DimensionedField<Type, volMesh>&
+ );
+
+ //- Construct and return a clone setting internal field reference
+ virtual tmp<fvPatchField<Type> > clone
+ (
+ const DimensionedField<Type, volMesh>& iF
+ ) const
+ {
+ return tmp<fvPatchField<Type> >
+ (
+ new semiImplicitOversetFvPatchField<Type>(*this, iF)
+ );
+ }
+
+
+ // Member functions
+
+ // Access
+
+ //- Return local reference cast into the cyclic AMI patch
+ const oversetFvPatch& oversetPatch() const
+ {
+ return oversetPatch_;
+ }
+
+
+ // Evaluation functions
+
+ //- Evaluate the patch field
+ virtual void evaluate(const Pstream::commsTypes commsType);
+
+ //- Initialise the evaluation of the patch field
+ virtual void initEvaluate(const Pstream::commsTypes commsType);
+
+ //- Update result field based on interface functionality
+ virtual void updateInterfaceMatrix
+ (
+ scalarField& result,
+ const bool add,
+ const scalarField& psiInternal,
+ const scalarField& coeffs,
+ const direction cmpt,
+ const Pstream::commsTypes commsType
+ ) const;
+
+ //- Update result field based on interface functionality
+ virtual void updateInterfaceMatrix
+ (
+ Field<Type>&,
+ const bool add,
+ const Field<Type>&,
+ const scalarField&,
+ const Pstream::commsTypes commsType
+ ) const
+ {
+ NotImplemented;
+ }
+
+
+ // I-O
+
+ //- Write
+ virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+# include "semiImplicitOversetFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.C b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.C
new file mode 100644
index 0000000000000000000000000000000000000000..1471e28ca46542c307ab22d146e6d41d8d69a798
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2016-2017 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 "semiImplicitOversetFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(semiImplicitOverset);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.H b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..e943538791499b24b75237ad34ccebef975c66fb
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetFvPatchFields_H
+#define semiImplicitOversetFvPatchFields_H
+
+#include "semiImplicitOversetFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(semiImplicitOverset);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C b/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C
new file mode 100644
index 0000000000000000000000000000000000000000..cf624a85ac6eff547e92002c56b4ae510e63525a
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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 "semiImplicitOversetGAMGInterfaceField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(semiImplicitOversetGAMGInterfaceField, 0);
+ addToRunTimeSelectionTable
+ (
+ GAMGInterfaceField,
+ semiImplicitOversetGAMGInterfaceField,
+ lduInterfaceField
+ );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::semiImplicitOversetGAMGInterfaceField::
+semiImplicitOversetGAMGInterfaceField
+(
+ const GAMGInterface& GAMGCp,
+ const lduInterfaceField& fineInterface
+)
+:
+ oversetGAMGInterfaceField(GAMGCp, fineInterface)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::semiImplicitOversetGAMGInterfaceField::
+~semiImplicitOversetGAMGInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::semiImplicitOversetGAMGInterfaceField::updateInterfaceMatrix
+(
+ scalarField& result,
+ const bool add,
+ const scalarField& psiInternal,
+ const scalarField& coeffs,
+ const direction cmpt,
+ const Pstream::commsTypes commsType
+) const
+{
+ // Add remote values
+ if (oversetInterface_.master())
+ {
+ const labelListList& stencil = oversetInterface_.stencil();
+
+ if (stencil.size() != psiInternal.size())
+ {
+ FatalErrorInFunction << "psiInternal:" << psiInternal.size()
+ << " stencil:" << stencil.size() << exit(FatalError);
+ }
+
+ const mapDistribute& map = oversetInterface_.cellInterpolationMap();
+ const List<scalarList>& wghts =
+ oversetInterface_.cellInterpolationWeights();
+ const labelList& cellIDs = oversetInterface_.interpolationCells();
+ const scalarList& factor = oversetInterface_.cellInterpolationWeight();
+ const scalarField& normalisation = oversetInterface_.normalisation();
+
+ scalarField work(psiInternal);
+ map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+
+ forAll(cellIDs, i)
+ {
+ label celli = cellIDs[i];
+
+ const scalarList& w = wghts[celli];
+ const labelList& nbrs = stencil[celli];
+ const scalar f = factor[celli];
+
+ scalar s(0.0);
+ forAll(nbrs, nbrI)
+ {
+ label slotI = nbrs[nbrI];
+ s += w[nbrI]*work[slotI];
+ }
+ if (add)
+ {
+ s = -1.0*s;
+ }
+
+ result[celli] += normalisation[celli]*f*s;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.H b/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.H
new file mode 100644
index 0000000000000000000000000000000000000000..c80bb384575b688deb970a5665cba0a03c521cd4
--- /dev/null
+++ b/src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2017 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/>.
+
+Class
+ Foam::semiImplicitOversetGAMGInterfaceField
+
+Description
+ GAMG agglomerated processor interface field.
+
+SourceFiles
+ processorGAMGInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetGAMGInterfaceField_H
+#define semiImplicitOversetGAMGInterfaceField_H
+
+#include "oversetGAMGInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class semiImplicitOversetGAMGInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+class semiImplicitOversetGAMGInterfaceField
+:
+ public oversetGAMGInterfaceField
+{
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ semiImplicitOversetGAMGInterfaceField
+ (
+ const semiImplicitOversetGAMGInterfaceField&
+ );
+
+ //- Disallow default bitwise assignment
+ void operator=(const semiImplicitOversetGAMGInterfaceField&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("semiImplicitOverset");
+
+
+ // Constructors
+
+ //- Construct from GAMG interface and fine level interface field
+ semiImplicitOversetGAMGInterfaceField
+ (
+ const GAMGInterface& GAMGCp,
+ const lduInterfaceField& fineInterface
+ );
+
+
+ //- Destructor
+ virtual ~semiImplicitOversetGAMGInterfaceField();
+
+
+ // Member Functions
+
+ //- Update result field based on interface functionality
+ virtual void updateInterfaceMatrix
+ (
+ scalarField& result,
+ const bool add,
+ const scalarField& psiInternal,
+ const scalarField& coeffs,
+ const direction cmpt,
+ const Pstream::commsTypes commsType
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSchemes b/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSchemes
index d31d118e8aa254accf6dfed3d899f1794d7e7c0b..2d4f3b5ee274849ee2f6277b024af248d687f9f9 100644
--- a/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSchemes
+++ b/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSchemes
@@ -56,7 +56,10 @@ snGradSchemes
oversetInterpolation
{
- method inverseDistance;
+ method inverseDistance;
+
+ searchBox (0 0 0)(0.01 0.01 0.01);
+ searchBoxDivisions (100 100 1);
}
fluxRequired
diff --git a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/decomposeParDict b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/decomposeParDict
index f1b18226f652df92338ddc7a533bc2ef1b91a1ef..602c82d2c25a5c2589f61729d55fea62c8bcdce8 100644
--- a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/decomposeParDict
+++ b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/decomposeParDict
@@ -16,11 +16,11 @@ FoamFile
numberOfSubdomains 8;
-method hierarchical;
+method scotch; //hierarchical;
coeffs
{
- n (8 1 1);
+ n (4 2 1);
//delta 0.001; // default=0.001
//order xyz; // default=xzy
}
diff --git a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/fvSolution b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/fvSolution
index 71cf71b4055b1515b1cc3ee9dc493bac665f7088..84de10a190be830b9f3dd55543d1bb7318032dba 100644
--- a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/fvSolution
+++ b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/fvSolution
@@ -60,6 +60,7 @@ solvers
SIMPLE
{
+ momentumPredictor true;
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
diff --git a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
index 04bfe4172aac70c8c36cf696d63b290cfb3af21b..f1aaf2b011a9244978117b7e2e8c5e48963b66f1 100644
--- a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
+++ b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution
@@ -29,8 +29,8 @@ solvers
"alpha.water.*"
{
- nAlphaCorr 2;
- nAlphaSubCycles 1;
+ nAlphaCorr 3;
+ nAlphaSubCycles 2;
cAlpha 1;
icAlpha 0;
@@ -82,8 +82,12 @@ PIMPLE
nCorrectors 2;
nNonOrthogonalCorrectors 0;
- ddtCorr yes;
- correctPhi no;
+ ddtCorr yes;
+ correctPhi no;
+ massFluxInterpolation no;
+
+ moveMeshOuterCorrectors no;
+ turbOnFinalIterOnly no;
oversetAdjustPhi no;