From da25ffc8ab94108f7d7c964f0e6749da9be002ce Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 16 Jan 2015 09:39:43 +0000
Subject: [PATCH] [reverse] reverse unintended part of r14070
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14071 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/box/boxfvelementgeometry.hh | 4 +-
.../implicit/common/implicitvelocityoutput.hh | 203 +++++-------------
2 files changed, 50 insertions(+), 157 deletions(-)
diff --git a/dumux/implicit/box/boxfvelementgeometry.hh b/dumux/implicit/box/boxfvelementgeometry.hh
index 1bd6fe5729..492aa42f87 100644
--- a/dumux/implicit/box/boxfvelementgeometry.hh
+++ b/dumux/implicit/box/boxfvelementgeometry.hh
@@ -429,8 +429,8 @@ class BoxFVElementGeometry
{0, 2, 0, 1, 3, 2}
};
static const int edgeToFacePyramid[2][8] = {
- {0, 2, 3, 0, 1, 3, 4, 2},
- {1, 0, 0, 4, 3, 2, 1, 4}
+ {1, 2, 3, 4, 1, 3, 4, 2},
+ {0, 0, 0, 0, 3, 2, 1, 4}
};
static const int edgeToFacePrism[2][9] = {
{1, 0, 2, 0, 3, 2, 4, 1, 4},
diff --git a/dumux/implicit/common/implicitvelocityoutput.hh b/dumux/implicit/common/implicitvelocityoutput.hh
index f3d7a88ee3..5de58cd233 100644
--- a/dumux/implicit/common/implicitvelocityoutput.hh
+++ b/dumux/implicit/common/implicitvelocityoutput.hh
@@ -19,7 +19,7 @@
/*!
* \file
*
- * \brief velocity output for implicit (porous media) models
+ * \brief Adaption of the BOX scheme to the two-phase two-component flow model.
*/
#ifndef DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
#define DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
@@ -47,7 +47,6 @@ class ImplicitVelocityOutput
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::ctype CoordScalar;
@@ -64,7 +63,6 @@ class ImplicitVelocityOutput
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
@@ -80,31 +78,41 @@ public:
{
// check, if velocity output can be used (works only for cubes so far)
velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
- if (velocityOutput_)
+ if (velocityOutput_ && isBox)
{
- if (isBox)
+ cellNum_.assign(problem_.gridView().size(dofCodim), 0);
+ }
+
+ ElementIterator eIt = problem_.gridView().template begin<0>();
+ ElementIterator eEndIt = problem_.gridView().template end<0>();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ if (eIt->geometry().type().isCube() == false)
+ {
+ velocityOutput_ = false;
+ }
+
+ if (velocityOutput_ && isBox)
{
- cellNum_.assign(problem_.gridView().size(dofCodim), 0);
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem_.gridView(), *eIt);
- ElementIterator eIt = problem_.gridView().template begin<0>();
- ElementIterator eEndIt = problem_.gridView().template end<0>();
- for (; eIt != eEndIt; ++eIt)
+ // transform vertex velocities from local to global coordinates
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
- FVElementGeometry fvGeometry;
- fvGeometry.update(problem_.gridView(), *eIt);
-
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
- int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
+ int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
#else
- int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
+ int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
#endif
- cellNum_[vIdxGlobal] += 1;
- }
+
+ cellNum_[vIdxGlobal] += 1;
}
}
}
+
+ if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
+ std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
}
bool enableOutput()
@@ -116,17 +124,16 @@ public:
void calculateVelocity(VelocityVector& velocity,
const ElementVolumeVariables& elemVolVars,
const FVElementGeometry& fvGeometry,
- const Element& element,
+ const Element& element,
int phaseIdx)
{
if (velocityOutput_)
{
- Dune::GeometryType geomType = element.geometry().type();
- const ReferenceElement &referenceElement
- = ReferenceElements::general(geomType);
+ // calculate vertex velocities
+ GlobalPosition tmpVelocity(0.0);
- const Dune::FieldVector<Scalar, dim>& localPos
- = referenceElement.position(0, 0);
+ const Dune::FieldVector<Scalar, dim>& localPos =
+ ReferenceElements::general(element.geometry().type()).position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed jacobianT2 =
@@ -134,8 +141,8 @@ public:
if (isBox)
{
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > ScvVelocities;
- ScvVelocities scvVelocities(fvGeometry.numScv);
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
+ SCVVelocities scvVelocities(pow(2,dim));
scvVelocities = 0;
for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++)
@@ -157,17 +164,17 @@ public:
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
- localNormal /= localNormal.two_norm();
+ // note only works for cubes
+ const Scalar localArea = pow(2,-(dim-1));
- // area of the subcontrolvolume face in the reference element
- Scalar localArea = scvfReferenceArea_(geomType, fIdx);
+ localNormal /= localNormal.two_norm();
- // get the volume flux divided by the area of the
- // subcontrolvolume face in the reference element
+ // Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
+ // face in the reference element.
Scalar flux = fluxVars.volumeFlux(phaseIdx) / localArea;
- // transform the volume flux into a velocity vector
- GlobalPosition tmpVelocity = localNormal;
+ // transform the normal Darcy velocity into a vector
+ tmpVelocity = localNormal;
tmpVelocity *= flux;
scvVelocities[fluxVars.face().i] += tmpVelocity;
@@ -193,11 +200,7 @@ public:
}
else
{
-#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
- std::vector<Scalar> scvfFluxes(element.subEntities(1), 0);
-#else
- std::vector<Scalar> scvfFluxes(element.template count<1>(), 0);
-#endif
+ Dune::FieldVector<Scalar, 2*dim> scvVelocities(0.0);
int fIdxInner = 0;
IntersectionIterator isEndIt = problem_.gridView().iend(element);
@@ -214,7 +217,8 @@ public:
fIdxInner,
elemVolVars);
- scvfFluxes[fIdx] += fluxVars.volumeFlux(phaseIdx);
+ Scalar flux = fluxVars.volumeFlux(phaseIdx);
+ scvVelocities[fIdx] += flux;
fIdxInner++;
}
@@ -226,68 +230,14 @@ public:
fIdx,
elemVolVars,true);
- scvfFluxes[fIdx] = fluxVars.volumeFlux(phaseIdx);
- }
- }
-
- // Correct boundary fluxes in case of Neumann conditions.
- // They are simply set to an average of the inner fluxes. In
- // this general setting, it would be very difficult to
- // calculate correct phase, i.e., volume, fluxes from arbitrary
- // Neumann conditions.
- if (element.hasBoundaryIntersections())
- {
- IntersectionIterator isEndIt = problem_.gridView().iend(element);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(element);
- isIt != isEndIt; ++isIt)
- {
- if (isIt->boundary())
- {
- BoundaryTypes bcTypes;
- problem_.boundaryTypes(bcTypes, *isIt);
-
- if (bcTypes.hasNeumann())
- {
- int fIdx = isIt->indexInInside();
-
- // cubes
- if (dim == 1 || geomType.isCube()){
- int fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
- scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
- }
- // simplices
- else if (geomType.isSimplex()) {
- scvfFluxes[fIdx] = 0;
- }
- }
- }
+ Scalar flux = fluxVars.volumeFlux(phaseIdx);
+ scvVelocities[fIdx] = flux;
}
}
- Dune::FieldVector < Scalar, dim > refVelocity;
- // cubes
- if (dim == 1 || geomType.isCube()){
- for (int i = 0; i < dim; i++)
- {
- refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
- }
- }
- // simplices
- else if (geomType.isSimplex()) {
- for (int dimIdx = 0; dimIdx < dim; dimIdx++)
- {
- refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
- for (int fIdx = 0; fIdx < dim + 1; fIdx++)
- {
- refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
- }
- }
- }
- // 3D prism and pyramids
- else {
- DUNE_THROW(Dune::NotImplemented,
- "velocity output for cell-centered and prism/pyramid");
- }
+ Dune::FieldVector < Scalar, dim > refVelocity(0);
+ for (int i = 0; i < dim; i++)
+ refVelocity[i] = 0.5 * (scvVelocities[2*i + 1] - scvVelocities[2*i]);
Dune::FieldVector<Scalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
@@ -305,63 +255,6 @@ public:
} // velocity output
}
-private:
- // The area of a subcontrolvolume face in a reference element.
- // The 3d non-cube values have been calculated with quadrilateralArea3D
- // of boxfvelementgeometry.hh.
- static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx)
- {
- if (dim == 1 || geomType.isCube())
- {
- return 1.0/(1 << (dim-1));
- }
- else if (geomType.isTriangle())
- {
- static const Scalar faceToArea[] = {0.372677996249965,
- 0.372677996249965,
- 0.235702260395516};
- return faceToArea[fIdx];
- }
- else if (geomType.isTetrahedron())
- {
- static const Scalar faceToArea[] = {0.102062072615966,
- 0.102062072615966,
- 0.058925565098879,
- 0.102062072615966,
- 0.058925565098879,
- 0.058925565098879};
- return faceToArea[fIdx];
- }
- else if (geomType.isPyramid())
- {
- static const Scalar faceToArea[] = {0.130437298687488,
- 0.130437298687488,
- 0.130437298687488,
- 0.130437298687488,
- 0.150923085635624,
- 0.1092906420717,
- 0.1092906420717,
- 0.0781735959970571};
- return faceToArea[fIdx];
- }
- else if (geomType.isPrism())
- {
- static const Scalar faceToArea[] = {0.166666666666667,
- 0.166666666666667,
- 0.166666666666667,
- 0.186338998124982,
- 0.186338998124982,
- 0.117851130197758,
- 0.186338998124982,
- 0.186338998124982,
- 0.117851130197758};
- return faceToArea[fIdx];
- }
- else {
- DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType");
- }
- }
-
protected:
const Problem& problem_;
bool velocityOutput_;
--
GitLab