ENH: Delete unused ot updated files
This commit is contained in:
sergio
@ -1,4 +0,0 @@
|
||||
bool ddtCorr
|
||||
(
|
||||
pimple.dict().getOrDefault("ddtCorr", true)
|
||||
);
|
||||
@ -1,9 +0,0 @@
|
||||
#include "readTimeControls.H"
|
||||
|
||||
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
|
||||
|
||||
checkMeshCourantNo =
|
||||
pimple.dict().getOrDefault("checkMeshCourantNo", false);
|
||||
|
||||
|
||||
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
|
||||
@ -1,9 +0,0 @@
|
||||
#include "readTimeControls.H"
|
||||
|
||||
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
|
||||
|
||||
checkMeshCourantNo =
|
||||
pimple.dict().getOrDefault("checkMeshCourantNo", false);
|
||||
|
||||
|
||||
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
|
||||
@ -1,26 +0,0 @@
|
||||
#include "createTimeControls.H"
|
||||
|
||||
bool correctPhi
|
||||
(
|
||||
pimple.dict().getOrDefault("correctPhi", false)
|
||||
);
|
||||
|
||||
bool checkMeshCourantNo
|
||||
(
|
||||
pimple.dict().getOrDefault("checkMeshCourantNo", false)
|
||||
);
|
||||
|
||||
bool massFluxInterpolation
|
||||
(
|
||||
pimple.dict().getOrDefault("massFluxInterpolation", false)
|
||||
);
|
||||
|
||||
bool adjustFringe
|
||||
(
|
||||
pimple.dict().getOrDefault("oversetAdjustPhi", false)
|
||||
);
|
||||
|
||||
bool ddtCorr
|
||||
(
|
||||
pimple.dict().getOrDefault("ddtCorr", true)
|
||||
);
|
||||
@ -1,273 +0,0 @@
|
||||
// Interpolation used
|
||||
interpolationCellPoint<vector> UInterpolator(HbyA);
|
||||
|
||||
// Determine faces on outside of interpolated cells
|
||||
bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
|
||||
bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
|
||||
|
||||
// Determine donor cells
|
||||
labelListList donorCell(mesh.nInternalFaces());
|
||||
|
||||
scalarListList weightCellCells(mesh.nInternalFaces());
|
||||
|
||||
// Interpolated HbyA faces
|
||||
vectorField UIntFaces(mesh.nInternalFaces(), Zero);
|
||||
|
||||
// Determine receptor neighbour cells
|
||||
labelList receptorNeigCell(mesh.nInternalFaces(), -1);
|
||||
|
||||
{
|
||||
const cellCellStencilObject& overlap = Stencil::New(mesh);
|
||||
const labelList& cellTypes = overlap.cellTypes();
|
||||
const labelIOList& zoneID = overlap.zoneID();
|
||||
|
||||
label nZones = gMax(zoneID)+1;
|
||||
PtrList<fvMeshSubset> meshParts(nZones);
|
||||
labelList nCellsPerZone(nZones, Zero);
|
||||
|
||||
// A mesh subset for each zone
|
||||
forAll(meshParts, zonei)
|
||||
{
|
||||
meshParts.set
|
||||
(
|
||||
zonei,
|
||||
// Select cells where the zoneID == zonei
|
||||
new fvMeshSubset(mesh, zonei, zoneID)
|
||||
);
|
||||
}
|
||||
|
||||
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
|
||||
{
|
||||
label ownType = cellTypes[mesh.faceOwner()[faceI]];
|
||||
label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
|
||||
if
|
||||
(
|
||||
ownType == cellCellStencil::INTERPOLATED
|
||||
&& neiType == cellCellStencil::CALCULATED
|
||||
)
|
||||
{
|
||||
isOwnerInterpolatedFace.set(faceI);
|
||||
|
||||
const vector& fc = mesh.faceCentres()[faceI];
|
||||
|
||||
for (label zoneI = 0; zoneI < nZones; zoneI++)
|
||||
{
|
||||
if (zoneI != zoneID[mesh.faceOwner()[faceI]])
|
||||
{
|
||||
const fvMesh& partMesh = meshParts[zoneI].subMesh();
|
||||
const labelList& cellMap = meshParts[zoneI].cellMap();
|
||||
label cellI = partMesh.findCell(fc);
|
||||
|
||||
if (cellI != -1)
|
||||
{
|
||||
// Determine weights
|
||||
labelList stencil(partMesh.cellCells()[cellI]);
|
||||
|
||||
stencil.append(cellI);
|
||||
|
||||
label st = stencil.size();
|
||||
|
||||
donorCell[faceI].setSize(st);
|
||||
|
||||
weightCellCells[faceI].setSize(st);
|
||||
|
||||
scalarField weights(st);
|
||||
|
||||
forAll(stencil, i)
|
||||
{
|
||||
scalar d = mag
|
||||
(
|
||||
partMesh.cellCentres()[stencil[i]]
|
||||
- fc
|
||||
);
|
||||
weights[i] = 1.0/d;
|
||||
donorCell[faceI][i] = cellMap[stencil[i]];
|
||||
}
|
||||
weights /= sum(weights);
|
||||
|
||||
weightCellCells[faceI] = weights;
|
||||
|
||||
forAll(stencil, i)
|
||||
{
|
||||
UIntFaces[faceI] +=
|
||||
weightCellCells[faceI][i]
|
||||
*UInterpolator.interpolate
|
||||
(
|
||||
fc,
|
||||
donorCell[faceI][i]
|
||||
);
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
|
||||
}
|
||||
else if
|
||||
(
|
||||
ownType == cellCellStencil::CALCULATED
|
||||
&& neiType == cellCellStencil::INTERPOLATED
|
||||
)
|
||||
{
|
||||
isNeiInterpolatedFace.set(faceI);
|
||||
|
||||
const vector& fc = mesh.faceCentres()[faceI];
|
||||
for (label zoneI = 0; zoneI < nZones; zoneI++)
|
||||
{
|
||||
if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
|
||||
{
|
||||
const fvMesh& partMesh = meshParts[zoneI].subMesh();
|
||||
const labelList& cellMap = meshParts[zoneI].cellMap();
|
||||
label cellI = partMesh.findCell(fc);
|
||||
|
||||
if (cellI != -1)
|
||||
{
|
||||
// Determine weights
|
||||
labelList stencil(partMesh.cellCells()[cellI]);
|
||||
|
||||
stencil.append(cellI);
|
||||
|
||||
label st = stencil.size();
|
||||
|
||||
donorCell[faceI].setSize(st);
|
||||
|
||||
weightCellCells[faceI].setSize(st);
|
||||
|
||||
scalarField weights(st);
|
||||
|
||||
forAll(stencil, i)
|
||||
{
|
||||
scalar d = mag
|
||||
(
|
||||
partMesh.cellCentres()[stencil[i]]
|
||||
- fc
|
||||
);
|
||||
weights[i] = 1.0/d;
|
||||
donorCell[faceI][i] = cellMap[stencil[i]];
|
||||
}
|
||||
weights /= sum(weights);
|
||||
|
||||
weightCellCells[faceI] = weights;
|
||||
|
||||
forAll(stencil, i)
|
||||
{
|
||||
UIntFaces[faceI] +=
|
||||
weightCellCells[faceI][i]
|
||||
*UInterpolator.interpolate
|
||||
(
|
||||
fc,
|
||||
donorCell[faceI][i]
|
||||
);
|
||||
}
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// contravariant U
|
||||
vectorField U1Contrav(mesh.nInternalFaces(), Zero);
|
||||
|
||||
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
|
||||
|
||||
forAll(isNeiInterpolatedFace, faceI)
|
||||
{
|
||||
label cellId = -1;
|
||||
if (isNeiInterpolatedFace.test(faceI))
|
||||
{
|
||||
cellId = mesh.faceNeighbour()[faceI];
|
||||
}
|
||||
else if (isOwnerInterpolatedFace.test(faceI))
|
||||
{
|
||||
cellId = mesh.faceOwner()[faceI];
|
||||
}
|
||||
|
||||
if (cellId != -1)
|
||||
{
|
||||
const vector& n = faceNormals[faceI];
|
||||
vector n1(Zero);
|
||||
|
||||
// 2-D cases
|
||||
if (mesh.nSolutionD() == 2)
|
||||
{
|
||||
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
|
||||
{
|
||||
if (mesh.geometricD()[cmpt] == -1)
|
||||
{
|
||||
switch (cmpt)
|
||||
{
|
||||
case vector::X:
|
||||
{
|
||||
n1 = vector(0, n.z(), -n.y());
|
||||
break;
|
||||
}
|
||||
|
||||
case vector::Y:
|
||||
{
|
||||
n1 = vector(n.z(), 0, -n.x());
|
||||
break;
|
||||
}
|
||||
|
||||
case vector::Z:
|
||||
{
|
||||
n1 = vector(n.y(), -n.x(), 0);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (mesh.nSolutionD() == 3)
|
||||
{
|
||||
//Determine which is the primary direction
|
||||
if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
|
||||
{
|
||||
n1 = vector(n.y(), -n.x(), 0);
|
||||
}
|
||||
else if (mag(n.y()) > mag(n.z()))
|
||||
{
|
||||
n1 = vector(0, n.z(), -n.y());
|
||||
}
|
||||
else
|
||||
{
|
||||
n1 = vector(-n.z(), 0, n.x());
|
||||
}
|
||||
}
|
||||
n1.normalise();
|
||||
|
||||
const vector n2 = normalised(n ^ n1);
|
||||
|
||||
tensor rot =
|
||||
tensor
|
||||
(
|
||||
n.x() ,n.y(), n.z(),
|
||||
n1.x() ,n1.y(), n1.z(),
|
||||
n2.x() ,n2.y(), n2.z()
|
||||
);
|
||||
|
||||
// tensor rot =
|
||||
// tensor
|
||||
// (
|
||||
// n & x ,n & y, n & z,
|
||||
// n1 & x ,n1 & y, n1 & z,
|
||||
// n2 & x ,n2 & y, n2 & z
|
||||
// );
|
||||
|
||||
U1Contrav[faceI].x() =
|
||||
2*transform(rot, UIntFaces[faceI]).x()
|
||||
- transform(rot, HbyA[receptorNeigCell[faceI]]).x();
|
||||
|
||||
U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
|
||||
|
||||
U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
|
||||
|
||||
HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
|
||||
}
|
||||
}
|
||||
@ -1,10 +0,0 @@
|
||||
#include "readTimeControls.H"
|
||||
|
||||
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
|
||||
|
||||
checkMeshCourantNo = pimple.dict().getOrDefault("checkMeshCourantNo", false);
|
||||
|
||||
massFluxInterpolation =
|
||||
pimple.dict().getOrDefault("massFluxInterpolation", false);
|
||||
|
||||
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
|
||||
@ -1,30 +0,0 @@
|
||||
bool correctPhi
|
||||
(
|
||||
pimple.dict().getOrDefault("correctPhi", true)
|
||||
);
|
||||
|
||||
bool checkMeshCourantNo
|
||||
(
|
||||
pimple.dict().getOrDefault("checkMeshCourantNo", false)
|
||||
);
|
||||
|
||||
bool moveMeshOuterCorrectors
|
||||
(
|
||||
pimple.dict().getOrDefault("moveMeshOuterCorrectors", false)
|
||||
);
|
||||
|
||||
|
||||
bool massFluxInterpolation
|
||||
(
|
||||
pimple.dict().getOrDefault("massFluxInterpolation", false)
|
||||
);
|
||||
|
||||
bool adjustFringe
|
||||
(
|
||||
pimple.dict().getOrDefault("oversetAdjustPhi", false)
|
||||
);
|
||||
|
||||
bool ddtCorr
|
||||
(
|
||||
pimple.dict().getOrDefault("ddtCorr", true)
|
||||
);
|
||||
@ -1,16 +0,0 @@
|
||||
#include "readTimeControls.H"
|
||||
|
||||
correctPhi = pimple.dict().getOrDefault("correctPhi", false);
|
||||
|
||||
checkMeshCourantNo =
|
||||
pimple.dict().getOrDefault("checkMeshCourantNo", false);
|
||||
|
||||
moveMeshOuterCorrectors =
|
||||
pimple.dict().getOrDefault("moveMeshOuterCorrectors", false);
|
||||
|
||||
massFluxInterpolation =
|
||||
pimple.dict().getOrDefault("massFluxInterpolation", false);
|
||||
|
||||
ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
|
||||
|
||||
adjustFringe = pimple.dict().getOrDefault("oversetAdjustPhi", false);
|
||||
@ -1,795 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2014-2020 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
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 "dynamicOversetFvMesh.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
#include "cellCellStencilObject.H"
|
||||
#include "zeroGradientFvPatchFields.H"
|
||||
#include "lduPrimitiveProcessorInterface.H"
|
||||
#include "globalIndex.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
|
||||
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
|
||||
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, doInit);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
bool Foam::dynamicOversetFvMesh::updateAddressing() const
|
||||
{
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
|
||||
// The (processor-local part of the) stencil determines the local
|
||||
// faces to add to the matrix. tbd: parallel
|
||||
const labelListList& stencil = overlap.cellStencil();
|
||||
|
||||
// Get the base addressing
|
||||
//const lduAddressing& baseAddr = dynamicMotionSolverFvMesh::lduAddr();
|
||||
const lduAddressing& baseAddr = dynamicFvMesh::lduAddr();
|
||||
|
||||
// Add to the base addressing
|
||||
labelList lowerAddr;
|
||||
labelList upperAddr;
|
||||
label nExtraFaces;
|
||||
|
||||
const globalIndex globalNumbering(baseAddr.size());
|
||||
labelListList localFaceCells;
|
||||
labelListList remoteFaceCells;
|
||||
|
||||
labelList globalCellIDs(overlap.cellInterpolationMap().constructSize());
|
||||
forAll(baseAddr, cellI)
|
||||
{
|
||||
globalCellIDs[cellI] = globalNumbering.toGlobal(cellI);
|
||||
}
|
||||
overlap.cellInterpolationMap().distribute(globalCellIDs);
|
||||
|
||||
|
||||
reverseFaceMap_ = fvMeshPrimitiveLduAddressing::addAddressing
|
||||
(
|
||||
baseAddr,
|
||||
stencil,
|
||||
nExtraFaces,
|
||||
lowerAddr,
|
||||
upperAddr,
|
||||
stencilFaces_,
|
||||
globalNumbering,
|
||||
globalCellIDs,
|
||||
localFaceCells,
|
||||
remoteFaceCells
|
||||
);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "dynamicOversetFvMesh::update() : extended addressing from"
|
||||
<< " nFaces:" << baseAddr.lowerAddr().size()
|
||||
<< " to nFaces:" << lowerAddr.size()
|
||||
<< " nExtraFaces:" << nExtraFaces << endl;
|
||||
}
|
||||
|
||||
|
||||
// Now for the tricky bits. We want to hand out processor faces according
|
||||
// to the localFaceCells/remoteFaceCells. Ultimately we need
|
||||
// per entry in stencil:
|
||||
// - the patch (or -1 for internal faces)
|
||||
// - the face (is either an internal face index or a patch face index)
|
||||
|
||||
stencilPatches_.setSize(stencilFaces_.size());
|
||||
|
||||
// Per processor to owner (local)/neighbour (remote)
|
||||
List<DynamicList<label>> procOwner(Pstream::nProcs());
|
||||
List<DynamicList<label>> dynProcNeighbour(Pstream::nProcs());
|
||||
forAll(stencil, celli)
|
||||
{
|
||||
const labelList& nbrs = stencil[celli];
|
||||
stencilPatches_[celli].setSize(nbrs.size());
|
||||
stencilPatches_[celli] = -1;
|
||||
|
||||
forAll(nbrs, nbri)
|
||||
{
|
||||
if (stencilFaces_[celli][nbri] == -1)
|
||||
{
|
||||
const label nbrCelli = nbrs[nbri];
|
||||
label globalNbr = globalCellIDs[nbrCelli];
|
||||
label proci = globalNumbering.whichProcID(globalNbr);
|
||||
label remoteCelli = globalNumbering.toLocal(proci, globalNbr);
|
||||
|
||||
// Overwrite the face to be a patch face
|
||||
stencilFaces_[celli][nbri] = procOwner[proci].size();
|
||||
stencilPatches_[celli][nbri] = proci;
|
||||
procOwner[proci].append(celli);
|
||||
dynProcNeighbour[proci].append(remoteCelli);
|
||||
|
||||
//Pout<< "From neighbour proc:" << proci
|
||||
// << " allocating patchFace:" << stencilFaces_[celli][nbri]
|
||||
// << " to get remote cell " << remoteCelli
|
||||
// << " onto local cell " << celli << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
labelListList procNeighbour(dynProcNeighbour.size());
|
||||
forAll(procNeighbour, i)
|
||||
{
|
||||
procNeighbour[i] = std::move(dynProcNeighbour[i]);
|
||||
}
|
||||
labelListList mySendCells;
|
||||
Pstream::exchange<labelList, label>(procNeighbour, mySendCells);
|
||||
|
||||
label nbri = 0;
|
||||
forAll(procOwner, proci)
|
||||
{
|
||||
if (procOwner[proci].size())
|
||||
{
|
||||
nbri++;
|
||||
}
|
||||
if (mySendCells[proci].size())
|
||||
{
|
||||
nbri++;
|
||||
}
|
||||
}
|
||||
remoteStencilInterfaces_.setSize(nbri);
|
||||
nbri = 0;
|
||||
|
||||
// E.g. if proc1 needs some data from proc2 and proc2 needs some data from
|
||||
// proc1. We first want the pair : proc1 receive and proc2 send
|
||||
// and then the pair : proc1 send, proc2 receive
|
||||
|
||||
|
||||
labelList procToInterface(Pstream::nProcs(), -1);
|
||||
|
||||
forAll(procOwner, proci)
|
||||
{
|
||||
if (proci < Pstream::myProcNo() && procOwner[proci].size())
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Adding interface " << nbri
|
||||
<< " to receive my " << procOwner[proci].size()
|
||||
<< " from " << proci << endl;
|
||||
}
|
||||
procToInterface[proci] = nbri;
|
||||
remoteStencilInterfaces_.set
|
||||
(
|
||||
nbri++,
|
||||
new lduPrimitiveProcessorInterface
|
||||
(
|
||||
procOwner[proci],
|
||||
Pstream::myProcNo(),
|
||||
proci,
|
||||
tensorField(0),
|
||||
Pstream::msgType()+2
|
||||
)
|
||||
);
|
||||
}
|
||||
else if (proci > Pstream::myProcNo() && mySendCells[proci].size())
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Adding interface " << nbri
|
||||
<< " to send my " << mySendCells[proci].size()
|
||||
<< " to " << proci << endl;
|
||||
}
|
||||
remoteStencilInterfaces_.set
|
||||
(
|
||||
nbri++,
|
||||
new lduPrimitiveProcessorInterface
|
||||
(
|
||||
mySendCells[proci],
|
||||
Pstream::myProcNo(),
|
||||
proci,
|
||||
tensorField(0),
|
||||
Pstream::msgType()+2
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
forAll(procOwner, proci)
|
||||
{
|
||||
if (proci > Pstream::myProcNo() && procOwner[proci].size())
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Adding interface " << nbri
|
||||
<< " to receive my " << procOwner[proci].size()
|
||||
<< " from " << proci << endl;
|
||||
}
|
||||
procToInterface[proci] = nbri;
|
||||
remoteStencilInterfaces_.set
|
||||
(
|
||||
nbri++,
|
||||
new lduPrimitiveProcessorInterface
|
||||
(
|
||||
procOwner[proci],
|
||||
Pstream::myProcNo(),
|
||||
proci,
|
||||
tensorField(0),
|
||||
Pstream::msgType()+3
|
||||
)
|
||||
);
|
||||
}
|
||||
else if (proci < Pstream::myProcNo() && mySendCells[proci].size())
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "Adding interface " << nbri
|
||||
<< " to send my " << mySendCells[proci].size()
|
||||
<< " to " << proci << endl;
|
||||
}
|
||||
remoteStencilInterfaces_.set
|
||||
(
|
||||
nbri++,
|
||||
new lduPrimitiveProcessorInterface
|
||||
(
|
||||
mySendCells[proci],
|
||||
Pstream::myProcNo(),
|
||||
proci,
|
||||
tensorField(0),
|
||||
Pstream::msgType()+3
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Rewrite stencilPatches now we know the actual interface (procToInterface)
|
||||
for (auto& patches : stencilPatches_)
|
||||
{
|
||||
for (auto& interface : patches)
|
||||
{
|
||||
if (interface != -1)
|
||||
{
|
||||
interface = procToInterface[interface]+boundary().size();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Get addressing and interfaces of all interfaces
|
||||
|
||||
|
||||
List<const labelUList*> patchAddr;
|
||||
{
|
||||
const fvBoundaryMesh& fvp = boundary();
|
||||
|
||||
patchAddr.setSize(fvp.size() + remoteStencilInterfaces_.size());
|
||||
|
||||
//allInterfaces_ = dynamicMotionSolverFvMesh::interfaces();
|
||||
allInterfaces_ = dynamicFvMesh::interfaces();
|
||||
allInterfaces_.setSize(patchAddr.size());
|
||||
|
||||
forAll(fvp, patchI)
|
||||
{
|
||||
patchAddr[patchI] = &fvp[patchI].faceCells();
|
||||
}
|
||||
forAll(remoteStencilInterfaces_, i)
|
||||
{
|
||||
label patchI = fvp.size()+i;
|
||||
const lduPrimitiveProcessorInterface& pp =
|
||||
remoteStencilInterfaces_[i];
|
||||
|
||||
//Pout<< "at patch:" << patchI
|
||||
// << " have procPatch:" << pp.type()
|
||||
// << " from:" << pp.myProcNo()
|
||||
// << " to:" << pp.neighbProcNo()
|
||||
// << " with fc:" << pp.faceCells().size() << endl;
|
||||
|
||||
patchAddr[patchI] = &pp.faceCells();
|
||||
allInterfaces_.set(patchI, &pp);
|
||||
}
|
||||
}
|
||||
const lduSchedule ps
|
||||
(
|
||||
lduPrimitiveMesh::nonBlockingSchedule<processorLduInterface>
|
||||
(
|
||||
allInterfaces_
|
||||
)
|
||||
);
|
||||
|
||||
lduPtr_.reset
|
||||
(
|
||||
new fvMeshPrimitiveLduAddressing
|
||||
(
|
||||
nCells(),
|
||||
std::move(lowerAddr),
|
||||
std::move(upperAddr),
|
||||
patchAddr,
|
||||
ps
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
// Check
|
||||
if (debug)
|
||||
{
|
||||
const lduAddressing& addr = lduPtr_(); //this->lduAddr();
|
||||
|
||||
Pout<< "Adapted addressing:"
|
||||
<< " lower:" << addr.lowerAddr().size()
|
||||
<< " upper:" << addr.upperAddr().size() << endl;
|
||||
|
||||
// Using lduAddressing::patch
|
||||
forAll(patchAddr, patchI)
|
||||
{
|
||||
Pout<< " " << patchI << "\tpatchAddr:"
|
||||
<< addr.patchAddr(patchI).size()
|
||||
<< endl;
|
||||
}
|
||||
|
||||
// Using interfaces
|
||||
const lduInterfacePtrsList& iFaces = allInterfaces_;
|
||||
Pout<< "Adapted interFaces:" << iFaces.size() << endl;
|
||||
forAll(iFaces, patchI)
|
||||
{
|
||||
if (iFaces.set(patchI))
|
||||
{
|
||||
Pout<< " " << patchI << "\tinterface:"
|
||||
<< iFaces[patchI].type() << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
|
||||
(
|
||||
const labelList& types,
|
||||
const labelList& nbrTypes,
|
||||
const scalarField& norm,
|
||||
const scalarField& nbrNorm,
|
||||
const label celli,
|
||||
bitSet& isFront
|
||||
) const
|
||||
{
|
||||
const labelList& own = faceOwner();
|
||||
const labelList& nei = faceNeighbour();
|
||||
const cell& cFaces = cells()[celli];
|
||||
|
||||
scalar avg = 0.0;
|
||||
label n = 0;
|
||||
label nFront = 0;
|
||||
for (const label facei : cFaces)
|
||||
{
|
||||
if (isInternalFace(facei))
|
||||
{
|
||||
label nbrCelli = (own[facei] == celli ? nei[facei] : own[facei]);
|
||||
if (norm[nbrCelli] == -GREAT)
|
||||
{
|
||||
// Invalid neighbour. Add to front
|
||||
if (isFront.set(facei))
|
||||
{
|
||||
nFront++;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Valid neighbour. Add to average
|
||||
avg += norm[nbrCelli];
|
||||
n++;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (nbrNorm[facei-nInternalFaces()] == -GREAT)
|
||||
{
|
||||
if (isFront.set(facei))
|
||||
{
|
||||
nFront++;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
avg += nbrNorm[facei-nInternalFaces()];
|
||||
n++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (n > 0)
|
||||
{
|
||||
return avg/n;
|
||||
}
|
||||
else
|
||||
{
|
||||
return norm[celli];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::dynamicOversetFvMesh::writeAgglomeration
|
||||
(
|
||||
const GAMGAgglomeration& agglom
|
||||
) const
|
||||
{
|
||||
labelList cellToCoarse(identity(nCells()));
|
||||
labelListList coarseToCell(invertOneToMany(nCells(), cellToCoarse));
|
||||
|
||||
// Write initial agglomeration
|
||||
{
|
||||
volScalarField scalarAgglomeration
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"agglomeration",
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar(dimless, Zero)
|
||||
);
|
||||
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
|
||||
forAll(fld, celli)
|
||||
{
|
||||
fld[celli] = cellToCoarse[celli];
|
||||
}
|
||||
fld /= max(fld);
|
||||
correctBoundaryConditions
|
||||
<
|
||||
volScalarField,
|
||||
oversetFvPatchField<scalar>
|
||||
>(scalarAgglomeration.boundaryFieldRef(), false);
|
||||
scalarAgglomeration.write();
|
||||
|
||||
Info<< "Writing initial cell distribution to "
|
||||
<< this->time().timeName() << endl;
|
||||
}
|
||||
|
||||
|
||||
for (label level = 0; level < agglom.size(); level++)
|
||||
{
|
||||
const labelList& addr = agglom.restrictAddressing(level);
|
||||
label coarseSize = max(addr)+1;
|
||||
|
||||
Info<< "Level : " << level << endl
|
||||
<< returnReduce(addr.size(), sumOp<label>()) << endl
|
||||
<< " current size : "
|
||||
<< returnReduce(addr.size(), sumOp<label>()) << endl
|
||||
<< " agglomerated size : "
|
||||
<< returnReduce(coarseSize, sumOp<label>()) << endl;
|
||||
|
||||
forAll(addr, fineI)
|
||||
{
|
||||
const labelList& cellLabels = coarseToCell[fineI];
|
||||
forAll(cellLabels, i)
|
||||
{
|
||||
cellToCoarse[cellLabels[i]] = addr[fineI];
|
||||
}
|
||||
}
|
||||
coarseToCell = invertOneToMany(coarseSize, cellToCoarse);
|
||||
|
||||
// Write agglomeration
|
||||
{
|
||||
volScalarField scalarAgglomeration
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"agglomeration_" + Foam::name(level),
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar(dimless, Zero)
|
||||
);
|
||||
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
|
||||
forAll(fld, celli)
|
||||
{
|
||||
fld[celli] = cellToCoarse[celli];
|
||||
}
|
||||
//if (normalise)
|
||||
//{
|
||||
// fld /= max(fld);
|
||||
//}
|
||||
correctBoundaryConditions
|
||||
<
|
||||
volScalarField,
|
||||
oversetFvPatchField<scalar>
|
||||
>(scalarAgglomeration.boundaryFieldRef(), false);
|
||||
scalarAgglomeration.write();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::dynamicOversetFvMesh::dynamicOversetFvMesh
|
||||
(
|
||||
const IOobject& io,
|
||||
const bool doInit
|
||||
)
|
||||
:
|
||||
dynamicMotionSolverListFvMesh(io, doInit)
|
||||
{
|
||||
if (doInit)
|
||||
{
|
||||
init(false); // do not initialise lower levels
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool Foam::dynamicOversetFvMesh::init(const bool doInit)
|
||||
{
|
||||
if (doInit)
|
||||
{
|
||||
dynamicMotionSolverListFvMesh::init(doInit);
|
||||
}
|
||||
|
||||
active_ = false;
|
||||
|
||||
// Load stencil (but do not update)
|
||||
(void)Stencil::New(*this, false);
|
||||
|
||||
// Assume something changed
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
|
||||
{
|
||||
if (!active_)
|
||||
{
|
||||
//return dynamicMotionSolverFvMesh::lduAddr();
|
||||
return dynamicFvMesh::lduAddr();
|
||||
}
|
||||
if (!lduPtr_)
|
||||
{
|
||||
// Build extended addressing
|
||||
updateAddressing();
|
||||
}
|
||||
return *lduPtr_;
|
||||
}
|
||||
|
||||
|
||||
Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
|
||||
{
|
||||
if (!active_)
|
||||
{
|
||||
//return dynamicMotionSolverFvMesh::interfaces();
|
||||
return dynamicFvMesh::interfaces();
|
||||
}
|
||||
if (!lduPtr_)
|
||||
{
|
||||
// Build extended addressing
|
||||
updateAddressing();
|
||||
}
|
||||
return allInterfaces_;
|
||||
}
|
||||
|
||||
|
||||
const Foam::fvMeshPrimitiveLduAddressing&
|
||||
Foam::dynamicOversetFvMesh::primitiveLduAddr() const
|
||||
{
|
||||
if (!lduPtr_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Extended addressing not allocated" << abort(FatalError);
|
||||
}
|
||||
|
||||
return *lduPtr_;
|
||||
}
|
||||
|
||||
|
||||
bool Foam::dynamicOversetFvMesh::update()
|
||||
{
|
||||
//if (dynamicMotionSolverFvMesh::update())
|
||||
if (dynamicMotionSolverListFvMesh::update())
|
||||
{
|
||||
// Calculate the local extra faces for the interpolation. Note: could
|
||||
// let demand-driven lduAddr() trigger it but just to make sure.
|
||||
updateAddressing();
|
||||
|
||||
// Addressing and/or weights have changed. Make interpolated cells
|
||||
// up to date with donors
|
||||
interpolateFields();
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
Foam::word Foam::dynamicOversetFvMesh::baseName(const word& name)
|
||||
{
|
||||
if (name.ends_with("_0"))
|
||||
{
|
||||
return baseName(name.substr(0, name.size()-2));
|
||||
}
|
||||
|
||||
return name;
|
||||
}
|
||||
|
||||
|
||||
bool Foam::dynamicOversetFvMesh::interpolateFields()
|
||||
{
|
||||
// Add the stencil suppression list
|
||||
wordHashSet suppressed(Stencil::New(*this).nonInterpolatedFields());
|
||||
|
||||
// Use whatever the solver has set up as suppression list
|
||||
const dictionary* dictPtr
|
||||
(
|
||||
this->schemesDict().findDict("oversetInterpolationSuppressed")
|
||||
);
|
||||
if (dictPtr)
|
||||
{
|
||||
suppressed.insert(dictPtr->toc());
|
||||
}
|
||||
|
||||
interpolate<volScalarField>(suppressed);
|
||||
interpolate<volVectorField>(suppressed);
|
||||
interpolate<volSphericalTensorField>(suppressed);
|
||||
interpolate<volSymmTensorField>(suppressed);
|
||||
interpolate<volTensorField>(suppressed);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
bool Foam::dynamicOversetFvMesh::writeObject
|
||||
(
|
||||
IOstreamOption streamOpt,
|
||||
const bool valid
|
||||
) const
|
||||
{
|
||||
//bool ok = dynamicMotionSolverFvMesh::writeObject(streamOpt, valid);
|
||||
bool ok = dynamicMotionSolverListFvMesh::writeObject(streamOpt, valid);
|
||||
|
||||
// For postprocessing : write cellTypes and zoneID
|
||||
{
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
|
||||
const labelUList& cellTypes = overlap.cellTypes();
|
||||
|
||||
volScalarField volTypes
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"cellTypes",
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar(dimless, Zero),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
||||
forAll(volTypes.internalField(), cellI)
|
||||
{
|
||||
volTypes[cellI] = cellTypes[cellI];
|
||||
}
|
||||
volTypes.correctBoundaryConditions();
|
||||
volTypes.writeObject(streamOpt, valid);
|
||||
}
|
||||
{
|
||||
volScalarField volZoneID
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"zoneID",
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar(dimless, Zero),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
const labelIOList& zoneID = overlap.zoneID();
|
||||
|
||||
forAll(zoneID, cellI)
|
||||
{
|
||||
volZoneID[cellI] = zoneID[cellI];
|
||||
}
|
||||
volZoneID.correctBoundaryConditions();
|
||||
volZoneID.writeObject(streamOpt, valid);
|
||||
}
|
||||
if (debug)
|
||||
{
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
const labelIOList& zoneID = overlap.zoneID();
|
||||
const labelListList& cellStencil = overlap.cellStencil();
|
||||
|
||||
// Get remote zones
|
||||
labelList donorZoneID(zoneID);
|
||||
overlap.cellInterpolationMap().distribute(donorZoneID);
|
||||
|
||||
// Get remote cellCentres
|
||||
pointField cc(C());
|
||||
overlap.cellInterpolationMap().distribute(cc);
|
||||
|
||||
volScalarField volDonorZoneID
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"donorZoneID",
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar("minOne", dimless, scalar(-1)),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
||||
forAll(cellStencil, cellI)
|
||||
{
|
||||
const labelList& stencil = cellStencil[cellI];
|
||||
if (stencil.size())
|
||||
{
|
||||
volDonorZoneID[cellI] = donorZoneID[stencil[0]];
|
||||
for (label i = 1; i < stencil.size(); i++)
|
||||
{
|
||||
if (donorZoneID[stencil[i]] != volDonorZoneID[cellI])
|
||||
{
|
||||
WarningInFunction << "Mixed donor meshes for cell "
|
||||
<< cellI << " at " << C()[cellI]
|
||||
<< " donors:" << UIndirectList<point>(cc, stencil)
|
||||
<< endl;
|
||||
volDonorZoneID[cellI] = -2;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//- Do not correctBoundaryConditions since re-interpolates!
|
||||
//volDonorZoneID.correctBoundaryConditions();
|
||||
volDonorZoneID.writeObject(streamOpt, valid);
|
||||
}
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,392 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2015-2020 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
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::dynamicOversetFvMesh
|
||||
|
||||
Description
|
||||
dynamicFvMesh with support for overset meshes.
|
||||
|
||||
SourceFiles
|
||||
dynamicOversetFvMesh.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef dynamicOversetFvMesh_H
|
||||
#define dynamicOversetFvMesh_H
|
||||
|
||||
#include "dynamicMotionSolverListFvMesh.H"
|
||||
#include "labelIOList.H"
|
||||
#include "fvMeshPrimitiveLduAddressing.H"
|
||||
#include "lduInterfaceFieldPtrsList.H"
|
||||
#include "volFields.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
class mapDistribute;
|
||||
class lduPrimitiveProcessorInterface;
|
||||
class GAMGAgglomeration;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class dynamicOversetFvMesh Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class dynamicOversetFvMesh
|
||||
:
|
||||
public dynamicMotionSolverListFvMesh
|
||||
{
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Select base addressing (false) or locally stored extended
|
||||
// lduAddressing (true)
|
||||
mutable bool active_;
|
||||
|
||||
//- Extended addressing (extended with local interpolation stencils)
|
||||
mutable autoPtr<fvMeshPrimitiveLduAddressing> lduPtr_;
|
||||
|
||||
//- Added (processor)lduInterfaces for remote bits of stencil.
|
||||
//PtrList<const lduInterface> remoteStencilInterfaces_;
|
||||
mutable PtrList<const lduPrimitiveProcessorInterface>
|
||||
remoteStencilInterfaces_;
|
||||
|
||||
//- Interfaces for above mesh. Contains both original and
|
||||
// above added processorLduInterfaces
|
||||
mutable lduInterfacePtrsList allInterfaces_;
|
||||
|
||||
//- Corresponding faces (in above lduPtr) to the stencil
|
||||
mutable labelListList stencilFaces_;
|
||||
|
||||
//- Corresponding patches (in above lduPtr) to the stencil
|
||||
mutable labelListList stencilPatches_;
|
||||
|
||||
//- From old to new face labels
|
||||
mutable labelList reverseFaceMap_;
|
||||
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
//- Calculate the extended lduAddressing
|
||||
virtual bool updateAddressing() const;
|
||||
|
||||
//- Debug: print matrix
|
||||
template<class Type>
|
||||
void write
|
||||
(
|
||||
Ostream&,
|
||||
const fvMatrix<Type>&,
|
||||
const lduAddressing&,
|
||||
const lduInterfacePtrsList&
|
||||
) const;
|
||||
|
||||
//- Explicit interpolation of acceptor cells from donor cells
|
||||
template<class T>
|
||||
void interpolate(Field<T>& psi) const;
|
||||
|
||||
//- Explicit interpolation of acceptor cells from donor cells with
|
||||
// boundary condition handling
|
||||
template<class GeoField>
|
||||
void interpolate(GeoField& psi) const;
|
||||
|
||||
//- Helper: strip off trailing _0
|
||||
static word baseName(const word& name);
|
||||
|
||||
//- Explicit interpolation of all registered fields. Excludes
|
||||
// selected fields (and their old-time fields)
|
||||
template<class GeoField>
|
||||
void interpolate(const wordHashSet& suppressed);
|
||||
|
||||
//- 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;
|
||||
|
||||
//- Determine normalisation for interpolation. This equals the
|
||||
// original diagonal except stabilised for zero diagonals (possible
|
||||
// in hole cells)
|
||||
template<class Type>
|
||||
tmp<scalarField> normalisation(const fvMatrix<Type>& m) const;
|
||||
|
||||
//- Add interpolation to matrix (coefficients)
|
||||
template<class Type>
|
||||
void addInterpolation(fvMatrix<Type>&, const scalarField& norm) const;
|
||||
|
||||
//- Solve given dictionary with settings
|
||||
template<class Type>
|
||||
SolverPerformance<Type> solve(fvMatrix<Type>&, const dictionary&) const;
|
||||
|
||||
//- Debug: correct coupled bc
|
||||
template<class GeoField>
|
||||
static void correctCoupledBoundaryConditions(GeoField& fld);
|
||||
|
||||
//- Average norm of valid neighbours
|
||||
scalar cellAverage
|
||||
(
|
||||
const labelList& types,
|
||||
const labelList& nbrTypes,
|
||||
const scalarField& norm,
|
||||
const scalarField& nbrNorm,
|
||||
const label celli,
|
||||
bitSet& isFront
|
||||
) const;
|
||||
|
||||
//- Debug: dump agglomeration
|
||||
void writeAgglomeration
|
||||
(
|
||||
const GAMGAgglomeration& agglom
|
||||
) const;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- No copy construct
|
||||
dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const dynamicOversetFvMesh&) = delete;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("dynamicOversetFvMesh");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from IOobject
|
||||
dynamicOversetFvMesh(const IOobject& io, const bool doInit=true);
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~dynamicOversetFvMesh();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
|
||||
// Extended addressing
|
||||
|
||||
//- Return extended ldu addressing
|
||||
const fvMeshPrimitiveLduAddressing& primitiveLduAddr() const;
|
||||
|
||||
//- Return ldu addressing. If active: is (extended)
|
||||
// primitiveLduAddr
|
||||
virtual const lduAddressing& lduAddr() const;
|
||||
|
||||
//- Return a list of pointers for each patch
|
||||
// with only those pointing to interfaces being set. If active:
|
||||
// return additional remoteStencilInterfaces_
|
||||
virtual lduInterfacePtrsList interfaces() 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
|
||||
|
||||
// Explicit interpolation
|
||||
|
||||
virtual void interpolate(scalarField& psi) const
|
||||
{
|
||||
interpolate<scalar>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(vectorField& psi) const
|
||||
{
|
||||
interpolate<vector>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(sphericalTensorField& psi) const
|
||||
{
|
||||
interpolate<sphericalTensor>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(symmTensorField& psi) const
|
||||
{
|
||||
interpolate<symmTensor>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(tensorField& psi) const
|
||||
{
|
||||
interpolate<tensor>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(volScalarField& psi) const
|
||||
{
|
||||
interpolate<volScalarField>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(volVectorField& psi) const
|
||||
{
|
||||
interpolate<volVectorField>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(volSphericalTensorField& psi) const
|
||||
{
|
||||
interpolate<volSphericalTensorField>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(volSymmTensorField& psi) const
|
||||
{
|
||||
interpolate<volSymmTensorField>(psi);
|
||||
}
|
||||
|
||||
virtual void interpolate(volTensorField& psi) const
|
||||
{
|
||||
interpolate<volTensorField>(psi);
|
||||
}
|
||||
|
||||
|
||||
// Implicit interpolation (matrix manipulation)
|
||||
|
||||
//- Solve returning the solution statistics given convergence
|
||||
// tolerance. Use the given solver controls
|
||||
virtual SolverPerformance<scalar> solve
|
||||
(
|
||||
fvMatrix<scalar>& m,
|
||||
const dictionary& dict
|
||||
) const
|
||||
{
|
||||
return solve<scalar>(m, dict);
|
||||
}
|
||||
|
||||
//- Solve returning the solution statistics given convergence
|
||||
// tolerance. Use the given solver controls
|
||||
virtual SolverPerformance<vector> solve
|
||||
(
|
||||
fvMatrix<vector>& m,
|
||||
const dictionary& dict
|
||||
) const
|
||||
{
|
||||
return solve<vector>(m, dict);
|
||||
}
|
||||
|
||||
//- Solve returning the solution statistics given convergence
|
||||
// tolerance. Use the given solver controls
|
||||
virtual SolverPerformance<symmTensor> solve
|
||||
(
|
||||
fvMatrix<symmTensor>& m,
|
||||
const dictionary& dict
|
||||
) const
|
||||
{
|
||||
return solve<symmTensor>(m, dict);
|
||||
}
|
||||
|
||||
//- Solve returning the solution statistics given convergence
|
||||
// tolerance. Use the given solver controls
|
||||
virtual SolverPerformance<tensor> solve
|
||||
(
|
||||
fvMatrix<tensor>& m,
|
||||
const dictionary& dict
|
||||
) const
|
||||
{
|
||||
return solve<tensor>(m, dict);
|
||||
}
|
||||
|
||||
|
||||
//- Initialise all non-demand-driven data
|
||||
virtual bool init(const bool doInit);
|
||||
|
||||
//- Update the mesh for both mesh motion and topology change
|
||||
virtual bool update();
|
||||
|
||||
//- Update fields when mesh is updated
|
||||
virtual bool interpolateFields();
|
||||
|
||||
//- Write using stream options
|
||||
virtual bool writeObject
|
||||
(
|
||||
IOstreamOption streamOpt,
|
||||
const bool valid
|
||||
) const;
|
||||
|
||||
//- Debug: check halo swap is ok
|
||||
template<class GeoField>
|
||||
static void checkCoupledBC(const GeoField& fld);
|
||||
|
||||
//- Correct boundary conditions of certain type (typeOnly = true)
|
||||
// or explicitly not of the type (typeOnly = false)
|
||||
template<class GeoField, class PatchType>
|
||||
static void correctBoundaryConditions
|
||||
(
|
||||
typename GeoField::Boundary& bfld,
|
||||
const bool typeOnly
|
||||
);
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
# include "dynamicOversetFvMeshTemplates.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,961 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2014-2019 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
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"
|
||||
#include "fvMatrix.H"
|
||||
#include "cellCellStencilObject.H"
|
||||
#include "oversetFvPatchField.H"
|
||||
#include "calculatedProcessorFvPatchField.H"
|
||||
#include "lduInterfaceFieldPtrsList.H"
|
||||
#include "processorFvPatch.H"
|
||||
#include "syncTools.H"
|
||||
|
||||
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class T>
|
||||
void Foam::dynamicOversetFvMesh::interpolate(Field<T>& psi) const
|
||||
{
|
||||
const cellCellStencil& overlap = Stencil::New(*this);
|
||||
const labelListList& stencil = overlap.cellStencil();
|
||||
|
||||
if (stencil.size() != nCells())
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
const mapDistribute& map = overlap.cellInterpolationMap();
|
||||
const List<scalarList>& wghts = overlap.cellInterpolationWeights();
|
||||
const labelList& cellIDs = overlap.interpolationCells();
|
||||
const scalarList& factor = overlap.cellInterpolationWeight();
|
||||
|
||||
Field<T> work(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];
|
||||
|
||||
T s(pTraits<T>::zero);
|
||||
forAll(nbrs, nbrI)
|
||||
{
|
||||
s += w[nbrI]*work[nbrs[nbrI]];
|
||||
}
|
||||
//Pout<< "Interpolated value:" << s << endl;
|
||||
//T oldPsi = psi[celli];
|
||||
psi[celli] = (1.0-f)*psi[celli] + f*s;
|
||||
//Pout<< "psi was:" << oldPsi << " now:" << psi[celli] << endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::dynamicOversetFvMesh::interpolate(GeoField& psi) const
|
||||
{
|
||||
interpolate(psi.primitiveFieldRef());
|
||||
psi.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::dynamicOversetFvMesh::interpolate(const wordHashSet& suppressed)
|
||||
{
|
||||
auto flds(this->lookupClass<GeoField>());
|
||||
for (auto fldPtr : flds)
|
||||
{
|
||||
const word& name = fldPtr->name();
|
||||
if (!suppressed.found(baseName(name)))
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "dynamicOversetFvMesh::interpolate: interpolating : "
|
||||
<< name << endl;
|
||||
}
|
||||
interpolate(fldPtr->primitiveFieldRef());
|
||||
}
|
||||
else
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "dynamicOversetFvMesh::interpolate: skipping : " << name
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField, class PatchType>
|
||||
void Foam::dynamicOversetFvMesh::correctBoundaryConditions
|
||||
(
|
||||
typename GeoField::Boundary& bfld,
|
||||
const bool typeOnly
|
||||
)
|
||||
{
|
||||
const label nReq = Pstream::nRequests();
|
||||
|
||||
forAll(bfld, patchi)
|
||||
{
|
||||
if (typeOnly == (isA<PatchType>(bfld[patchi]) != nullptr))
|
||||
{
|
||||
bfld[patchi].initEvaluate(Pstream::defaultCommsType);
|
||||
}
|
||||
}
|
||||
|
||||
// Block for any outstanding requests
|
||||
if (Pstream::parRun())
|
||||
{
|
||||
Pstream::waitRequests(nReq);
|
||||
}
|
||||
|
||||
forAll(bfld, patchi)
|
||||
{
|
||||
if (typeOnly == (isA<PatchType>(bfld[patchi]) != nullptr))
|
||||
{
|
||||
bfld[patchi].evaluate(Pstream::defaultCommsType);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
|
||||
(
|
||||
const fvMatrix<Type>& m
|
||||
) const
|
||||
{
|
||||
// Determine normalisation. This is normally the original diagonal plus
|
||||
// remote contributions. This needs to be stabilised for hole cells
|
||||
// which can have a zero diagonal. Assume that if any component has
|
||||
// a non-zero diagonal the cell does not need stabilisation.
|
||||
tmp<scalarField> tnorm(tmp<scalarField>::New(m.diag()));
|
||||
scalarField& norm = tnorm.ref();
|
||||
|
||||
// Add remote coeffs to duplicate behaviour of fvMatrix::addBoundaryDiag
|
||||
const FieldField<Field, Type>& internalCoeffs = m.internalCoeffs();
|
||||
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
|
||||
{
|
||||
forAll(internalCoeffs, patchi)
|
||||
{
|
||||
const labelUList& fc = lduAddr().patchAddr(patchi);
|
||||
const Field<Type>& intCoeffs = internalCoeffs[patchi];
|
||||
const scalarField cmptCoeffs(intCoeffs.component(cmpt));
|
||||
forAll(fc, i)
|
||||
{
|
||||
norm[fc[i]] += cmptCoeffs[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Count number of problematic cells
|
||||
label nZeroDiag = 0;
|
||||
forAll(norm, celli)
|
||||
{
|
||||
const scalar& n = norm[celli];
|
||||
if (magSqr(n) < sqr(SMALL))
|
||||
{
|
||||
//Pout<< "For field " << m.psi().name()
|
||||
// << " have diagonal " << n << " for cell " << celli
|
||||
// << " at:" << cellCentres()[celli] << endl;
|
||||
nZeroDiag++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(nZeroDiag, sumOp<label>());
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "For field " << m.psi().name() << " have zero diagonals for "
|
||||
<< nZeroDiag << " cells" << endl;
|
||||
}
|
||||
|
||||
if (nZeroDiag > 0)
|
||||
{
|
||||
// Walk out the norm across hole cells
|
||||
|
||||
const labelList& own = faceOwner();
|
||||
const labelList& nei = faceNeighbour();
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
const labelUList& types = overlap.cellTypes();
|
||||
|
||||
label nHoles = 0;
|
||||
scalarField extrapolatedNorm(norm);
|
||||
forAll(types, celli)
|
||||
{
|
||||
if (types[celli] == cellCellStencil::HOLE)
|
||||
{
|
||||
extrapolatedNorm[celli] = -GREAT;
|
||||
nHoles++;
|
||||
}
|
||||
}
|
||||
|
||||
bitSet isFront(nFaces());
|
||||
for (label facei = 0; facei < nInternalFaces(); facei++)
|
||||
{
|
||||
label ownType = types[own[facei]];
|
||||
label neiType = types[nei[facei]];
|
||||
if
|
||||
(
|
||||
(ownType == cellCellStencil::HOLE)
|
||||
!= (neiType == cellCellStencil::HOLE)
|
||||
)
|
||||
{
|
||||
isFront.set(facei);
|
||||
}
|
||||
}
|
||||
labelList nbrTypes;
|
||||
syncTools::swapBoundaryCellList(*this, types, nbrTypes);
|
||||
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
|
||||
{
|
||||
label ownType = types[own[facei]];
|
||||
label neiType = nbrTypes[facei-nInternalFaces()];
|
||||
if
|
||||
(
|
||||
(ownType == cellCellStencil::HOLE)
|
||||
!= (neiType == cellCellStencil::HOLE)
|
||||
)
|
||||
{
|
||||
isFront.set(facei);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
while (true)
|
||||
{
|
||||
scalarField nbrNorm;
|
||||
syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
|
||||
|
||||
bitSet newIsFront(nFaces());
|
||||
scalarField newNorm(extrapolatedNorm);
|
||||
|
||||
label nChanged = 0;
|
||||
for (const label facei : isFront)
|
||||
{
|
||||
if (extrapolatedNorm[own[facei]] == -GREAT)
|
||||
{
|
||||
// Average owner cell, add faces to newFront
|
||||
newNorm[own[facei]] = cellAverage
|
||||
(
|
||||
types,
|
||||
nbrTypes,
|
||||
extrapolatedNorm,
|
||||
nbrNorm,
|
||||
own[facei],
|
||||
newIsFront
|
||||
);
|
||||
nChanged++;
|
||||
}
|
||||
if
|
||||
(
|
||||
isInternalFace(facei)
|
||||
&& extrapolatedNorm[nei[facei]] == -GREAT
|
||||
)
|
||||
{
|
||||
// Average nei cell, add faces to newFront
|
||||
newNorm[nei[facei]] = cellAverage
|
||||
(
|
||||
types,
|
||||
nbrTypes,
|
||||
extrapolatedNorm,
|
||||
nbrNorm,
|
||||
nei[facei],
|
||||
newIsFront
|
||||
);
|
||||
nChanged++;
|
||||
}
|
||||
}
|
||||
|
||||
reduce(nChanged, sumOp<label>());
|
||||
if (nChanged == 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
// Transfer new front
|
||||
extrapolatedNorm.transfer(newNorm);
|
||||
isFront.transfer(newIsFront);
|
||||
syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
|
||||
}
|
||||
|
||||
|
||||
forAll(norm, celli)
|
||||
{
|
||||
scalar& n = norm[celli];
|
||||
if (magSqr(n) < sqr(SMALL))
|
||||
{
|
||||
//Pout<< "For field " << m.psi().name()
|
||||
// << " for cell " << celli
|
||||
// << " at:" << cellCentres()[celli]
|
||||
// << " have norm " << n
|
||||
// << " have extrapolated norm " << extrapolatedNorm[celli]
|
||||
// << endl;
|
||||
// Override the norm
|
||||
n = extrapolatedNorm[celli];
|
||||
}
|
||||
}
|
||||
}
|
||||
return tnorm;
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::dynamicOversetFvMesh::addInterpolation
|
||||
(
|
||||
fvMatrix<Type>& m,
|
||||
const scalarField& normalisation
|
||||
) const
|
||||
{
|
||||
const cellCellStencilObject& overlap = Stencil::New(*this);
|
||||
const List<scalarList>& wghts = overlap.cellInterpolationWeights();
|
||||
const labelListList& stencil = overlap.cellStencil();
|
||||
const labelList& cellIDs = overlap.interpolationCells();
|
||||
const scalarList& factor = overlap.cellInterpolationWeight();
|
||||
const labelUList& types = overlap.cellTypes();
|
||||
|
||||
|
||||
// Force asymmetric matrix (if it wasn't already)
|
||||
scalarField& lower = m.lower();
|
||||
scalarField& upper = m.upper();
|
||||
Field<Type>& source = m.source();
|
||||
scalarField& diag = m.diag();
|
||||
|
||||
|
||||
// Get the addressing. Note that the addressing is now extended with
|
||||
// any interpolation faces.
|
||||
const lduAddressing& addr = lduAddr();
|
||||
const labelUList& upperAddr = addr.upperAddr();
|
||||
const labelUList& lowerAddr = addr.lowerAddr();
|
||||
const labelUList& ownerStartAddr = addr.ownerStartAddr();
|
||||
const labelUList& losortAddr = addr.losortAddr();
|
||||
const lduInterfacePtrsList& interfaces = allInterfaces_;
|
||||
|
||||
if (!isA<fvMeshPrimitiveLduAddressing>(addr))
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Problem : addressing is not fvMeshPrimitiveLduAddressing"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// 1. Adapt lduMatrix for additional faces and new ordering
|
||||
upper.setSize(upperAddr.size(), 0.0);
|
||||
inplaceReorder(reverseFaceMap_, upper);
|
||||
lower.setSize(lowerAddr.size(), 0.0);
|
||||
inplaceReorder(reverseFaceMap_, lower);
|
||||
|
||||
|
||||
//const label nOldInterfaces = dynamicMotionSolverFvMesh::interfaces().size();
|
||||
const label nOldInterfaces = dynamicFvMesh::interfaces().size();
|
||||
|
||||
|
||||
if (interfaces.size() > nOldInterfaces)
|
||||
{
|
||||
// Extend matrix coefficients
|
||||
m.internalCoeffs().setSize(interfaces.size());
|
||||
m.boundaryCoeffs().setSize(interfaces.size());
|
||||
|
||||
// 1b. Adapt for additional interfaces
|
||||
for
|
||||
(
|
||||
label patchi = nOldInterfaces;
|
||||
patchi < interfaces.size();
|
||||
patchi++
|
||||
)
|
||||
{
|
||||
const labelUList& fc = interfaces[patchi].faceCells();
|
||||
m.internalCoeffs().set(patchi, new Field<Type>(fc.size(), Zero));
|
||||
m.boundaryCoeffs().set(patchi, new Field<Type>(fc.size(), Zero));
|
||||
}
|
||||
|
||||
// 1c. Adapt field for additional interfaceFields (note: solver uses
|
||||
// GeometricField::scalarInterfaces() to get hold of interfaces)
|
||||
typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
|
||||
|
||||
typename GeoField::Boundary& bfld =
|
||||
const_cast<GeoField&>(m.psi()).boundaryFieldRef();
|
||||
|
||||
bfld.setSize(interfaces.size());
|
||||
|
||||
|
||||
// This gets quite interesting: we do not want to add additional
|
||||
// fvPatches (since direct correspondence to polyMesh) so instead
|
||||
// add a reference to an existing processor patch
|
||||
label addPatchi = 0;
|
||||
for (label patchi = 0; patchi < nOldInterfaces; patchi++)
|
||||
{
|
||||
if (isA<processorFvPatch>(bfld[patchi].patch()))
|
||||
{
|
||||
addPatchi = patchi;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
for
|
||||
(
|
||||
label patchi = nOldInterfaces;
|
||||
patchi < interfaces.size();
|
||||
patchi++
|
||||
)
|
||||
{
|
||||
bfld.set
|
||||
(
|
||||
patchi,
|
||||
new calculatedProcessorFvPatchField<Type>
|
||||
(
|
||||
interfaces[patchi],
|
||||
bfld[addPatchi].patch(), // dummy processorFvPatch
|
||||
m.psi()
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 2. Adapt fvMatrix level: faceFluxCorrectionPtr
|
||||
// Question: do we need to do this?
|
||||
// This seems to be set/used only by the gaussLaplacianScheme and
|
||||
// fvMatrix:correction, both of which are outside the linear solver.
|
||||
|
||||
|
||||
|
||||
// Clear out existing connections on cells to be interpolated
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
// Note: could avoid doing the zeroing of the new faces since these
|
||||
// are set to zero anyway.
|
||||
|
||||
forAll(upperAddr, facei)
|
||||
{
|
||||
if (types[upperAddr[facei]] == cellCellStencil::INTERPOLATED)
|
||||
{
|
||||
// Disconnect upper from lower
|
||||
label celli = upperAddr[facei];
|
||||
lower[facei] *= 1.0-factor[celli];
|
||||
}
|
||||
if (types[lowerAddr[facei]] == cellCellStencil::INTERPOLATED)
|
||||
{
|
||||
// Disconnect lower from upper
|
||||
label celli = lowerAddr[facei];
|
||||
upper[facei] *= 1.0-factor[celli];
|
||||
}
|
||||
}
|
||||
|
||||
for (label patchi = 0; patchi < nOldInterfaces; ++patchi)
|
||||
{
|
||||
const labelUList& fc = addr.patchAddr(patchi);
|
||||
Field<Type>& intCoeffs = m.internalCoeffs()[patchi];
|
||||
Field<Type>& bouCoeffs = m.boundaryCoeffs()[patchi];
|
||||
forAll(fc, i)
|
||||
{
|
||||
label celli = fc[i];
|
||||
{
|
||||
if (types[celli] == cellCellStencil::INTERPOLATED)
|
||||
{
|
||||
scalar f = factor[celli];
|
||||
intCoeffs[i] *= 1.0-f;
|
||||
bouCoeffs[i] *= 1.0-f;
|
||||
}
|
||||
else if (types[celli] == cellCellStencil::HOLE)
|
||||
{
|
||||
intCoeffs[i] = pTraits<Type>::zero;
|
||||
bouCoeffs[i] = pTraits<Type>::zero;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Modify matrix
|
||||
// ~~~~~~~~~~~~~
|
||||
|
||||
// Do hole cells. Note: maybe put into interpolationCells() loop above?
|
||||
forAll(types, celli)
|
||||
{
|
||||
if (types[celli] == cellCellStencil::HOLE)
|
||||
{
|
||||
label startLabel = ownerStartAddr[celli];
|
||||
label endLabel = ownerStartAddr[celli + 1];
|
||||
|
||||
for (label facei = startLabel; facei < endLabel; facei++)
|
||||
{
|
||||
upper[facei] = 0.0;
|
||||
}
|
||||
|
||||
startLabel = addr.losortStartAddr()[celli];
|
||||
endLabel = addr.losortStartAddr()[celli + 1];
|
||||
|
||||
for (label i = startLabel; i < endLabel; i++)
|
||||
{
|
||||
label facei = losortAddr[i];
|
||||
lower[facei] = 0.0;
|
||||
}
|
||||
|
||||
diag[celli] = normalisation[celli];
|
||||
source[celli] = normalisation[celli]*m.psi()[celli];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//const globalIndex globalNumbering(V().size());
|
||||
//labelList globalCellIDs(overlap.cellInterpolationMap().constructSize());
|
||||
//forAll(V(), cellI)
|
||||
//{
|
||||
// globalCellIDs[cellI] = globalNumbering.toGlobal(cellI);
|
||||
//}
|
||||
//overlap.cellInterpolationMap().distribute(globalCellIDs);
|
||||
|
||||
|
||||
forAll(cellIDs, i)
|
||||
{
|
||||
label celli = cellIDs[i];
|
||||
|
||||
const scalar f = factor[celli];
|
||||
const scalarList& w = wghts[celli];
|
||||
const labelList& nbrs = stencil[celli];
|
||||
const labelList& nbrFaces = stencilFaces_[celli];
|
||||
const labelList& nbrPatches = stencilPatches_[celli];
|
||||
|
||||
if (types[celli] == cellCellStencil::HOLE)
|
||||
{
|
||||
FatalErrorInFunction << "Found HOLE cell " << celli
|
||||
<< " at:" << C()[celli]
|
||||
<< " . Should this be in interpolationCells()????"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Create interpolation stencil
|
||||
|
||||
diag[celli] *= (1.0-f);
|
||||
source[celli] *= (1.0-f);
|
||||
|
||||
forAll(nbrs, nbri)
|
||||
{
|
||||
label patchi = nbrPatches[nbri];
|
||||
label facei = nbrFaces[nbri];
|
||||
|
||||
if (patchi == -1)
|
||||
{
|
||||
label nbrCelli = nbrs[nbri];
|
||||
|
||||
// Add the coefficients
|
||||
const scalar s = normalisation[celli]*f*w[nbri];
|
||||
|
||||
scalar& u = upper[facei];
|
||||
scalar& l = lower[facei];
|
||||
if (celli < nbrCelli)
|
||||
{
|
||||
diag[celli] += s;
|
||||
u += -s;
|
||||
}
|
||||
else
|
||||
{
|
||||
diag[celli] += s;
|
||||
l += -s;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Patch face. Store in boundaryCoeffs. Note sign change.
|
||||
//const label globalCelli = globalCellIDs[nbrs[nbri]];
|
||||
//const label proci =
|
||||
// globalNumbering.whichProcID(globalCelli);
|
||||
//const label remoteCelli =
|
||||
// globalNumbering.toLocal(proci, globalCelli);
|
||||
//
|
||||
//Pout<< "for cell:" << celli
|
||||
// << " need weight from remote slot:" << nbrs[nbri]
|
||||
// << " proc:" << proci << " remote cell:" << remoteCelli
|
||||
// << " patch:" << patchi
|
||||
// << " patchFace:" << facei
|
||||
// << " weight:" << w[nbri]
|
||||
// << endl;
|
||||
|
||||
const scalar s = normalisation[celli]*f*w[nbri];
|
||||
m.boundaryCoeffs()[patchi][facei] += pTraits<Type>::one*s;
|
||||
m.internalCoeffs()[patchi][facei] += pTraits<Type>::one*s;
|
||||
|
||||
// Note: do NOT add to diagonal - this is in the
|
||||
// internalCoeffs and gets added to the diagonal
|
||||
// inside fvMatrix::solve
|
||||
}
|
||||
}
|
||||
|
||||
//if (mag(diag[celli]) < SMALL)
|
||||
//{
|
||||
// Pout<< "for cell:" << celli
|
||||
// << " at:" << this->C()[celli]
|
||||
// << " diag:" << diag[celli] << endl;
|
||||
//
|
||||
// forAll(nbrs, nbri)
|
||||
// {
|
||||
// label patchi = nbrPatches[nbri];
|
||||
// label facei = nbrFaces[nbri];
|
||||
//
|
||||
// const label globalCelli = globalCellIDs[nbrs[nbri]];
|
||||
// const label proci =
|
||||
// globalNumbering.whichProcID(globalCelli);
|
||||
// const label remoteCelli =
|
||||
// globalNumbering.toLocal(proci, globalCelli);
|
||||
//
|
||||
// Pout<< " need weight from slot:" << nbrs[nbri]
|
||||
// << " proc:" << proci << " remote cell:"
|
||||
// << remoteCelli
|
||||
// << " patch:" << patchi
|
||||
// << " patchFace:" << facei
|
||||
// << " weight:" << w[nbri]
|
||||
// << endl;
|
||||
// }
|
||||
// Pout<< endl;
|
||||
//}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
|
||||
(
|
||||
fvMatrix<Type>& m,
|
||||
const dictionary& dict
|
||||
) const
|
||||
{
|
||||
typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
|
||||
// Check we're running with bcs that can handle implicit matrix manipulation
|
||||
typename GeoField::Boundary& bpsi =
|
||||
const_cast<GeoField&>(m.psi()).boundaryFieldRef();
|
||||
|
||||
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);
|
||||
return dynamicFvMesh::solve(m, dict);
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "dynamicOversetFvMesh::solve() :"
|
||||
<< " adjusting matrix for interpolation for field "
|
||||
<< m.psi().name() << endl;
|
||||
}
|
||||
|
||||
// Calculate stabilised diagonal as normalisation for interpolation
|
||||
const scalarField norm(normalisation(m));
|
||||
|
||||
if (debug)
|
||||
{
|
||||
volScalarField scale
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
m.psi().name() + "_normalisation",
|
||||
this->time().timeName(),
|
||||
*this,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE,
|
||||
false
|
||||
),
|
||||
*this,
|
||||
dimensionedScalar(dimless, Zero)
|
||||
);
|
||||
scale.ref().field() = norm;
|
||||
correctBoundaryConditions
|
||||
<
|
||||
volScalarField,
|
||||
oversetFvPatchField<scalar>
|
||||
>(scale.boundaryFieldRef(), false);
|
||||
scale.write();
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "dynamicOversetFvMesh::solve() :"
|
||||
<< " writing matrix normalisation for field " << m.psi().name()
|
||||
<< " to " << scale.name() << endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Switch to extended addressing (requires mesh::update() having been
|
||||
// called)
|
||||
active(true);
|
||||
|
||||
// Adapt matrix
|
||||
scalarField oldUpper(m.upper());
|
||||
scalarField oldLower(m.lower());
|
||||
FieldField<Field, Type> oldInt(m.internalCoeffs());
|
||||
FieldField<Field, Type> oldBou(m.boundaryCoeffs());
|
||||
const label oldSize = bpsi.size();
|
||||
|
||||
addInterpolation(m, norm);
|
||||
|
||||
// Swap psi values so added patches have patchNeighbourField
|
||||
correctBoundaryConditions<GeoField, calculatedProcessorFvPatchField<Type>>
|
||||
(
|
||||
bpsi,
|
||||
true
|
||||
);
|
||||
|
||||
|
||||
// Print a bit
|
||||
//write(Pout, m, lduAddr(), interfaces());
|
||||
//{
|
||||
// const fvSolution& sol = static_cast<const fvSolution&>(*this);
|
||||
// const dictionary& pDict = sol.subDict("solvers").subDict("p");
|
||||
// writeAgglomeration(GAMGAgglomeration::New(m, pDict));
|
||||
//}
|
||||
|
||||
// Use lower level solver
|
||||
//SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
|
||||
SolverPerformance<Type> s(dynamicFvMesh::solve(m, dict));
|
||||
|
||||
// Restore boundary field
|
||||
bpsi.setSize(oldSize);
|
||||
|
||||
// Restore matrix
|
||||
m.upper().transfer(oldUpper);
|
||||
m.lower().transfer(oldLower);
|
||||
m.internalCoeffs().transfer(oldInt);
|
||||
m.boundaryCoeffs().transfer(oldBou);
|
||||
|
||||
// Switch to original addressing
|
||||
active(false);
|
||||
|
||||
return s;
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
void Foam::dynamicOversetFvMesh::write
|
||||
(
|
||||
Ostream& os,
|
||||
const fvMatrix<Type>& m,
|
||||
const lduAddressing& addr,
|
||||
const lduInterfacePtrsList& interfaces
|
||||
) const
|
||||
{
|
||||
os << "** Matrix **" << endl;
|
||||
const labelUList& upperAddr = addr.upperAddr();
|
||||
const labelUList& lowerAddr = addr.lowerAddr();
|
||||
const labelUList& ownerStartAddr = addr.ownerStartAddr();
|
||||
const labelUList& losortAddr = addr.losortAddr();
|
||||
|
||||
const scalarField& lower = m.lower();
|
||||
const scalarField& upper = m.upper();
|
||||
const Field<Type>& source = m.source();
|
||||
const scalarField& diag = m.diag();
|
||||
|
||||
|
||||
// Invert patch addressing
|
||||
labelListList cellToPatch(addr.size());
|
||||
labelListList cellToPatchFace(addr.size());
|
||||
{
|
||||
forAll(interfaces, patchi)
|
||||
{
|
||||
if (interfaces.set(patchi))
|
||||
{
|
||||
const labelUList& fc = interfaces[patchi].faceCells();
|
||||
|
||||
forAll(fc, i)
|
||||
{
|
||||
cellToPatch[fc[i]].append(patchi);
|
||||
cellToPatchFace[fc[i]].append(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
forAll(source, celli)
|
||||
{
|
||||
os << "cell:" << celli << " diag:" << diag[celli]
|
||||
<< " source:" << source[celli] << endl;
|
||||
|
||||
label startLabel = ownerStartAddr[celli];
|
||||
label endLabel = ownerStartAddr[celli + 1];
|
||||
|
||||
for (label facei = startLabel; facei < endLabel; facei++)
|
||||
{
|
||||
os << " face:" << facei
|
||||
<< " nbr:" << upperAddr[facei] << " coeff:" << upper[facei]
|
||||
<< endl;
|
||||
}
|
||||
|
||||
startLabel = addr.losortStartAddr()[celli];
|
||||
endLabel = addr.losortStartAddr()[celli + 1];
|
||||
|
||||
for (label i = startLabel; i < endLabel; i++)
|
||||
{
|
||||
label facei = losortAddr[i];
|
||||
os << " face:" << facei
|
||||
<< " nbr:" << lowerAddr[facei] << " coeff:" << lower[facei]
|
||||
<< endl;
|
||||
}
|
||||
|
||||
forAll(cellToPatch[celli], i)
|
||||
{
|
||||
label patchi = cellToPatch[celli][i];
|
||||
label patchFacei = cellToPatchFace[celli][i];
|
||||
|
||||
os << " patch:" << patchi
|
||||
<< " patchface:" << patchFacei
|
||||
<< " intcoeff:" << m.internalCoeffs()[patchi][patchFacei]
|
||||
<< " boucoeff:" << m.boundaryCoeffs()[patchi][patchFacei]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
forAll(m.internalCoeffs(), patchi)
|
||||
{
|
||||
if (m.internalCoeffs().set(patchi))
|
||||
{
|
||||
const labelUList& fc = addr.patchAddr(patchi);
|
||||
|
||||
os << "patch:" << patchi
|
||||
//<< " type:" << interfaces[patchi].type()
|
||||
<< " size:" << fc.size() << endl;
|
||||
if
|
||||
(
|
||||
interfaces.set(patchi)
|
||||
&& isA<processorLduInterface>(interfaces[patchi])
|
||||
)
|
||||
{
|
||||
const processorLduInterface& ppp =
|
||||
refCast<const processorLduInterface>(interfaces[patchi]);
|
||||
os << "(processor with my:" << ppp.myProcNo()
|
||||
<< " nbr:" << ppp.neighbProcNo()
|
||||
<< ")" << endl;
|
||||
}
|
||||
|
||||
forAll(fc, i)
|
||||
{
|
||||
os << " cell:" << fc[i]
|
||||
<< " int:" << m.internalCoeffs()[patchi][i]
|
||||
<< " boun:" << m.boundaryCoeffs()[patchi][i]
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
lduInterfaceFieldPtrsList interfaceFields =
|
||||
m.psi().boundaryField().scalarInterfaces();
|
||||
forAll(interfaceFields, inti)
|
||||
{
|
||||
if (interfaceFields.set(inti))
|
||||
{
|
||||
os << "interface:" << inti
|
||||
<< " if type:" << interfaceFields[inti].interface().type()
|
||||
<< " fld type:" << interfaceFields[inti].type() << endl;
|
||||
}
|
||||
}
|
||||
|
||||
os << "** End of Matrix **" << endl;
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
|
||||
{
|
||||
typename GeoField::Boundary& bfld = fld.boundaryFieldRef();
|
||||
|
||||
const label nReq = Pstream::nRequests();
|
||||
|
||||
forAll(bfld, patchi)
|
||||
{
|
||||
if (bfld[patchi].coupled())
|
||||
{
|
||||
//Pout<< "initEval of " << bfld[patchi].patch().name() << endl;
|
||||
bfld[patchi].initEvaluate(Pstream::defaultCommsType);
|
||||
}
|
||||
}
|
||||
|
||||
// Block for any outstanding requests
|
||||
if (Pstream::parRun())
|
||||
{
|
||||
Pstream::waitRequests(nReq);
|
||||
}
|
||||
|
||||
forAll(bfld, patchi)
|
||||
{
|
||||
if (bfld[patchi].coupled())
|
||||
{
|
||||
//Pout<< "eval of " << bfld[patchi].patch().name() << endl;
|
||||
bfld[patchi].evaluate(Pstream::defaultCommsType);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class GeoField>
|
||||
void Foam::dynamicOversetFvMesh::checkCoupledBC(const GeoField& fld)
|
||||
{
|
||||
Pout<< "** starting checking of " << fld.name() << endl;
|
||||
|
||||
GeoField fldCorr(fld.name()+"_correct", fld);
|
||||
correctCoupledBoundaryConditions(fldCorr);
|
||||
|
||||
const typename GeoField::Boundary& bfld = fld.boundaryField();
|
||||
const typename GeoField::Boundary& bfldCorr = fldCorr.boundaryField();
|
||||
|
||||
forAll(bfld, patchi)
|
||||
{
|
||||
const typename GeoField::Patch& pfld = bfld[patchi];
|
||||
const typename GeoField::Patch& pfldCorr = bfldCorr[patchi];
|
||||
|
||||
Pout<< "Patch:" << pfld.patch().name() << endl;
|
||||
|
||||
forAll(pfld, i)
|
||||
{
|
||||
if (pfld[i] != pfldCorr[i])
|
||||
{
|
||||
Pout<< " " << i << " orig:" << pfld[i]
|
||||
<< " corrected:" << pfldCorr[i] << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
Pout<< "** end of " << fld.name() << endl;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user