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

[mpfao] remove interior boundary support

The support for interior boundaries will be handled explicitly in the
facet coupling implementation in the multidimension module. Its support
in the standard mpfa implementation causes an overhead which we don't want
here since interior boundaries do not have to be available for standard models.
parent 65b9a722
......@@ -114,10 +114,6 @@ public:
bool scvfTouchesBoundary(const SubControlVolumeFace& scvf) const
{ return boundaryVertices_[scvf.vertexIndex()]; }
//! returns true if the vertex touches the domain boundary
bool isDomainBoundaryVertex(const IndexType vIdxGlobal) const
{ return outerBoundaryVertices_[vIdxGlobal]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update(const Problem& problem)
{
......@@ -130,7 +126,6 @@ public:
scvfs_.reserve(getNumScvf_());
scvfIndicesOfScv_.resize(numScvs);
boundaryVertices_.resize(gridView_.size(dim), false);
outerBoundaryVertices_.resize(gridView_.size(dim), false);
elementMap_.resize(numScvs);
// the quadrature point to be used on the scvf
......@@ -190,7 +185,7 @@ public:
const auto vIdxGlobal = problem.vertexMapper().subIndex(element, vIdxLocal, dim);
if (boundary)
outerBoundaryVertices_[vIdxGlobal] = true;
boundaryVertices_[vIdxGlobal] = true;
// make the scv face
scvfIndexSet.push_back(scvfIdx);
......@@ -214,7 +209,7 @@ public:
}
// Initialize the interaction volume seeds, this will also initialize the vector of boundary vertices
globalInteractionVolumeSeeds_.update(problem, boundaryVertices_);
globalInteractionVolumeSeeds_.update(problem);
}
/*!
......@@ -272,7 +267,6 @@ public:
std::vector<SubControlVolumeFace> scvfs_;
std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<bool> boundaryVertices_;
std::vector<bool> outerBoundaryVertices_;
IndexType numBoundaryScvf_;
GlobalInteractionVolumeSeeds globalInteractionVolumeSeeds_;
};
......
......@@ -56,8 +56,7 @@ public:
CCMpfaGlobalInteractionVolumeSeeds(const GridView gridView) : gridView_(gridView) {}
// initializes the interaction volumes or the seeds
template<typename BoolVector>
void update(const Problem& problem, BoolVector& vertexTouchesBoundary)
void update(const Problem& problem)
{
problemPtr_ = &problem;
seeds_.clear();
......@@ -68,7 +67,7 @@ public:
scvfIndexMap_.resize(numScvf, -1);
// detect and handle the boundary first
initializeBoundarySeeds_(vertexTouchesBoundary);
initializeBoundarySeeds_();
initializeInteriorSeeds_();
}
......@@ -81,20 +80,8 @@ public:
return boundarySeeds_[scvfIndexMap_[scvf.index()]][eqIdx];
}
//! returns whether or not an scvf is on an interior or outer boundary
bool isScvfOnInteriorBoundary(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf)
{
for (LocalIndexType eqIdx = 0; eqIdx < numEq; ++eqIdx)
if (Helper::getMpfaFaceType(problem_(), element, scvf, eqIdx) != MpfaFaceTypes::interior)
return true;
return false;
}
private:
template<typename BoolVector>
void initializeBoundarySeeds_(BoolVector& vertexTouchesBoundary)
void initializeBoundarySeeds_()
{
auto numBoundaryScvf = problem_().model().globalFvGeometry().numBoundaryScvf();
boundarySeeds_.reserve(numBoundaryScvf);
......@@ -106,17 +93,10 @@ private:
fvGeometry.bind(element);
for (const auto& scvf : scvfs(fvGeometry))
{
// skip the rest if we already handled this face
if (scvfIndexMap_[scvf.index()] != -1)
// skip the rest if we already handled this face or if face is not on boundary
if (scvfIndexMap_[scvf.index()] != -1 || !scvf.boundary())
continue;
// also, skip the rest if this face is not on a boundary
if (!scvf.boundary() && !isScvfOnInteriorBoundary(problem_(), element, scvf))
continue;
// the vertex connected to this scvf touches an interior or outer boundary
vertexTouchesBoundary[scvf.vertexIndex()] = true;
// container to store the interaction volume seeds
std::vector<BoundaryInteractionVolumeSeed> seedVector;
seedVector.reserve(numEq);
......
......@@ -97,18 +97,18 @@ private:
const SubControlVolumeFace& scvf,
const LocalIndexType eqIdx)
{
// The vertex index around which we construct the interaction volume
auto vIdxGlobal = scvf.vertexIndex();
bool onBoundary = problem.model().globalFvGeometry().isDomainBoundaryVertex(vIdxGlobal);
// Check whether or not we are touching the boundary here
bool onBoundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf);
// Get the two scv faces in the first scv
auto scvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, element, fvGeometry);
auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
// The global index of the first scv of the interaction region
auto scvIdx0 = scvf.insideScvIdx();
// rotate counter clockwise and create the entities
performRotation_(problem, scvfVector, scvSeeds, scvfSeeds, scvIdx0, eqIdx);
if (onBoundary)
{
// the local scvf index of the second local face of the first local scv
......@@ -185,13 +185,6 @@ private:
GlobalIndexSet({globalScvfIdx, scvfVector[1]->index()}),
faceType));
// create duplicate face if it is an interior boundary,
if (faceType != MpfaFaceTypes::interior)
scvfSeeds.emplace_back(ScvfSeed(curScvf,
LocalIndexSet({0, insideLocalScvIdx}),
GlobalIndexSet({scvfVector[1]->index(), globalScvfIdx}),
faceType));
// rotation loop is finished
finished = true; return;
}
......@@ -217,17 +210,6 @@ private:
faceType));
localScvfIdx++;
// create duplicate face if it is an interior boundary,
if (faceType != MpfaFaceTypes::interior)
{
// face from "outside" to "inside", index sets change
scvfSeeds.emplace_back(ScvfSeed(commonScvf,
LocalIndexSet({outsideLocalScvIdx, insideLocalScvIdx}),
GlobalIndexSet({commonGlobalScvfIdx, globalScvfIdx}),
faceType));
localScvfIdx++;
}
// create index set storing the two local scvf indices
LocalIndexSet localScvfs(2);
localScvfs[commonFaceCoordIdx] = localScvfIdx-1;
......@@ -383,8 +365,7 @@ private:
{
if (scvSeed.globalIndex() == outsideGlobalScvIdx)
{
outsideExists = true;
break;
outsideExists = true; break;
}
// keep track of local index
outsideLocalScvIdx++;
......@@ -426,22 +407,8 @@ private:
// boundary scvf seeds have no outside scvf
if (!scvfSeed.boundary() && scvfSeed.outsideGlobalScvfIndex() == globalScvfIndex)
{
auto faceType = scvfSeed.faceType();
// on interior boundaries we have to create a new face
if (faceType != MpfaFaceTypes::interior)
{
// set the local scvfIndex of the face that is about to created
actualScvSeed.setLocalScvfIndex(coordDir, scvfSeeds.size());
// some data on the face
LocalIndexSet scvIndicesLocal({actualLocalScvIdx, scvfSeed.insideLocalScvIndex()});
GlobalIndexSet scvfIndicesGlobal({globalScvfIndex, scvfSeed.insideGlobalScvfIndex()});
scvfSeeds.emplace_back(*scvfVector[coordDir], std::move(scvIndicesLocal), std::move(scvfIndicesGlobal), faceType);
}
// pass local scvf index to local scv
else
actualScvSeed.setLocalScvfIndex(coordDir, localScvfIdx);
actualScvSeed.setLocalScvfIndex(coordDir, localScvfIdx);
// we found the corresponding face
found = true; break;
......
......@@ -141,8 +141,7 @@ public:
: problemRef_(problem),
fvGeometryRef_(fvGeometry),
elemVolVarsRef_(elemVolVars),
onBoundary_(seed.onBoundary()),
hasInteriorBoundary_(false)
onBoundary_(seed.onBoundary())
{
// create local sub control entities from the seed
createLocalEntities_(seed);
......@@ -167,21 +166,14 @@ public:
auto faceType = localScvf.faceType();
// eventually add vol var index and corresponding position
if (faceType == MpfaFaceTypes::dirichlet || faceType == MpfaFaceTypes::interiorDirichlet)
if (faceType == MpfaFaceTypes::dirichlet)
{
volVarsStencil_.push_back(localScvf.outsideGlobalScvIndex());
volVarsPositions_.push_back(localScvf.ip());
dirichletFaceIndexSet_.push_back(localScvfIdx++);
}
// fill local scvf type index sets accordingly
if (faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet)
fluxFaceIndexSet_.push_back(localScvfIdx++);
else
dirichletFaceIndexSet_.push_back(localScvfIdx++);
// if face is on an interior boundary, save this information
if (!hasInteriorBoundary_ && (faceType == MpfaFaceTypes::interiorDirichlet || faceType == MpfaFaceTypes::interiorNeumann))
hasInteriorBoundary_ = true;
fluxFaceIndexSet_.push_back(localScvfIdx++);
}
neumannFluxes_ = DynamicVector(fluxFaceIndexSet_.size(), 0.0);
......@@ -227,7 +219,7 @@ public:
{
const auto& localScvf = localScvf_(localFluxFaceIdx);
const auto faceType = localScvf.faceType();
if (faceType == MpfaFaceTypes::neumann || faceType == MpfaFaceTypes::interiorNeumann)
if (faceType == MpfaFaceTypes::neumann)
{
// boundary neumann fluxes stay zero when tpfa boundary handling is on
if (GET_PROP_VALUE(TypeTag, UseTpfaBoundary) && faceType == MpfaFaceTypes::neumann)
......@@ -238,13 +230,10 @@ public:
auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx];
neumannFlux *= globalScvf.area();
// if we are not on an interior neumann boundary, we have to recover -k*gradh
if (faceType == MpfaFaceTypes::neumann)
{
const auto& insideScv = fvGeometry_().scv(globalScvf.insideScvIdx());
const auto& volVars = elemVolVars_()[insideScv];
neumannFlux /= upwindFactor(volVars);
}
// recover -k*gradh
const auto& insideScv = fvGeometry_().scv(globalScvf.insideScvIdx());
const auto& volVars = elemVolVars_()[insideScv];
neumannFlux /= upwindFactor(volVars);
neumannFluxes_[fluxFaceIdx] = neumannFlux;
}
......@@ -256,9 +245,6 @@ public:
bool onBoundary() const
{ return onBoundary_; }
bool hasInteriorBoundary() const
{ return hasInteriorBoundary_; }
const GlobalIdxSet& stencil() const
{ return stencil_; }
......@@ -280,13 +266,6 @@ public:
globalScvfs.push_back(localScvf.outsideGlobalScvfIndex());
}
// if interior boundaries are present we have to make the entries unique
if (hasInteriorBoundary())
{
std::sort(globalScvfs.begin(), globalScvfs.end());
globalScvfs.erase(std::unique(globalScvfs.begin(), globalScvfs.end()), globalScvfs.end());
}
return globalScvfs;
}
......@@ -377,7 +356,6 @@ private:
DynamicMatrix& C,
DynamicMatrix& D)
{
const Scalar xi = GET_PROP_VALUE(TypeTag, Xi);
const std::size_t numLocalScvs = localScvs_.size();
// loop over the local faces
......@@ -385,7 +363,7 @@ private:
for (const auto& localScvf : localScvfs_)
{
auto faceType = localScvf.faceType();
bool hasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
bool hasUnknown = faceType != MpfaFaceTypes::dirichlet;
LocalIdxType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1;
// get diffusion tensor in "positive" sub volume
......@@ -405,7 +383,7 @@ private:
auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
auto curFaceType = curLocalScvf.faceType();
bool curFaceHasUnknown = curFaceType != MpfaFaceTypes::dirichlet && curFaceType != MpfaFaceTypes::interiorDirichlet;
bool curFaceHasUnknown = curFaceType != MpfaFaceTypes::dirichlet;
// First, add the entries associated with face pressures (unkown or dirichlet)
if (curFaceHasUnknown)
......@@ -415,24 +393,14 @@ private:
C[rowIdx][curIdxInFluxFaces] += posWijk[localDir];
if (hasUnknown)
{
if (faceType == MpfaFaceTypes::interiorNeumann)
{
A[idxInFluxFaces][curIdxInFluxFaces] += xi*posWijk[localDir];
// If facet coupling active, eventually add entry coming from the lowdim domain
if (GET_PROP_VALUE(TypeTag, FacetCoupling) && curIdxInFluxFaces == idxInFluxFaces)
DUNE_THROW(Dune::NotImplemented, "Facet coupling");
}
else
A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
}
A[idxInFluxFaces][curIdxInFluxFaces] += posWijk[localDir];
}
else
{
// the current face is a Dirichlet face and creates entries in D & eventually B
auto curIdxInDiriFaces = this->findLocalIndex(dirichletScvfIndexSet_(), curLocalScvfIdx);
D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
D[rowIdx][numLocalScvs + curIdxInDiriFaces] += posWijk[localDir];
if (hasUnknown)
B[idxInFluxFaces][numLocalScvs + curIdxInDiriFaces] -= posWijk[localDir];
}
......@@ -441,16 +409,11 @@ private:
D[rowIdx][posLocalScvIdx] -= posWijk[localDir];
if (hasUnknown)
{
if (faceType == MpfaFaceTypes::interiorNeumann)
B[idxInFluxFaces][posLocalScvIdx] += xi*posWijk[localDir];
else
B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
}
B[idxInFluxFaces][posLocalScvIdx] += posWijk[localDir];
}
// If not on a boundary or interior dirichlet face, add entries for the "negative" scv
if (faceType == MpfaFaceTypes::interior || faceType == MpfaFaceTypes::interiorNeumann)
if (faceType == MpfaFaceTypes::interior)
{
const auto negLocalScvIdx = localScvf.outsideLocalScvIndex();
const auto& negLocalScv = localScv_(negLocalScvIdx);
......@@ -474,11 +437,7 @@ private:
{
// we need the index of the current local scvf in the flux face indices
auto curIdxInFluxFaces = this->findLocalIndex(fluxScvfIndexSet_(), curLocalScvfIdx);
if (faceType == MpfaFaceTypes::interiorNeumann)
A[idxInFluxFaces][curIdxInFluxFaces] -= (1.0-xi)*negWijk[localDir];
else
A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
A[idxInFluxFaces][curIdxInFluxFaces] -= negWijk[localDir];
}
else
{
......@@ -488,10 +447,7 @@ private:
}
// add entries to matrix B
if (faceType == MpfaFaceTypes::interiorNeumann)
B[idxInFluxFaces][negLocalScvIdx] -= (1.0-xi)*negWijk[localDir];
else
B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
B[idxInFluxFaces][negLocalScvIdx] -= negWijk[localDir];
}
}
// go to the next face
......@@ -589,7 +545,6 @@ private:
const ElementVolumeVariables& elemVolVarsRef_;
bool onBoundary_;
bool hasInteriorBoundary_;
LocalIdxType eqIdx_;
std::vector<Element> localElements_;
......
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