diff --git a/dumux/boxmodels/common/boxfvelementgeometry.hh b/dumux/boxmodels/common/boxfvelementgeometry.hh
index 35a71dcd5aec79aa1efa0921c02f0add2ad43d27..197f0c2f978c04fd2a81cf8ded1f241ff950b01f 100644
--- a/dumux/boxmodels/common/boxfvelementgeometry.hh
+++ b/dumux/boxmodels/common/boxfvelementgeometry.hh
@@ -39,6 +39,7 @@ namespace Properties
{
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(ImplicitUseTwoPointFlux);
}
/////////////////////
@@ -613,8 +614,14 @@ public:
numEdges = referenceElement.size(dim-1);
numFaces = (dim<3)?0:referenceElement.size(1);
numSCV = numVertices;
- numFAP = numVertices;
-
+
+ bool useTwoPointFlux
+ = GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, UseTwoPointFlux);
+ if (useTwoPointFlux)
+ numFAP = 2;
+ else
+ numFAP = numVertices;
+
// corners:
for (int vert = 0; vert < numVertices; vert++) {
subContVol[vert].local = referenceElement.position(vert, dim);
@@ -690,26 +697,41 @@ public:
elementGlobal, faceCoord[leftFace]);
}
- // get the global integration point and the Jacobian inverse
- subContVolFace[k].ipGlobal = geometry.global(ipLocal);
- Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT = geometry.jacobianInverseTransposed(ipLocal);
-
-
- // std::cout << "SCV Face " << k << ", i = " << i << ", j = " << j
- // << ", ipLocal = " << ipLocal << ", ipGlobal = " << subContVolFace[k].ipGlobal << ", normal = " << subContVolFace[k].normal
- // << std::endl;
-
- // calculate the shape function gradients
- //typedef Dune::FieldVector< Dune::FieldVector< CoordScalar, dim >, 1 > ShapeJacobian;
- typedef Dune::FieldVector< Scalar, 1 > ShapeValue;
- std::vector<ShapeJacobian> localJac;
- std::vector<ShapeValue> shapeVal;
- localFiniteElement.localBasis().evaluateJacobian(subContVolFace[k].ipLocal, localJac);
- localFiniteElement.localBasis().evaluateFunction(subContVolFace[k].ipLocal, shapeVal);
- for (int vert = 0; vert < numVertices; vert++) {
- jacInvT.mv(localJac[vert][0], subContVolFace[k].grad[vert]);
- subContVolFace[k].shapeValue[vert] = Scalar(shapeVal[vert]);
- subContVolFace[k].fapIndices[vert] = vert;
+ if (useTwoPointFlux)
+ {
+ GlobalPosition distVec = subContVol[i].global;
+ distVec -= subContVol[j].global;
+ distVec /= distVec.two_norm2();
+
+ // gradients using a two-point flux approximation
+ for (int idx = 0; idx < 2; idx++)
+ {
+ subContVolFace[k].grad[idx] = distVec;
+ subContVolFace[k].shapeValue[idx] = 0.5;
+ }
+ subContVolFace[k].grad[1] *= -1.0;
+
+ subContVolFace[k].fapIndices[0] = subContVolFace[k].i;
+ subContVolFace[k].fapIndices[1] = subContVolFace[k].j;
+ }
+ else
+ {
+ // get the global integration point and the Jacobian inverse
+ subContVolFace[k].ipGlobal = geometry.global(ipLocal);
+ Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT = geometry.jacobianInverseTransposed(ipLocal);
+
+ // calculate the shape function gradients
+ //typedef Dune::FieldVector< Dune::FieldVector< CoordScalar, dim >, 1 > ShapeJacobian;
+ typedef Dune::FieldVector< Scalar, 1 > ShapeValue;
+ std::vector<ShapeJacobian> localJac;
+ std::vector<ShapeValue> shapeVal;
+ localFiniteElement.localBasis().evaluateJacobian(subContVolFace[k].ipLocal, localJac);
+ localFiniteElement.localBasis().evaluateFunction(subContVolFace[k].ipLocal, shapeVal);
+ for (int vert = 0; vert < numVertices; vert++) {
+ jacInvT.mv(localJac[vert][0], subContVolFace[k].grad[vert]);
+ subContVolFace[k].shapeValue[vert] = Scalar(shapeVal[vert]);
+ subContVolFace[k].fapIndices[vert] = vert;
+ }
}
} // end loop over edges / sub control volume faces
@@ -773,68 +795,59 @@ public:
boundaryFace[bfIdx].shapeValue[vert] = Scalar(shapeVal[vert]);
boundaryFace[bfIdx].fapIndices[vert] = vert;
}
-
- // std::cout << "boundary face " << face << ", vert = " << vertInElement << ", ipLocal = "
- // << boundaryFace[bfIdx].ipLocal << ", ipGlobal = " << boundaryFace[bfIdx].ipGlobal
- // << ", area = " << boundaryFace[bfIdx].area << std::endl;
-
}
}
-
- // calculate gradients at the center of the scv
- for (int scvIdx = 0; scvIdx < numVertices; scvIdx++)
- if (dim == 2)
+ // calculate gradients at the center of the scv
+ for (int scvIdx = 0; scvIdx < numVertices; scvIdx++)
+ if (dim == 2)
+ {
+ switch (scvIdx)
{
-// if (!subContVol[scvIdx].inner)
- {
- switch (scvIdx)
- {
- case 0:
- if (numVertices == 4) {
- subContVol[scvIdx].localCenter[0] = 0.25;
- subContVol[scvIdx].localCenter[1] = 0.25;
- }
- else {
- subContVol[scvIdx].localCenter[0] = 1.0/6.0;
- subContVol[scvIdx].localCenter[1] = 1.0/6.0;
- }
- break;
- case 1:
- if (numVertices == 4) {
- subContVol[scvIdx].localCenter[0] = 0.75;
- subContVol[scvIdx].localCenter[1] = 0.25;
- }
- else {
- subContVol[scvIdx].localCenter[0] = 4.0/6.0;
- subContVol[scvIdx].localCenter[1] = 1.0/6.0;
- }
- break;
- case 2:
- if (numVertices == 4) {
- subContVol[scvIdx].localCenter[0] = 0.25;
- subContVol[scvIdx].localCenter[1] = 0.75;
- }
- else {
- subContVol[scvIdx].localCenter[0] = 1.0/6.0;
- subContVol[scvIdx].localCenter[1] = 4.0/6.0;
- }
- break;
- case 3:
- subContVol[scvIdx].localCenter[0] = 0.75;
- subContVol[scvIdx].localCenter[1] = 0.75;
- break;
- }
+ case 0:
+ if (numVertices == 4) {
+ subContVol[scvIdx].localCenter[0] = 0.25;
+ subContVol[scvIdx].localCenter[1] = 0.25;
+ }
+ else {
+ subContVol[scvIdx].localCenter[0] = 1.0/6.0;
+ subContVol[scvIdx].localCenter[1] = 1.0/6.0;
+ }
+ break;
+ case 1:
+ if (numVertices == 4) {
+ subContVol[scvIdx].localCenter[0] = 0.75;
+ subContVol[scvIdx].localCenter[1] = 0.25;
+ }
+ else {
+ subContVol[scvIdx].localCenter[0] = 4.0/6.0;
+ subContVol[scvIdx].localCenter[1] = 1.0/6.0;
+ }
+ break;
+ case 2:
+ if (numVertices == 4) {
+ subContVol[scvIdx].localCenter[0] = 0.25;
+ subContVol[scvIdx].localCenter[1] = 0.75;
+ }
+ else {
+ subContVol[scvIdx].localCenter[0] = 1.0/6.0;
+ subContVol[scvIdx].localCenter[1] = 4.0/6.0;
}
+ break;
+ case 3:
+ subContVol[scvIdx].localCenter[0] = 0.75;
+ subContVol[scvIdx].localCenter[1] = 0.75;
+ break;
+ }
- std::vector<ShapeJacobian> localJac;
- localFiniteElement.localBasis().evaluateJacobian(subContVol[scvIdx].localCenter, localJac);
+ std::vector<ShapeJacobian> localJac;
+ localFiniteElement.localBasis().evaluateJacobian(subContVol[scvIdx].localCenter, localJac);
- Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT =
- geometry.jacobianInverseTransposed(subContVol[scvIdx].localCenter);
- for (int vert = 0; vert < numVertices; vert++)
- jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
- }
+ Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT =
+ geometry.jacobianInverseTransposed(subContVol[scvIdx].localCenter);
+ for (int vert = 0; vert < numVertices; vert++)
+ jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
+ }
}
};
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index 3dc6801476b4b9ce6e5da24fdcd5eb6368792d65..ed5e67f21f277846c446f585dbf562dc88483e72 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -125,6 +125,9 @@ NEW_PROP_TAG(NumericDifferenceMethod);//DEPRECATED
NEW_PROP_TAG(ImplicitEnableHints);
NEW_PROP_TAG(EnableHints);//DEPRECATED
+//! indicates whether two-point flux should be used
+NEW_PROP_TAG(ImplicitUseTwoPointFlux);
+
// mappers from local to global indices
//! maper for vertices
diff --git a/dumux/boxmodels/common/boxpropertydefaults.hh b/dumux/boxmodels/common/boxpropertydefaults.hh
index 105482e12be55ec433f5931993d6de6c8e5f8db2..6222431015e997b7f5cf00ff8c1bc2df8cede65f 100644
--- a/dumux/boxmodels/common/boxpropertydefaults.hh
+++ b/dumux/boxmodels/common/boxpropertydefaults.hh
@@ -177,6 +177,9 @@ SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);//DEPRECATED
SET_BOOL_PROP(BoxModel, ImplicitEnablePartialReassemble, GET_PROP_VALUE(TypeTag, EnablePartialReassemble));
SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);//DEPRECATED
+// disable two-point-flux by default
+SET_BOOL_PROP(BoxModel, ImplicitUseTwoPointFlux, false);
+
//! Set the type of a global jacobian matrix from the solution types
SET_PROP(BoxModel, JacobianMatrix)
{