Commit db8cf19d authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[facet][mortar] first draft of coupling manager

parent d72e8293
This diff is collapsed.
......@@ -81,87 +81,7 @@ public:
if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) )
DUNE_THROW(Dune::NotImplemented, "Xi != 1.0 cannot be used with the Box-Facet-Coupling scheme");
// get some references for convenience
const auto& fluxVarCache = elemFluxVarCache[scvf];
const auto& shapeValues = fluxVarCache.shapeValues();
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
// evaluate user-defined interior boundary types
const auto bcTypes = problem.interiorBoundaryTypes(element, scvf);
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
// on interior Neumann boundaries, evaluate the flux using the facet permeability
if (bcTypes.hasOnlyNeumann())
{
// interpolate pressure/density to scvf integration point
Scalar p = 0.0;
Scalar rho = 0.0;
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
p += volVars.pressure(phaseIdx)*shapeValues[scv.indexInElement()][0];
rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
}
// compute tpfa flux from integration point to facet centerline
const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
using std::sqrt;
// If this is a surface grid, use the square root of the facet extrusion factor
// as an approximate average distance from scvf ip to facet center
using std::sqrt;
const auto a = facetVolVars.extrusionFactor();
auto gradP = scvf.unitOuterNormal();
gradP *= dim == dimWorld ? 0.5*a : 0.5*sqrt(a);
gradP /= gradP.two_norm2();
gradP *= (facetVolVars.pressure(phaseIdx) - p);
if (enableGravity)
gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center()));
// apply facet permeability and return the flux
return -1.0*scvf.area()
*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), gradP);
}
// on interior Dirichlet boundaries use the facet pressure and evaluate flux
else if (bcTypes.hasOnlyDirichlet())
{
// create vector with nodal pressures
std::vector<Scalar> pressures(element.subEntities(dim));
for (const auto& scv : scvs(fvGeometry))
pressures[scv.localDofIndex()] = elemVolVars[scv].pressure(phaseIdx);
// substitute with facet pressures for those scvs touching this facet
for (const auto& scvfJ : scvfs(fvGeometry))
if (scvfJ.interiorBoundary() && scvfJ.facetIndexInElement() == scvf.facetIndexInElement())
pressures[ fvGeometry.scv(scvfJ.insideScvIdx()).localDofIndex() ]
= problem.couplingManager().getLowDimVolVars(element, scvfJ).pressure(phaseIdx);
// evaluate gradP - rho*g at integration point
Scalar rho(0.0);
Dune::FieldVector<Scalar, dimWorld> gradP(0.0);
for (const auto& scv : scvs(fvGeometry))
{
rho += elemVolVars[scv].density(phaseIdx)*shapeValues[scv.indexInElement()][0];
gradP.axpy(pressures[scv.localDofIndex()], fluxVarCache.gradN(scv.indexInElement()));
}
if (enableGravity)
gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center()));
// apply matrix permeability and return the flux
return -1.0*scvf.area()
*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), gradP);
}
// mixed boundary types are not supported
else
DUNE_THROW(Dune::NotImplemented, "Mixed boundary types are not supported");
return problem.couplingManager().getMortarFlux(element, scvf);
}
// compute transmissibilities ti for analytical jacobians
......
......@@ -53,30 +53,31 @@ public:
const auto& scvf = fluxVars.scvFace();
const auto& insideScv = fluxVars.fvGeometry().scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
return flux*upwindTerm(insideVolVars);
// check if this is an interior boundary
if (scvf.interiorBoundary())
{
const auto& cm = fluxVars.problem().couplingManager();
const auto& outsideVolVars = cm.getLowDimVolVars(fluxVars.element(), scvf);
if (std::signbit(flux)) // if sign of flux is negative
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
else
return flux*(upwindWeight*upwindTerm(insideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
}
else
{
const auto& outsideScv = fluxVars.fvGeometry().scv(scvf.outsideScvIdx());
const auto& outsideVolVars = elemVolVars[outsideScv];
if (std::signbit(flux)) // if sign of flux is negative
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
else
return flux*(upwindWeight*upwindTerm(insideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
}
// // check if this is an interior boundary
// if (scvf.interiorBoundary())
// {
// const auto& cm = fluxVars.problem().couplingManager();
// const auto& outsideVolVars = cm.getLowDimVolVars(fluxVars.element(), scvf);
// if (std::signbit(flux)) // if sign of flux is negative
// return flux*(upwindWeight*upwindTerm(outsideVolVars)
// + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
// else
// return flux*(upwindWeight*upwindTerm(insideVolVars)
// + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
// }
// else
// {
// const auto& outsideScv = fluxVars.fvGeometry().scv(scvf.outsideScvIdx());
// const auto& outsideVolVars = elemVolVars[outsideScv];
// if (std::signbit(flux)) // if sign of flux is negative
// return flux*(upwindWeight*upwindTerm(outsideVolVars)
// + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
// else
// return flux*(upwindWeight*upwindTerm(insideVolVars)
// + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
// }
}
};
......
......@@ -32,6 +32,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/assembly/felocalresidual.hh>
#include <dumux/discretization/fem/ipdata.hh>
namespace Dumux {
......
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