Commit 880adc19 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

port mpfa to new structure for network grids

parent ae380a02
......@@ -60,7 +60,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FluxVarsCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using ElementFluxVarsCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
......@@ -76,9 +76,11 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const unsigned int phaseIdx,
const FluxVarsCache& fluxVarsCache)
const ElementFluxVarsCache& elemFluxVarsCache)
{
const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil(phaseIdx);
const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions(phaseIdx);
const auto& tij = fluxVarsCache.advectionTij(phaseIdx);
......
......@@ -75,6 +75,13 @@ public:
return globalFvGeometry().scvf(scvfIdx);
}
//! Get an element sub control volume face with a global scvf index
//! We separate element and neighbor scvfs to speed up mapping
const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvIdx = 0) const
{
return globalFvGeometry().flipScvf(scvfIdx, outsideScvIdx);
}
//! iterator range for sub control volumes. Iterates over
//! all scvs of the bound element (not including neighbor scvs)
//! This is a free function found by means of ADL
......@@ -187,6 +194,29 @@ public:
return neighborScvfs_[findLocalIndex(scvfIdx, neighborScvfIndices_)];
}
//! Get an element sub control volume face with a global scvf index
//! We separate element and neighbor scvfs to speed up mapping
const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvIdx = 0) const
{
auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
if (it != scvfIndices_.end())
{
assert(flipScvfIndices_[std::distance(scvfIndices_.begin(), it)][outsideScvIdx] != -1);
return neighborScvfs_[ flipScvfIndices_[std::distance(scvfIndices_.begin(), it)][outsideScvIdx] ];
}
else
{
const auto neighborScvfIdxLocal = findLocalIndex(scvfIdx, neighborScvfIndices_);
auto&& scvf = neighborScvfs_[neighborScvfIdxLocal];
assert(neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] != -1);
if (scvf.outsideScvIdx(outsideScvIdx) == scvIndices_[0])
return scvfs_[ neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] ];
else
return neighborScvfs_[ neighborFlipScvfIndices_[neighborScvfIdxLocal][outsideScvIdx] ];
}
}
//! iterator range for sub control volumes. Iterates over
//! all scvs of the bound element (not including neighbor scvs)
//! This is a free function found by means of ADL
......@@ -260,6 +290,10 @@ public:
neighborScvIndices_.shrink_to_fit();
neighborScvfIndices_.shrink_to_fit();
neighborScvfs_.shrink_to_fit();
// set up the scvf flip indices for network grids
if (dim < dimWorld)
makeFlipIndexSet();
}
//! Binding of an element preparing the geometries only inside the element
......@@ -448,6 +482,97 @@ private:
}
}
void makeFlipIndexSet()
{
const auto numInsideScvfs = scvfIndices_.size();
const auto numNeighborScvfs = neighborScvfIndices_.size();
// first, handle the interior scvfs (flip scvf will be in outside scvfs)
flipScvfIndices_.resize(numInsideScvfs);
for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
{
auto&& scvf = scvfs_[insideFace];
if (scvf.boundary())
continue;
const auto numOutsideScvs = scvf.numOutsideScvs();
const auto vIdxGlobal = scvf.vertexIndex();
const auto insideScvIdx = scvf.insideScvIdx();
flipScvfIndices_[insideFace].resize(numOutsideScvs, -1);
for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
{
const auto outsideScvIdx = scvf.outsideScvIdx(nIdx);
for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
{
auto&& outsideScvf = neighborScvfs_[outsideFace];
// check if outside face is the flip face
if (outsideScvf.insideScvIdx() == outsideScvIdx &&
outsideScvf.vertexIndex() == vIdxGlobal &&
MpfaHelper::contains(outsideScvf.outsideScvIndices(), insideScvIdx))
{
flipScvfIndices_[insideFace][nIdx] = outsideFace;
// there is always only one flip face in an outside element
break;
}
}
}
}
// Now we look for the flip indices of the outside faces
neighborFlipScvfIndices_.resize(numNeighborScvfs);
for (unsigned int outsideFace = 0; outsideFace < numNeighborScvfs; ++outsideFace)
{
auto&& scvf = neighborScvfs_[outsideFace];
if (scvf.boundary())
continue;
const auto numOutsideScvs = scvf.numOutsideScvs();
const auto vIdxGlobal = scvf.vertexIndex();
const auto insideScvIdx = scvf.insideScvIdx();
neighborFlipScvfIndices_[outsideFace].resize(numOutsideScvs, -1);
for (unsigned int nIdx = 0; nIdx < numOutsideScvs; ++nIdx)
{
const auto neighborScvIdx = scvf.outsideScvIdx(nIdx);
// if the neighbor scv idx is the index of the bound element,
// then the flip scvf will be within the inside scvfs
if (neighborScvIdx == scvIndices_[0])
{
for (unsigned int insideFace = 0; insideFace < numInsideScvfs; ++insideFace)
{
auto&& insideScvf = scvfs_[insideFace];
// check if the face is the flip face
if (insideScvf.vertexIndex() == vIdxGlobal &&
MpfaHelper::contains(insideScvf.outsideScvIndices(), insideScvIdx))
{
neighborFlipScvfIndices_[outsideFace][nIdx] = insideFace;
// there is always only one flip face in an outside element
break;
}
}
}
else
{
for (unsigned int outsideFace2 = 0; outsideFace2 < numNeighborScvfs; ++outsideFace2)
{
const auto& outsideScvf = neighborScvfs_[outsideFace2];
// check if outside face is the flip face
if (outsideScvf.insideScvIdx() == neighborScvIdx &&
outsideScvf.vertexIndex() == vIdxGlobal &&
MpfaHelper::contains(outsideScvf.outsideScvIndices(), insideScvIdx))
{
neighborFlipScvfIndices_[outsideFace][nIdx] = outsideFace2;
// there is always only one flip face in an outside element
break;
}
}
}
}
}
}
const IndexType findLocalIndex(const IndexType idx,
const std::vector<IndexType>& indices) const
{
......@@ -484,6 +609,10 @@ private:
std::vector<IndexType> neighborScvfIndices_;
std::vector<SubControlVolume> neighborScvs_;
std::vector<SubControlVolumeFace> neighborScvfs_;
// flip index sets
std::vector< std::vector<IndexType> > flipScvfIndices_;
std::vector< std::vector<IndexType> > neighborFlipScvfIndices_;
};
} // end namespace
......
......@@ -285,6 +285,39 @@ public:
// in parallel problems we might have reserved more scvfs than we actually use
scvfs_.shrink_to_fit();
// Make the flip index set for network and surface grids
if (dim < dimWorld)
{
flipScvfIndices_.resize(scvfs_.size());
for (auto&& scvf : scvfs_)
{
if (scvf.boundary())
continue;
const auto numOutsideScvs = scvf.numOutsideScvs();
const auto vIdxGlobal = scvf.vertexIndex();
const auto insideScvIdx = scvf.insideScvIdx();
flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
for (unsigned int i = 0; i < numOutsideScvs; ++i)
{
const auto outsideScvIdx = scvf.outsideScvIdx(i);
for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
{
const auto& outsideScvf = this->scvf(outsideScvfIndex);
if (outsideScvf.vertexIndex() == vIdxGlobal &&
MpfaHelper::contains(outsideScvf.outsideScvIndices(), insideScvIdx))
{
flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
// there is always only one flip face in an outside element
break;
}
}
}
}
}
// Initialize the interaction volume seeds
globalInteractionVolumeSeeds_.update(problem, boundaryVertices_);
}
......@@ -305,6 +338,11 @@ public:
const SubControlVolumeFace& scvf(IndexType scvfIdx) const
{ return scvfs_[scvfIdx]; }
const SubControlVolumeFace& flipScvf(IndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
{
return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]];
}
//! Get the sub control volume face indices of an scv by global index
const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
{ return scvfIndicesOfScv_[scvIdx]; }
......@@ -331,6 +369,8 @@ private:
std::vector<bool> boundaryVertices_;
std::vector<bool> branchingVertices_;
IndexType numBoundaryScvf_;
// needed for embedded surface and network grids (dim < dimWorld)
std::vector<std::vector<IndexType>> flipScvfIndices_;
// the global interaction volume seeds
GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_;
......
......@@ -117,12 +117,16 @@ public:
{ return insideScvIdx_; }
//! index of the outside sub control volume for spatial param evaluation
//! returns in undefined behaviour is boundary is true
//! is not uniquely defined for grids where dim < dimWorld and shouldn't be used in that case.
IndexType outsideScvIdx() const
//! returns in undefined behaviour if boundary is true or index exceeds numOutsideScvs
IndexType outsideScvIdx(int i = 0) const
{
assert(outsideScvIndices_.size() == 1 && "outside scv index not uniquely defined");
return outsideScvIndices_[0];
return outsideScvIndices_[i];
}
//! The number of outside scvs connection via this scv face
std::size_t numOutsideScvs() const
{
return outsideScvIndices_.size();
}
//! returns the outside scv indices (can be more than one index for dim < dimWorld)
......
......@@ -78,10 +78,10 @@ protected:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf,
const ElementBoundaryTypes& bcTypes,
const FluxVariablesCache& fluxVarsCache)
const ElementFluxVariablesCache& elemFluxVarsCache)
{
if (!scvf.boundary() || !useTpfaBoundary)
return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
else
return PrimaryVariables(0.0);
......@@ -91,7 +91,7 @@ protected:
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache)
const ElementFluxVariablesCache& elemFluxVarsCache)
{
ElementBoundaryTypes bcTypes;
bcTypes.update(this->problem(), element, fvGeometry);
......@@ -100,11 +100,11 @@ protected:
this->residual_ = 0;
if (!scvf.boundary() || !useTpfaBoundary)
return this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, bcTypes, fluxVarsCache);
return this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache);
else
{
auto bcTypes = this->problem().boundaryTypes(element, scvf);
return this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, fluxVarsCache);
return this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache);
}
}
......@@ -122,7 +122,7 @@ protected:
if (scvf.boundary())
{
auto bcTypes = this->problem().boundaryTypes(element, scvf);
this->residual_[0] += this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache[scvf]);
this->residual_[0] += this->asImp_().evalBoundaryFluxes_(element, fvGeometry, elemVolVars, scvf, bcTypes, elemFluxVarsCache);
}
}
}
......@@ -136,11 +136,11 @@ protected:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf,
const BoundaryTypes& bcTypes,
const FluxVariablesCache& fluxVarsCache)
const ElementFluxVariablesCache& elemFluxVarsCache)
{
// return the boundary fluxes
if (bcTypes.hasOnlyDirichlet())
return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
else if (bcTypes.hasOnlyNeumann())
return this->asImp_().evalNeumannSegment_(element, fvGeometry, elemVolVars, scvf, bcTypes);
else if (bcTypes.hasOutflow())
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment