Skip to content
Snippets Groups Projects
Commit 7d0d7d4d authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

BoxFVElementGeometry: introduce possibility of two-point flux approximation 

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9004 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 104fe99d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dumux/boxmodels/common/boxfvelementgeometry.hh
+ 90
− 77
View file @ 7d0d7d4d
...@@ -39,6 +39,7 @@ namespace Properties ...@@ -39,6 +39,7 @@ namespace Properties
{ {
NEW_PROP_TAG(GridView); NEW_PROP_TAG(GridView);
NEW_PROP_TAG(Scalar); NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(ImplicitUseTwoPointFlux);
} }
///////////////////// /////////////////////
...@@ -613,8 +614,14 @@ public: ...@@ -613,8 +614,14 @@ public:
numEdges = referenceElement.size(dim-1); numEdges = referenceElement.size(dim-1);
numFaces = (dim<3)?0:referenceElement.size(1); numFaces = (dim<3)?0:referenceElement.size(1);
numSCV = numVertices; numSCV = numVertices;
numFAP = numVertices;
bool useTwoPointFlux
= GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, UseTwoPointFlux);
if (useTwoPointFlux)
numFAP = 2;
else
numFAP = numVertices;
// corners: // corners:
for (int vert = 0; vert < numVertices; vert++) { for (int vert = 0; vert < numVertices; vert++) {
subContVol[vert].local = referenceElement.position(vert, dim); subContVol[vert].local = referenceElement.position(vert, dim);
...@@ -690,26 +697,41 @@ public: ...@@ -690,26 +697,41 @@ public:
elementGlobal, faceCoord[leftFace]); elementGlobal, faceCoord[leftFace]);
} }
// get the global integration point and the Jacobian inverse if (useTwoPointFlux)
subContVolFace[k].ipGlobal = geometry.global(ipLocal); {
Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT = geometry.jacobianInverseTransposed(ipLocal); GlobalPosition distVec = subContVol[i].global;
distVec -= subContVol[j].global;
distVec /= distVec.two_norm2();
// std::cout << "SCV Face " << k << ", i = " << i << ", j = " << j
// << ", ipLocal = " << ipLocal << ", ipGlobal = " << subContVolFace[k].ipGlobal << ", normal = " << subContVolFace[k].normal // gradients using a two-point flux approximation
// << std::endl; for (int idx = 0; idx < 2; idx++)
{
// calculate the shape function gradients subContVolFace[k].grad[idx] = distVec;
//typedef Dune::FieldVector< Dune::FieldVector< CoordScalar, dim >, 1 > ShapeJacobian; subContVolFace[k].shapeValue[idx] = 0.5;
typedef Dune::FieldVector< Scalar, 1 > ShapeValue; }
std::vector<ShapeJacobian> localJac; subContVolFace[k].grad[1] *= -1.0;
std::vector<ShapeValue> shapeVal;
localFiniteElement.localBasis().evaluateJacobian(subContVolFace[k].ipLocal, localJac); subContVolFace[k].fapIndices[0] = subContVolFace[k].i;
localFiniteElement.localBasis().evaluateFunction(subContVolFace[k].ipLocal, shapeVal); subContVolFace[k].fapIndices[1] = subContVolFace[k].j;
for (int vert = 0; vert < numVertices; vert++) { }
jacInvT.mv(localJac[vert][0], subContVolFace[k].grad[vert]); else
subContVolFace[k].shapeValue[vert] = Scalar(shapeVal[vert]); {
subContVolFace[k].fapIndices[vert] = vert; // 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 } // end loop over edges / sub control volume faces
...@@ -773,68 +795,59 @@ public: ...@@ -773,68 +795,59 @@ public:
boundaryFace[bfIdx].shapeValue[vert] = Scalar(shapeVal[vert]); boundaryFace[bfIdx].shapeValue[vert] = Scalar(shapeVal[vert]);
boundaryFace[bfIdx].fapIndices[vert] = 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
// calculate gradients at the center of the scv for (int scvIdx = 0; scvIdx < numVertices; scvIdx++)
for (int scvIdx = 0; scvIdx < numVertices; scvIdx++) if (dim == 2)
if (dim == 2) {
switch (scvIdx)
{ {
// if (!subContVol[scvIdx].inner) case 0:
{ if (numVertices == 4) {
switch (scvIdx) subContVol[scvIdx].localCenter[0] = 0.25;
{ subContVol[scvIdx].localCenter[1] = 0.25;
case 0: }
if (numVertices == 4) { else {
subContVol[scvIdx].localCenter[0] = 0.25; subContVol[scvIdx].localCenter[0] = 1.0/6.0;
subContVol[scvIdx].localCenter[1] = 0.25; subContVol[scvIdx].localCenter[1] = 1.0/6.0;
} }
else { break;
subContVol[scvIdx].localCenter[0] = 1.0/6.0; case 1:
subContVol[scvIdx].localCenter[1] = 1.0/6.0; if (numVertices == 4) {
} subContVol[scvIdx].localCenter[0] = 0.75;
break; subContVol[scvIdx].localCenter[1] = 0.25;
case 1: }
if (numVertices == 4) { else {
subContVol[scvIdx].localCenter[0] = 0.75; subContVol[scvIdx].localCenter[0] = 4.0/6.0;
subContVol[scvIdx].localCenter[1] = 0.25; subContVol[scvIdx].localCenter[1] = 1.0/6.0;
} }
else { break;
subContVol[scvIdx].localCenter[0] = 4.0/6.0; case 2:
subContVol[scvIdx].localCenter[1] = 1.0/6.0; if (numVertices == 4) {
} subContVol[scvIdx].localCenter[0] = 0.25;
break; subContVol[scvIdx].localCenter[1] = 0.75;
case 2: }
if (numVertices == 4) { else {
subContVol[scvIdx].localCenter[0] = 0.25; subContVol[scvIdx].localCenter[0] = 1.0/6.0;
subContVol[scvIdx].localCenter[1] = 0.75; subContVol[scvIdx].localCenter[1] = 4.0/6.0;
}
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;
}
} }
break;
case 3:
subContVol[scvIdx].localCenter[0] = 0.75;
subContVol[scvIdx].localCenter[1] = 0.75;
break;
}
std::vector<ShapeJacobian> localJac; std::vector<ShapeJacobian> localJac;
localFiniteElement.localBasis().evaluateJacobian(subContVol[scvIdx].localCenter, localJac); localFiniteElement.localBasis().evaluateJacobian(subContVol[scvIdx].localCenter, localJac);
Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT = Dune::FieldMatrix<CoordScalar,dim,dim> jacInvT =
geometry.jacobianInverseTransposed(subContVol[scvIdx].localCenter); geometry.jacobianInverseTransposed(subContVol[scvIdx].localCenter);
for (int vert = 0; vert < numVertices; vert++) for (int vert = 0; vert < numVertices; vert++)
jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]); jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
} }
} }
}; };
......
dumux/boxmodels/common/boxproperties.hh
+ 3
− 0
View file @ 7d0d7d4d
...@@ -125,6 +125,9 @@ NEW_PROP_TAG(NumericDifferenceMethod);//DEPRECATED ...@@ -125,6 +125,9 @@ NEW_PROP_TAG(NumericDifferenceMethod);//DEPRECATED
NEW_PROP_TAG(ImplicitEnableHints); NEW_PROP_TAG(ImplicitEnableHints);
NEW_PROP_TAG(EnableHints);//DEPRECATED NEW_PROP_TAG(EnableHints);//DEPRECATED
//! indicates whether two-point flux should be used
NEW_PROP_TAG(ImplicitUseTwoPointFlux);
// mappers from local to global indices // mappers from local to global indices
//! maper for vertices //! maper for vertices
......
dumux/boxmodels/common/boxpropertydefaults.hh
+ 3
− 0
View file @ 7d0d7d4d
...@@ -177,6 +177,9 @@ SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);//DEPRECATED ...@@ -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, ImplicitEnablePartialReassemble, GET_PROP_VALUE(TypeTag, EnablePartialReassemble));
SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);//DEPRECATED 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 the type of a global jacobian matrix from the solution types
SET_PROP(BoxModel, JacobianMatrix) SET_PROP(BoxModel, JacobianMatrix)
{ {
......
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