Commit 45a86958 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

added caching to all local assemblers

parent b01c33f7
Pipeline #1270 passed with stage
in 20 minutes and 58 seconds
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/LocalOperator.hpp> #include <amdis/LocalOperator.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/Utility.hpp> #include <amdis/common/Utility.hpp>
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -53,30 +54,34 @@ namespace AMDiS ...@@ -53,30 +54,34 @@ namespace AMDiS
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType; using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType; using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
using RangeType = typename RowLocalBasisType::Traits::RangeType;
using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode); auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
for (auto const& qp : quad) {
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight(); const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
thread_local std::vector<RangeType> shapeValues; // the values of the shape functions on the reference element at the quadrature point
rowLocalFE.localBasis().evaluateFunction(local, shapeValues); auto const& shapeValues = shapeValuesCache[iq];
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients; using WorldVector = FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> colGradients;
colGradients.resize(shapeGradients.size()); colGradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i) for (std::size_t i = 0; i < colGradients.size(); ++i)
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -50,32 +51,37 @@ namespace AMDiS ...@@ -50,32 +51,37 @@ namespace AMDiS
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType; using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType; using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
using RangeType = typename RowLocalBasisType::Traits::RangeType;
using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode); auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
for (auto const& qp : quad) {
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = context.integrationElement(qp.position()) * qp.weight(); const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
const auto b = Super::coefficient(local); const auto b = Super::coefficient(local);
thread_local std::vector<RangeType> shapeValues; // the values of the shape functions on the reference element at the quadrature point
rowLocalFE.localBasis().evaluateFunction(local, shapeValues); auto const& shapeValues = shapeValuesCache[iq];
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients; using WorldVector = FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> colGradients;
colGradients.resize(shapeGradients.size()); colGradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i) for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], colGradients[i]); jacobian.mv(shapeGradients[i][0], colGradients[i]);
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -54,29 +55,32 @@ namespace AMDiS ...@@ -54,29 +55,32 @@ namespace AMDiS
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType; using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType; using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
using RangeType = typename RowLocalBasisType::Traits::RangeType;
using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode); auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
for (auto const& qp : quad) {
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
assert(jacobian.M() == Context::dim); assert(jacobian.M() == Context::dim);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = context.integrationElement(qp.position()) * qp.weight(); const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
const auto c = Super::coefficient(local); const auto c = Super::coefficient(local);
thread_local std::vector<RangeType> shapeValues; // the values of the shape functions on the reference element at the quadrature point
rowLocalFE.localBasis().evaluateFunction(local, shapeValues); auto const& shapeValues = shapeValuesCache[iq];
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<RangeFieldType> colPartial; thread_local std::vector<RangeFieldType> colPartial;
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -53,31 +54,36 @@ namespace AMDiS ...@@ -53,31 +54,36 @@ namespace AMDiS
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType; using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType; using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
using RangeType = typename RowLocalBasisType::Traits::RangeType;
using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
using JacobianType = typename ColLocalBasisType::Traits::JacobianType;
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode); auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
for (auto const& qp : quad) {
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight(); const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
thread_local std::vector<RangeType> shapeValues; // the values of the shape functions on the reference element at the quadrature point
rowLocalFE.localBasis().evaluateFunction(local, shapeValues); auto const& shapeValues = shapeValuesCache[iq];
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
colLocalFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > colGradients; using WorldVector = FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> colGradients;
colGradients.resize(shapeGradients.size()); colGradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i) for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], colGradients[i]); jacobian.mv(shapeGradients[i][0], colGradients[i]);
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -73,34 +74,38 @@ namespace AMDiS ...@@ -73,34 +74,38 @@ namespace AMDiS
using RowFieldType = typename RowLocalBasisType::Traits::RangeFieldType; using RowFieldType = typename RowLocalBasisType::Traits::RangeFieldType;
using ColFieldType = typename ColLocalBasisType::Traits::RangeFieldType; using ColFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
using RowJacobianType = typename RowLocalBasisType::Traits::JacobianType;
using ColJacobianType = typename ColLocalBasisType::Traits::JacobianType;
for (auto const& qp : quad) { LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& rowGradientsCache = rowCache.evaluateJacobianAtQp(context, quad);
auto const& colGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight(); const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<RowJacobianType> rowShapeGradients; auto const& rowShapeGradients = rowGradientsCache[iq];
rowLocalFE.localBasis().evaluateJacobian(local, rowShapeGradients); auto const& colShapeGradients = colGradientsCache[iq];
thread_local std::vector<ColJacobianType> colShapeGradients;
colLocalFE.localBasis().evaluateJacobian(local, colShapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RowFieldType,Context::dow> > rowGradients; using RowWorldVector = FieldVector<RowFieldType,Context::dow>;
thread_local std::vector<RowWorldVector> rowGradients;
rowGradients.resize(rowShapeGradients.size()); rowGradients.resize(rowShapeGradients.size());
for (std::size_t i = 0; i < rowGradients.size(); ++i) for (std::size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowShapeGradients[i][0], rowGradients[i]); jacobian.mv(rowShapeGradients[i][0], rowGradients[i]);
thread_local std::vector<Dune::FieldVector<ColFieldType,Context::dow> > colGradients; using ColWorldVector = FieldVector<ColFieldType,Context::dow>;
thread_local std::vector<ColWorldVector> colGradients;
colGradients.resize(colShapeGradients.size()); colGradients.resize(colShapeGradients.size());
for (std::size_t i = 0; i < colGradients.size(); ++i) for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colShapeGradients[i][0], colGradients[i]); jacobian.mv(colShapeGradients[i][0], colGradients[i]);
...@@ -130,25 +135,28 @@ namespace AMDiS ...@@ -130,25 +135,28 @@ namespace AMDiS
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
using JacobianType = typename LocalBasisType::Traits::JacobianType;
for (auto const& qp : quad) { LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = Super::coefficient(local) * context.integrationElement(qp.position()) * qp.weight(); const auto factor = Super::coefficient(local) * context.integrationElement(quad[iq].position()) * quad[iq].weight();
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients; using WorldVector = FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients;
gradients.resize(shapeGradients.size()); gradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i) for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]); jacobian.mv(shapeGradients[i][0], gradients[i]);
......
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/GridFunctionOperator.hpp> #include <amdis/GridFunctionOperator.hpp>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/ValueCategory.hpp> #include <amdis/common/ValueCategory.hpp>
...@@ -72,26 +73,29 @@ namespace AMDiS ...@@ -72,26 +73,29 @@ namespace AMDiS
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
using JacobianType = typename LocalBasisType::Traits::JacobianType;
for (auto const& qp : quad) { LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = context.integrationElement(qp.position()) * qp.weight(); const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
const auto exprValue = Super::coefficient(local); const auto exprValue = Super::coefficient(local);
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients; using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients;
gradients.resize(shapeGradients.size()); gradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i) for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]); jacobian.mv(shapeGradients[i][0], gradients[i]);
...@@ -135,7 +139,8 @@ namespace AMDiS ...@@ -135,7 +139,8 @@ namespace AMDiS
auto const& shapeGradients = shapeGradientsCache[iq]; auto const& shapeGradients = shapeGradientsCache[iq];
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients; using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients;
gradients.resize(shapeGradients.size()); gradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i) for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]); jacobian.mv(shapeGradients[i][0], gradients[i]);
...@@ -168,25 +173,26 @@ namespace AMDiS ...@@ -168,25 +173,26 @@ namespace AMDiS
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
using JacobianType = typename LocalBasisType::Traits::JacobianType;
for (auto const& qp : quad) { LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(qp.position()); decltype(auto) local = context.local(quad[iq].position());
// The transposed inverse Jacobian of the map from the reference element to the element // The transposed inverse Jacobian of the map from the reference element to the element
const auto jacobian = context.geometry().jacobianInverseTransposed(local); const auto jacobian = context.geometry().jacobianInverseTransposed(local);
// The multiplicative factor in the integral transformation formula // The multiplicative factor in the integral transformation formula
const auto factor = context.integrationElement(qp.position()) * qp.weight(); const auto factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
const auto exprValue = Super::coefficient(local); const auto exprValue = Super::coefficient(local);
// The gradients of the shape functions on the reference element // The gradients of the shape functions on the reference element
thread_local std::vector<JacobianType> shapeGradients; auto const& shapeGradients = shapeGradientsCache[iq];
localFE.localBasis().evaluateJacobian(local, shapeGradients);
// Compute the shape function gradients on the real element // Compute the shape function gradients on the real element
thread_local std::vector<Dune::FieldVector<RangeFieldType,Context::dow> > gradients; using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients;
gradients.resize(shapeGradients.size()); gradients.resize(shapeGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i) for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]); jacobian.mv(shapeGradients[i][0], gradients[i]);
......