Commit bfde81bd authored by Nicolas Schwenck's avatar Nicolas Schwenck
Browse files

fixed wrong velocity output calculation for non-cube elements in decoupled models

(except MPFA which anyway only works on quadrilateral grids)


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13584 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent e27f3776
......@@ -123,7 +123,15 @@ public:
CellData& cellData = problem_.variables().cellData(globalIdx);
Dune::FieldVector<Scalar, 2*dim> flux(0);
const typename Element::Geometry& geometry = eIt->geometry();
// get corresponding reference element
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::ReferenceElement< Scalar , dim > & refElement =
ReferenceElements::general( geometry.type() );
const int numberOfFaces=refElement.size(1);
std::vector<Scalar> flux(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator
isEndIt = problem_.gridView().iend(*eIt);
......@@ -136,19 +144,31 @@ public:
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
Dune::FieldVector<Scalar, dim> refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
const typename Element::Geometry& geometry = eIt->geometry();
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::FieldVector<Scalar, dim>& localPos
= ReferenceElements::general(geometry.type()).position(0, 0);
const Dune::FieldVector<Scalar, dim>& localPos
= refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed& jacobianT =
const typename Element::Geometry::JacobianTransposed& jacobianT =
geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
......
......@@ -198,8 +198,16 @@ public:
CellData& cellData = problem_.variables().cellData(globalIdx);
Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
const typename Element::Geometry& geometry = eIt->geometry();
// get corresponding reference element
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::ReferenceElement< Scalar , dim > & refElement =
ReferenceElements::general( geometry.type() );
const int numberOfFaces=refElement.size(1);
std::vector<Scalar> fluxW(numberOfFaces,0);
std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
......@@ -213,31 +221,61 @@ public:
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const Dune::FieldVector<Scalar, dim> localPos =
ReferenceElements::general(eIt->geometry().type()).position(0, 0);
refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const JacobianTransposed jacobianT =
eIt->geometry().jacobianTransposed(localPos);
geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
Dune::FieldVector < Scalar, dim > elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
velocity[globalIdx] = elementVelocity;
refVelocity = 0;
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
velocitySecondPhase[globalIdx] = elementVelocity;
}
......
......@@ -280,8 +280,15 @@ public:
(*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
}
Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
const typename Element::Geometry& geometry = eIt->geometry();
// get corresponding reference element
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::ReferenceElement< Scalar , dim > & refElement =
ReferenceElements::general( geometry.type() );
const int numberOfFaces=refElement.size(1);
std::vector<Scalar> fluxW(numberOfFaces,0);
std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
......@@ -295,29 +302,59 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity;
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
const DimVector& localPos = ReferenceElements::general(eIt->geometry().type()).position(0, 0);
DimVector refVelocity(0.0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const JacobianTransposed jacobianT = eIt->geometry().jacobianTransposed(localPos);
const JacobianTransposed jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
(*velocityWetting)[globalIdx] = elementVelocity;
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
(*velocityNonwetting)[globalIdx] = elementVelocity;
}
......
......@@ -292,8 +292,15 @@ public:
(*pressureSecond)[globalIdx] = cellData.pressure(wPhaseIdx);
}
Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
const typename Element::Geometry& geometry = eIt->geometry();
// get corresponding reference element
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::ReferenceElement< Scalar , dim > & refElement =
ReferenceElements::general( geometry.type() );
const int numberOfFaces=refElement.size(1);
std::vector<Scalar> fluxW(numberOfFaces,0);
std::vector<Scalar> fluxNw(numberOfFaces,0);
// run through all intersections with neighbors and boundary
IntersectionIterator isEndIt = problem_.gridView().iend(*eIt);
......@@ -307,29 +314,59 @@ public:
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
DimVector refVelocity;
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
const DimVector& localPos = ReferenceElements::general(eIt->geometry().type()).position(0, 0);
Dune::FieldVector < Scalar, dim > refVelocity(0);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxW[1] - fluxW[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxW[0] - fluxW[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxW[2*k+1] - fluxW[2*k]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const JacobianTransposed jacobianT = eIt->geometry().jacobianTransposed(localPos);
const JacobianTransposed jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
(*velocityWetting)[globalIdx] = elementVelocity;
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*fluxNw[1] - fluxNw[2]);
refVelocity[1] = 1.0/3.0*(2.0*fluxNw[0] - fluxNw[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.type().isCube() ){
for (int k = 0; k < dim; k++)
refVelocity[k] = 0.5 * (fluxNw[2*k+1] - fluxNw[2*k]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
(*velocityNonwetting)[globalIdx] = elementVelocity;
}
......
......@@ -211,8 +211,16 @@ public:
break;
}
Dune::FieldVector < Scalar, 2 * dim > flux(0);
// Berechne Verfeinerungsindikator an allen Zellen
const typename Element::Geometry& geometry = eIt->geometry();
// get corresponding reference element
typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
const Dune::ReferenceElement< Scalar , dim > & refElement =
ReferenceElements::general( geometry.type() );
const int numberOfFaces=refElement.size(1);
std::vector<Scalar> flux(numberOfFaces,0);
// Calculate refinement indicator on all cells
IntersectionIterator isItend = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItend; ++isIt)
{
......@@ -331,19 +339,34 @@ public:
if (useFluxInd_ || usePercentileFlux_)
{
DimVector refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
//2D simplices
if ( dim==2 && geometry.corners() == 3 ) {
refVelocity[0] = 1.0/3.0*(2.0*flux[1] - flux[2]);
refVelocity[1] = 1.0/3.0*(2.0*flux[0] - flux[2]);
}
//3D tetrahedra
else if ( dim==3 && geometry.corners() == 4 ){
DUNE_THROW(Dune::NotImplemented, "velocity output for tetrahedra not implemented");
}
//2D and 3D cubes
else if ( refElement.typ().isCube() ){
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (flux[2*i + 1] - flux[2*i]);
}
//3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
}
const DimVector& localPos =
ReferenceElements::general(eIt->geometry().type()).position(0, 0);
const DimVector& localPos = refElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const DimMatrix jacobianT = eIt->geometry().jacobianTransposed(localPos);
const DimMatrix jacobianT = geometry.jacobianTransposed(localPos);
// calculate the element velocity by the Piola transformation
DimVector elementVelocity(0);
jacobianT.umtv(refVelocity, elementVelocity);
elementVelocity /= eIt->geometry().integrationElement(localPos);
elementVelocity /= geometry.integrationElement(localPos);
Scalar velNorm = elementVelocity.two_norm();
indicatorVectorFlux_[globalIdxI] = std::max(velNorm, indicatorVectorFlux_[globalIdxI]);
......
Markdown is supported
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