Skip to content
GitLab
Commit c089668d authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Some improvements in the assemblers

parent fff818b2
Show whitespace changes
Inline Side-by-side
src/amdis/assembler/FirstOrderTestDivTrialvec.hpp
View file @ c089668d
......@@ -48,8 +48,11 @@ namespace AMDiS
static const std::size_t CHILDREN = ColNode::CHILDREN;
static_assert( CHILDREN == Context::dow, "" );
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -59,8 +62,6 @@ namespace AMDiS
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
......@@ -87,10 +88,10 @@ namespace AMDiS
for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
const auto value = factor * shapeValues[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
for (std::size_t k = 0; k < CHILDREN; ++k) {
const auto local_kj = colNode.child(k).localIndex(j);
elementMatrix[local_i][local_kj] += value * colGradients[j][k];
......
src/amdis/assembler/FirstOrderTestGradTrial.hpp
View file @ c089668d
......@@ -45,8 +45,11 @@ namespace AMDiS
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -56,8 +59,6 @@ namespace AMDiS
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
......@@ -85,10 +86,10 @@ namespace AMDiS
for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto local_j = colNode.localIndex(j);
const auto value = factor * (b * colGradients[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * shapeValues[i];
}
......
src/amdis/assembler/FirstOrderTestPartialTrial.hpp
View file @ c089668d
......@@ -49,8 +49,11 @@ namespace AMDiS
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -60,8 +63,6 @@ namespace AMDiS
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
......@@ -91,10 +92,10 @@ namespace AMDiS
colPartial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto local_j = colNode.localIndex(j);
const auto value = factor * (c * colPartial[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * shapeValues[i];
}
......
src/amdis/assembler/FirstOrderTestvecGradTrial.hpp
View file @ c089668d
......@@ -48,8 +48,11 @@ namespace AMDiS
static const std::size_t CHILDREN = RowNode::CHILDREN;
static_assert( CHILDREN == Context::dow, "" );
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -59,8 +62,6 @@ namespace AMDiS
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQp(context, quad);
auto const& shapeGradientsCache = colCache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
......@@ -87,9 +88,9 @@ namespace AMDiS
for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto value = factor * shapeValues[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto local_j = colNode.localIndex(j);
for (std::size_t k = 0; k < CHILDREN; ++k) {
const auto local_ki = rowNode.child(k).localIndex(i);
......
src/amdis/assembler/SecondOrderDivTestvecDivTrialvec.hpp
View file @ c089668d
......@@ -68,6 +68,8 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.child(0).finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -109,8 +111,8 @@ namespace AMDiS
for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colShapeGradients[i][0], colGradients[i]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t i = 0; i < rowSize; ++i) {
for (std::size_t j = 0; j < colSize; ++j) {
for (std::size_t k = 0; k < CHILDREN; ++k) {
const auto local_ki = rowNode.child(k).localIndex(i);
const auto local_kj = colNode.child(k).localIndex(j);
......@@ -132,6 +134,7 @@ namespace AMDiS
static const std::size_t CHILDREN = Node::CHILDREN;
auto const& localFE = node.child(0).finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
......@@ -161,12 +164,12 @@ namespace AMDiS
jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_ki = node.child(k).localIndex(i);
const auto value = factor * gradients[i][k];
elementMatrix[local_ki][local_ki] += value * gradients[i][k];
for (std::size_t j = i+1; j < localFE.size(); ++j) {
for (std::size_t j = i+1; j < size; ++j) {
const auto local_kj = node.child(k).localIndex(j);
elementMatrix[local_ki][local_kj] += value * gradients[j][k];
elementMatrix[local_kj][local_ki] += value * gradients[j][k];
......
src/amdis/assembler/SecondOrderGradTestGradTrial.hpp
View file @ c089668d
......@@ -70,6 +70,7 @@ namespace AMDiS
ElementMatrix& elementMatrix)
{
auto const& localFE = rowNode.finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
......@@ -99,9 +100,9 @@ namespace AMDiS
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = rowNode.localIndex(i);
for (std::size_t j = 0; j < localFE.size(); ++j) {
for (std::size_t j = 0; j < size; ++j) {
const auto local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
}
......@@ -118,6 +119,7 @@ namespace AMDiS
tag::scalar)
{
auto const& localFE = node.finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
......@@ -145,12 +147,12 @@ namespace AMDiS
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = node.localIndex(i);
elementMatrix[local_i][local_i] += factor * (gradients[i] * gradients[i]);
for (std::size_t j = i+1; j < localFE.size(); ++j) {
for (std::size_t j = i+1; j < size; ++j) {
const auto local_j = node.localIndex(j);
const auto value = factor * (gradients[i] * gradients[j]);
......@@ -170,6 +172,7 @@ namespace AMDiS
tag::matrix)
{
auto const& localFE = node.finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
......@@ -197,9 +200,9 @@ namespace AMDiS
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(shapeGradients[i][0], gradients[i]);
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = node.localIndex(i);
for (std::size_t j = 0; j < localFE.size(); ++j) {
for (std::size_t j = 0; j < size; ++j) {
const auto local_j = node.localIndex(j);
elementMatrix[local_i][local_j] += eval(exprValue, factor, gradients[i], gradients[j]);
}
......
src/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp
View file @ c089668d
......@@ -51,8 +51,11 @@ namespace AMDiS
static_assert(RowNode::isLeaf && ColNode::isLeaf,
"Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -60,8 +63,6 @@ namespace AMDiS
using RowFieldType = typename RowLocalBasisType::Traits::RangeFieldType;
using ColFieldType = typename ColLocalBasisType::Traits::RangeFieldType;
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
LocalBasisCache<RowLocalBasisType> rowCache(rowLocalFE.localBasis());
LocalBasisCache<ColLocalBasisType> colCache(colLocalFE.localBasis());
......@@ -99,10 +100,10 @@ namespace AMDiS
colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
}
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto local_j = colNode.localIndex(j);
const auto value = factor * (c * colPartial[j]);
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * rowPartial[i];
}
......
src/amdis/assembler/StokesOperator.hpp
View file @ c089668d
......@@ -48,6 +48,8 @@ namespace AMDiS
using namespace Dune::Indices;
auto const& quad = this->getQuadratureRule(context.type(), tree, tree);
auto const& velocityLocalFE = tree.child(_0,0).finiteElement();
std::size_t numVelocityLocalFE = velocityLocalFE.size();
......@@ -60,10 +62,6 @@ namespace AMDiS
LocalBasisCache<VelocityLocalBasisType> velocityCache(velocityLocalFE.localBasis());
LocalBasisCache<PressureLocalBasisType> pressureCache(pressureLocalFE.localBasis());
using RangeFieldType = typename VelocityLocalBasisType::Traits::RangeFieldType;
auto const& quad = this->getQuadratureRule(context.type(), tree, tree);
auto const& shapeGradientsCache = velocityCache.evaluateJacobianAtQp(context, quad);
auto const& pressureValuesCache = pressureCache.evaluateFunctionAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
......@@ -80,6 +78,7 @@ namespace AMDiS
auto const& shapeGradients = shapeGradientsCache[iq];
// Compute the shape function gradients on the real element
using RangeFieldType = typename VelocityLocalBasisType::Traits::RangeFieldType;
using WorldVector = Dune::FieldVector<RangeFieldType,Context::dow>;
thread_local std::vector<WorldVector> gradients;
gradients.resize(shapeGradients.size());
......
src/amdis/assembler/ZeroOrderTest.hpp
View file @ c089668d
......@@ -42,12 +42,13 @@ namespace AMDiS
{
static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only");
auto const& quad = this->getQuadratureRule(context.type(), node);
auto const& localFE = node.finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
auto const& quad = this->getQuadratureRule(context.type(), node);
auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context,quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
......@@ -58,7 +59,7 @@ namespace AMDiS
* quad[iq].weight();
auto const& shapeValues = shapeValuesCache[iq];
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = node.localIndex(i);
elementVector[local_i] += factor * shapeValues[i];
}
......
src/amdis/assembler/ZeroOrderTestTrial.hpp
View file @ c089668d
......@@ -64,6 +64,8 @@ namespace AMDiS
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -82,11 +84,11 @@ namespace AMDiS
auto const& rowShapeValues = rowShapeValuesCache[iq];
auto const& colShapeValues = colShapeValuesCache[iq];
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
const auto value = factor * rowShapeValues[i];
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto local_j = colNode.localIndex(j);
elementMatrix[local_i][local_j] += value * colShapeValues[j];
}
......@@ -107,6 +109,7 @@ namespace AMDiS
"Operator can be applied to Leaf-Nodes only.");
auto const& localFE = node.finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
......@@ -121,13 +124,13 @@ namespace AMDiS
auto const& shapeValues = shapeValuesCache[iq];
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto local_i = node.localIndex(i);
const auto value = factor * shapeValues[i];
elementMatrix[local_i][local_i] += value * shapeValues[i];
for (std::size_t j = i+1; j < localFE.size(); ++j) {
for (std::size_t j = i+1; j < size; ++j) {
const auto local_j = node.localIndex(j);
elementMatrix[local_i][local_j] += value * shapeValues[j];
......
src/amdis/assembler/ZeroOrderTestTrialvec.hpp
View file @ c089668d
......@@ -47,9 +47,11 @@ namespace AMDiS
static const std::size_t CHILDREN = ColNode::CHILDREN;
static_assert( std::is_same<typename GridFct::Range, Dune::FieldVector<double, CHILDREN>>::value, "" );
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
auto const& rowLocalFE = rowNode.finiteElement();
auto const& colLocalFE = colNode.child(0).finiteElement();
auto const& quad = this->getQuadratureRule(context.type(), rowNode, colNode);
std::size_t rowSize = rowLocalFE.size();
std::size_t colSize = colLocalFE.size();
using RowLocalBasisType = typename std::decay_t<decltype(rowLocalFE)>::Traits::LocalBasisType;
using ColLocalBasisType = typename std::decay_t<decltype(colLocalFE)>::Traits::LocalBasisType;
......@@ -69,10 +71,10 @@ namespace AMDiS
auto const& rowShapeValues = rowShapeValuesCache[iq];
auto const& colShapeValues = colShapeValuesCache[iq];
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t j = 0; j < colSize; ++j) {
const auto value = b * (factor * rowShapeValues[i] * colShapeValues[j]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
......
src/amdis/assembler/ZeroOrderTestvec.hpp
View file @ c089668d
......@@ -45,8 +45,9 @@ namespace AMDiS
static const std::size_t CHILDREN = Node::CHILDREN;
auto const& localFE = node.child(0).finiteElement();
auto const& quad = this->getQuadratureRule(context.type(), node);
auto const& localFE = node.child(0).finiteElement();
std::size_t size = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType;
LocalBasisCache<LocalBasisType> cache(localFE.localBasis());
......@@ -62,7 +63,7 @@ namespace AMDiS
auto const& shapeValues = shapeValuesCache[iq];
for (std::size_t i = 0; i < localFE.size(); ++i) {
for (std::size_t i = 0; i < size; ++i) {
const auto value = exprValue * (factor * shapeValues[i]);
for (std::size_t k = 0; k < CHILDREN; ++k) {
const auto local_ki = node.child(k).localIndex(i);
......
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