Commit ec5a8c79 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[implicit] store Geometry if it is re-used

Since grid implementations like Yasp calculate the Geometries on the
fly, every call to entity.geometry() may involve computation. This patch
stores a Geometry when it can be used later on in the same block. Only
for the implicit models, decoupled contains a lot more.

For simple models this actually pays off: for 1p in 3d on a larger grid,
the running time could be reduced from 15.5 to 14 s.

Reviewed by Timo.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14284 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 5e11ae86
......@@ -177,8 +177,9 @@ namespace Dumux
Scalar stabilizationTerm(0.0);
if(withStabilization_){
// calculate distance h between nodes i and j
DimVector hVec = this->element_().geometry().corner(fluxVars.face().j)
- this->element_().geometry().corner(fluxVars.face().i);
const auto geometry = this->element_().geometry();
DimVector hVec = geometry.corner(fluxVars.face().j)
- geometry.corner(fluxVars.face().i);
Scalar h = hVec.two_norm();
stabilizationTerm = (h * h) /
(4 * (fluxVars.lambda()
......
......@@ -219,7 +219,8 @@ public:
Traits::LocalBasisType::Traits::RangeType RT_P;
// select quadrature rule for the element geometry type and with the order=qorder
Dune::GeometryType geomType = eg.geometry().type();
const auto geometry = eg.geometry();
Dune::GeometryType geomType = geometry.type();
const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(geomType,qorder);
// loop over quadrature points
......@@ -232,7 +233,7 @@ public:
// get inverse transposed jacobian for quadrature point
const JacobianInverseTransposed jacobian = eg.geometry().jacobianInverseTransposed(it->position());
const JacobianInverseTransposed jacobian = geometry.jacobianInverseTransposed(it->position());
// calculate shape function gradients at the quadrature point in global coordinates. This is done
// by multiplying the reference element shape functions with the inverse transposed jacobian
......@@ -297,7 +298,7 @@ public:
RT_P pn = pw + MaterialLaw::pc(materialParams, sw);
RT_P pEff;
const GlobalPosition& globalPos = eg.geometry().global(it->position());
const GlobalPosition& globalPos = geometry.global(it->position());
// calculate change in effective pressure with respect to initial conditions pInit (pInit is negativ)
pEff = pw*sw + pn*sn + model_.problem().pInit(globalPos, it->position(), eg.entity());
......@@ -339,7 +340,7 @@ public:
RF rhoDiff = volVars.density(nPhaseIdx) - volVars.density(wPhaseIdx);
// geometric weight need for quadrature rule evaluation (numerical integration)
RF qWeight = it->weight() * eg.geometry().integrationElement(it->position());
RF qWeight = it->weight() * geometry.integrationElement(it->position());
// evaluate basis functions
std::vector<RT_V> vBasis(dispSize);
......@@ -399,7 +400,7 @@ public:
// position of quadrature point in local coordinates of element
DimVector local = isIt->geometryInInside().global(it->position());
GlobalPosition globalPos = eg.geometry().global(local);
GlobalPosition globalPos = geometry.global(local);
// evaluate boundary condition type
BoundaryTypes boundaryTypes;
......@@ -465,7 +466,7 @@ public:
// this doesn't work: DimVector local = isIt->geometryInInside().global(face_refElement.position(j,codim-1));
DimVector local = refElement.template geometry<1>(fIdx).global(face_refElement.position(j, codim-1));
GlobalPosition globalPos = eg.geometry().global(local);
GlobalPosition globalPos = geometry.global(local);
// evaluate boundary condition type
BoundaryTypes boundaryTypes;
......
......@@ -429,8 +429,10 @@ public:
effPorosity[eIdx] +=elemVolVars[scvIdx].effPorosity / numScv;
};
const GlobalPosition& cellCenter = eIt->geometry().center();
const GlobalPosition& cellCenterLocal = eIt->geometry().local(cellCenter);
const auto geometry = eIt->geometry();
const GlobalPosition& cellCenter = geometry.center();
const GlobalPosition& cellCenterLocal = geometry.local(cellCenter);
deltaEffPressure[eIdx] = effectivePressure[eIdx] + this->problem().pInit(cellCenter, cellCenterLocal, *eIt);
// determin changes in effective stress from current solution
......@@ -439,7 +441,7 @@ public:
displacementLFS.child(0).finiteElement().localBasis().evaluateJacobian(cellCenterLocal, vRefShapeGradient);
// get jacobian to transform the gradient to physical element
const JacobianInverseTransposed jacInvT = eIt->geometry().jacobianInverseTransposed(cellCenterLocal);
const JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(cellCenterLocal);
std::vector < Dune::FieldVector<RF, dim> > vShapeGradient(dispSize);
for (size_t i = 0; i < dispSize; i++) {
vShapeGradient[i] = 0.0;
......
......@@ -331,8 +331,9 @@ protected:
}
else {
// use two-point gradients
const GlobalPosition &globalPosI = element.geometry().corner(face().i);
const GlobalPosition &globalPosJ = element.geometry().corner(face().j);
const auto geometry = element.geometry();
const GlobalPosition &globalPosI = geometry.corner(face().i);
const GlobalPosition &globalPosJ = geometry.corner(face().j);
tmp = globalPosI;
tmp -= globalPosJ;
Scalar dist = tmp.two_norm();
......
......@@ -201,11 +201,13 @@ private:
int i = this->face().i;
int j = this->face().j;
const auto geometry = element.geometry();
// compute sum of pressure gradients for each phase
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
{
const GlobalPosition localIdx_i = element.geometry().corner(i);
const GlobalPosition localIdx_j = element.geometry().corner(j);
const GlobalPosition localIdx_i = geometry.corner(i);
const GlobalPosition localIdx_j = geometry.corner(j);
isFracture_ = problem.spatialParams().isEdgeFracture(element, fIdx);
......
......@@ -130,7 +130,8 @@ public:
/*
* Calculate the fracture volume fraction wf = 0.5 * Fwidth * 0.5 * Length
*/
Dune::GeometryType geomType = elem.geometry().type();
const auto geometry = elem.geometry();
Dune::GeometryType geomType = geometry.type();
const ReferenceElement &refElement = ReferenceElements::general(geomType);
Scalar vol; //subcontrol volume
......@@ -147,8 +148,8 @@ public:
{
Scalar fracture_width = this->problem_().spatialParams().fractureWidth(elem, fIdx);
const GlobalPosition global_i = elem.geometry().corner(i);
const GlobalPosition global_j = elem.geometry().corner(j);
const GlobalPosition global_i = geometry.corner(i);
const GlobalPosition global_j = geometry.corner(j);
GlobalPosition diff_ij = global_j;
diff_ij -=global_i;
Scalar fracture_length = 0.5*diff_ij.two_norm();
......@@ -165,7 +166,7 @@ public:
storageFracture[nPhaseIdx] = 0.0;
storageMatrix[wPhaseIdx] = 0.0;
storageMatrix[nPhaseIdx] = 0.0;
// const GlobalPosition &globalPos = elem.geometry().corner(scvIdx);
// const GlobalPosition &globalPos = geometry.corner(scvIdx);
Scalar dsm_dsf = volVars.dsm_dsf();
if (!this->problem_().useInterfaceCondition())
......
......@@ -258,7 +258,8 @@ public:
int scvIdx)
{
Scalar volSCVFracture;
Dune::GeometryType geomType = element.geometry().type();
const auto geometry = element.geometry();
Dune::GeometryType geomType = geometry.type();
const ReferenceElement &refElement = ReferenceElements::general(geomType);
for (int fIdx=0; fIdx<refElement.size(1); fIdx++)
......@@ -272,8 +273,8 @@ public:
{
Scalar fracture_width = problem.spatialParams().fractureWidth();
const GlobalPosition global_i = element.geometry().corner(i);
const GlobalPosition global_j = element.geometry().corner(j);
const GlobalPosition global_i = geometry.corner(i);
const GlobalPosition global_j = geometry.corner(j);
GlobalPosition diff_ij = global_j;
diff_ij -=global_i;
//fracture length in the subcontrol volume is half d_ij
......
......@@ -100,7 +100,7 @@ public:
void updateInner(const Element& element)
{
const Geometry& geometry = element.geometry();
const Geometry geometry = element.geometry();
elementVolume = geometry.volume();
elementGlobal = geometry.center();
......@@ -126,7 +126,7 @@ public:
{
updateInner(element);
const Geometry& geometry = element.geometry();
const Geometry geometry = element.geometry();
bool onBoundary = false;
......@@ -134,6 +134,8 @@ public:
IntersectionIterator isEndIt = gridView.iend(element);
for (IntersectionIterator isIt = gridView.ibegin(element); isIt != isEndIt; ++isIt)
{
const auto isGeometry = isIt->geometry();
// neighbor information and inner cvf data:
if (isIt->neighbor())
{
......@@ -147,11 +149,12 @@ public:
scvFace.i = 0;
scvFace.j = scvfIdx + 1;
scvFace.ipGlobal = isIt->geometry().center();
scvFace.ipGlobal = isGeometry.center();
scvFace.ipLocal = geometry.local(scvFace.ipGlobal);
scvFace.normal = isIt->centerUnitOuterNormal();
scvFace.normal *= isIt->geometry().volume();
scvFace.area = isIt->geometry().volume();
Scalar volume = isGeometry.volume();
scvFace.normal *= volume;
scvFace.area = volume;
GlobalPosition distVec = elementGlobal
- neighbors[scvfIdx+1]->geometry().center();
......@@ -179,11 +182,12 @@ public:
int bfIdx = isIt->indexInInside();
SubControlVolumeFace& bFace = boundaryFace[bfIdx];
bFace.ipGlobal = isIt->geometry().center();
bFace.ipGlobal = isGeometry.center();
bFace.ipLocal = geometry.local(bFace.ipGlobal);
bFace.normal = isIt->centerUnitOuterNormal();
bFace.normal *= isIt->geometry().volume();
bFace.area = isIt->geometry().volume();
Scalar volume = isGeometry.volume();
bFace.normal *= volume;
bFace.area = volume;
bFace.i = 0;
bFace.j = 0;
......
......@@ -148,7 +148,9 @@ public:
{
if (velocityOutput_)
{
Dune::GeometryType geomType = element.geometry().type();
const auto geometry = element.geometry();
Dune::GeometryType geomType = geometry.type();
const ReferenceElement &referenceElement
= ReferenceElements::general(geomType);
......@@ -157,7 +159,7 @@ public:
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed jacobianT2 =
element.geometry().jacobianTransposed(localPos);
geometry.jacobianTransposed(localPos);
if (isBox)
{
......@@ -172,7 +174,7 @@ public:
// Transformation of the global normal vector to normal vector in the reference element
const typename Element::Geometry::JacobianTransposed jacobianT1 =
element.geometry().jacobianTransposed(localPosIP);
geometry.jacobianTransposed(localPosIP);
FluxVariables fluxVars(problem_,
element,
......@@ -213,7 +215,7 @@ public:
Dune::FieldVector<CoordScalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(scvVelocities[scvIdx], scvVelocity);
scvVelocity /= element.geometry().integrationElement(localPos)*cellNum_[vIdxGlobal];
scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
// add up the wetting phase subcontrolvolume velocities for each vertex
velocity[vIdxGlobal] += scvVelocity;
}
......@@ -319,7 +321,7 @@ public:
Dune::FieldVector<Scalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
scvVelocity /= element.geometry().integrationElement(localPos);
scvVelocity /= geometry.integrationElement(localPos);
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
int eIdxGlobal = problem_.elementMapper().index(element);
......
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