Commit 79f6307c authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[discretization] throw if evalGradients is called with non-uniform phase states

Evaluating gradients inside an element is only correct if all
element vertices have the same phase state, i.e. the same set of
primary variables. Determine if all vertex phase states are the
same. If not, throw. Enable the possibility to ignore the phase
states and enforce the gradient evaluation.
parent 02c18be7
......@@ -25,10 +25,15 @@
#define DUMUX_DISCRETIZATION_EVAL_GRADIENTS_HH
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/typetraits/state.hh>
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/discretization/box/elementsolution.hh>
#include <dumux/discretization/cellcentered/elementsolution.hh>
#include "evalsolution.hh"
namespace Dumux {
/*!
......@@ -50,40 +55,52 @@ auto evalGradients(const Element& element,
const typename Element::Geometry& geometry,
const typename FVElementGeometry::FVGridGeometry& fvGridGeometry,
const BoxElementSolution<FVElementGeometry, PrimaryVariables>& elemSol,
const typename Element::Geometry::GlobalCoordinate& globalPos)
const typename Element::Geometry::GlobalCoordinate& globalPos,
bool ignoreState = false)
{
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// determine if all states are the same at all vertices
using HasState = decltype(isValid(Detail::hasState())(elemSol[0]));
bool allStatesEqual = ignoreState || Detail::allStatesEqual(elemSol, HasState{});
// evaluate gradients using the local finite element basis
const auto& localBasis = fvGridGeometry.feCache().get(geometry.type()).localBasis();
if (allStatesEqual)
{
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
// evaluate the shape function gradients at the scv center
using ShapeJacobian = typename std::decay_t< decltype(localBasis) >::Traits::JacobianType;
const auto localPos = geometry.local(globalPos);
std::vector< ShapeJacobian > shapeJacobian;
localBasis.evaluateJacobian(localPos, shapeJacobian);
// evaluate gradients using the local finite element basis
const auto& localBasis = fvGridGeometry.feCache().get(geometry.type()).localBasis();
// the inverse transposed of the jacobian matrix
const auto jacInvT = geometry.jacobianInverseTransposed(localPos);
// evaluate the shape function gradients at the scv center
using ShapeJacobian = typename std::decay_t< decltype(localBasis) >::Traits::JacobianType;
const auto localPos = geometry.local(globalPos);
std::vector< ShapeJacobian > shapeJacobian;
localBasis.evaluateJacobian(localPos, shapeJacobian);
// interpolate the gradients
Dune::FieldVector<GlobalPosition, PrimaryVariables::dimension> result( GlobalPosition(0.0) );
for (int i = 0; i < element.subEntities(Element::Geometry::mydimension); ++i)
{
// the global shape function gradient
GlobalPosition gradN;
jacInvT.mv(shapeJacobian[i][0], gradN);
// the inverse transposed of the jacobian matrix
const auto jacInvT = geometry.jacobianInverseTransposed(localPos);
// add gradient to global privar gradients
for (unsigned int pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
// interpolate the gradients
Dune::FieldVector<GlobalPosition, PrimaryVariables::dimension> result( GlobalPosition(0.0) );
for (int i = 0; i < element.subEntities(Element::Geometry::mydimension); ++i)
{
GlobalPosition tmp(gradN);
tmp *= elemSol[i][pvIdx];
result[pvIdx] += tmp;
// the global shape function gradient
GlobalPosition gradN;
jacInvT.mv(shapeJacobian[i][0], gradN);
// add gradient to global privar gradients
for (unsigned int pvIdx = 0; pvIdx < PrimaryVariables::dimension; ++pvIdx)
{
GlobalPosition tmp(gradN);
tmp *= elemSol[i][pvIdx];
result[pvIdx] += tmp;
}
}
}
return result;
return result;
}
else
{
DUNE_THROW(Dune::NotImplemented, "Element vertices have different phase states. Enforce calculation by setting ignoreState to true.");
}
}
/*!
......
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