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

optimized assembler, i.e. ommit assembling of empty blocks

parent c4446c8c
Hide whitespace changes
Inline Side-by-side
doc/doxygen/Doxylocal
View file @ 5d5e8c7d
# This file contains local changes to the doxygen configuration
# please us '+=' to add file/directories to the lists
FILE_PATTERNS += *.hpp \
*.cpp
HIDE_SCOPE_NAMES = YES
HIDE_UNDOC_CLASSES = NO
INTERNAL_DOCS = NO
MARKDOWN_SUPPORT = YES
EXCLUDE_SYMBOLS = AMDiS::Impl \
AMDiS::traits::Impl \
AMDiS::detail \
itl::details
PREDEFINED += HAVE_UMFPACK \
HAVE_ALBERTA \
HAVE_UG \
AMDIS_BACKEND_MTL
# The INPUT tag can be used to specify the files and/or directories that contain
# documented source files. You may enter file names like "myfile.cpp" or
# directories like "/usr/src/myproject". Separate the files or directories
# with spaces.
INPUT += @top_srcdir@/dune/
INPUT += @top_srcdir@/dune/amdis \
@top_srcdir@/dune/amdis/common \
@top_srcdir@/dune/amdis/utility \
@top_srcdir@/dune/amdis/linear_algebra \
@top_srcdir@/dune/amdis/linear_algebra/mtl
# see e.g. dune-grid for the examples of mainpage and modules
# INPUT += @srcdir@/mainpage \
# @srcdir@/modules
#INPUT += @srcdir@/mainpage \
# @srcdir@/modules
# The EXCLUDE tag can be used to specify files and/or directories that should
# excluded from the INPUT source files. This way you can easily exclude a
# subdirectory from a directory tree whose root is specified with the INPUT tag.
# EXCLUDE += @top_srcdir@/dune/amdis/test
EXCLUDE += @top_srcdir@/dune/amdis/test \
@top_srcdir@/dune/amdis/linear_algebra/istl
# The EXAMPLE_PATH tag can be used to specify one or more files or
# directories that contain example code fragments that are included (see
# the \include command).
# EXAMPLE_PATH += @top_srcdir@/src
EXAMPLE_PATH += @top_srcdir@/src
# The IMAGE_PATH tag can be used to specify one or more files or
# directories that contain image that are included in the documentation (see
......
dune-amdis.pc.in
View file @ 5d5e8c7d
......@@ -9,7 +9,7 @@ DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: amdis module
URL: http://dune-project.org/
URL: http://www.github.com/spraetor
Requires: dune-common dune-geometry dune-localfunctions dune-istl dune-typetree dune-grid dune-functions
Libs: -L${libdir}
Cflags: -I${includedir}
dune/amdis/DirichletBC.hpp
View file @ 5d5e8c7d
......@@ -10,13 +10,28 @@
namespace AMDiS
{
/// Implements a boundary condition of Dirichlet-type.
/**
* By calling the methods \ref init() and \ref finish before and after
* assembling the system-matrix, respectively, dirichlet boundary conditions
* can be applied to the matrix and system vector. Therefore a predicate
* functions indicates the DOFs where values should be enforced and a second
* functor provided in the constructor is responsible for determining the
* values to be set at the DOFs.
*
* In the \ref finish method the matrix is called with \ref applyDirichletBC
* to erase the corresponding rows and columns for the DOF indices. This
* application of boundary conditions can be symmetric if the matrix does
* support this symmetric modification. As a result, this method returns a list
* of columns values, that should be subtracted from the rhs.
**/
template <class WorldVector>
class DirichletBC
{
public:
template <class Predicate, class Values,
class = std::enable_if_t< concepts::Functor<Predicate, bool(WorldVector)> &&
concepts::Functor<Values, double(WorldVector)> > >
class = std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> &&
Concepts::Functor<Values, double(WorldVector)> > >
DirichletBC(Predicate&& predicate, Values&& values)
: predicate(std::forward<Predicate>(predicate))
, values(std::forward<Values>(values))
......
dune/amdis/Math.hpp
View file @ 5d5e8c7d
/** \defgroup Common Common
*/
#pragma once
#include <string>
......@@ -8,26 +5,40 @@
#include <cmath>
#include <cfloat>
#include <boost/math/special_functions/pow.hpp>
#include <dune/amdis/common/ScalarTypes.hpp>
namespace AMDiS
{
namespace math
namespace Math
{
template <class T>
std::enable_if_t<std::is_arithmetic<T>::value, T>
constexpr abs(T a)
/// Implementation of the absolute value \f$|a|\f$ of arithmetic types.
template <class T,
class = std::enable_if_t<Concepts::Arithmetic<T>> >
constexpr T abs(T a)
{
return a >= 0 ? a : -a;
}
template <class T>
std::enable_if_t<std::is_arithmetic<T>::value, T>
constexpr sqr(T a)
/// Implementation of the square \f$ a^2 \f$ of arithmetic types.
template <class T,
class = std::enable_if_t<Concepts::Arithmetic<T>> >
constexpr auto sqr(T a)
{
return a*a;
}
template <size_t p, class T,
class = std::enable_if_t<Concepts::Arithmetic<T>> >
constexpr auto pow(T v)
{
return boost::math::pow<p>(v);
}
/// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type
/// supporting the `>` relation.
template <class T0, class T1>
constexpr auto min(T0 a, T1 b)
{
......@@ -35,13 +46,15 @@ namespace AMDiS
}
/// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type
/// supporting the `>` relation.
template <class T0, class T1>
constexpr auto max(T0 a, T1 b)
{
return a > b ? a : b;
}
} // end namespace math
} // end namespace Math
template <class T> inline void nullify(T& a)
......
dune/amdis/Mesh.hpp
View file @ 5d5e8c7d
......@@ -129,6 +129,9 @@ namespace AMDiS
} else {
AMDIS_ERROR_EXIT("Construction of UGGrid without filename not yet implemented!");
}
AMDIS_ERROR_EXIT("No way to construct UG-Grid found");
return {};
}
};
#endif
......
dune/amdis/Operator.hpp
View file @ 5d5e8c7d
#pragma once
#include <list>
#include <vector>
#include "OperatorTerm.hpp"
#include "TermGenerator.hpp"
#include <dune/amdis/OperatorTermBase.hpp>
#include <dune/amdis/terms/TermGenerator.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
#include <dune/amdis/utility/GetDegree.hpp>
namespace AMDiS
{
......@@ -19,6 +23,11 @@ namespace AMDiS
{
using Self = Operator;
using OperatorTermType = OperatorTerm<MeshView>;
using ElementVector = Impl::ElementVector;
using ElementMatrix = Impl::ElementMatrix;
using IdType = typename MeshView::Grid::LocalIdSet::IdType;
public:
/// Add coefficients for zero-order operator < c * u, v >.
......@@ -115,52 +124,46 @@ namespace AMDiS
int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double* factor = NULL);
template <class RowView, class ElementVector>
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor = NULL);
protected:
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor);
ElementMatrix& elementMatrix);
template <class RowView, class ElementVector>
template <class RowView>
void assembleZeroOrderTerms(RowView const& rowView,
ElementVector& elementVector,
double factor);
ElementVector& elementVector);
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor);
ElementMatrix& elementMatrix);
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor);
ElementMatrix& elementMatrix);
template <class RowView, class ElementVector>
template <class RowView>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ElementVector& elementVector,
double factor);
ElementVector& elementVector);
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor);
ElementMatrix& elementMatrix);
private:
template <class Term>
......@@ -179,6 +182,7 @@ namespace AMDiS
template <class Term, size_t I, size_t J>
Self& addSOTImpl(Term const& term, const index_<I>, const index_<J>);
template <class... Args>
static shared_ptr<Self> create(index_<0>, Args&&... args)
{
......@@ -197,6 +201,9 @@ namespace AMDiS
return sot(std::forward<Args>(args)...);
}
private:
template <class LocalView>
IdType getElementId(LocalView const& localView);
private:
/// List of all second order terms
......@@ -213,6 +220,12 @@ namespace AMDiS
int psiDegree = 1;
int phiDegree = 1;
IdType lastMatrixId = 0;
IdType lastVectorId = 0;
ElementMatrix cachedElementMatrix;
ElementVector cachedElementVector;
};
} // end namespace AMDiS
......
dune/amdis/Operator.inc.hpp
View file @ 5d5e8c7d
#pragma once
#include <vector>
namespace AMDiS
{
template <class MeshView>
......@@ -11,14 +9,18 @@ namespace AMDiS
{
AMDIS_FUNCNAME("Operator::init()");
// const auto& rowFE = rowView.tree().finiteElement();
// const auto& colFE = colView.tree().finiteElement();
using IdType = typename Operator<MeshView>::IdType;
lastMatrixId = std::numeric_limits<IdType>::max();
lastVectorId = std::numeric_limits<IdType>::max();
// auto const& rowFE = rowView.tree().finiteElement();
// auto const& colFE = colView.tree().finiteElement();
// psiDegree = rowFE.localBasis().order();
// phiDegree = colFE.localBasis().order();
psiDegree = GetDegree<RowFeSpace>::value;
phiDegree = GetDegree<ColFeSpace>::value;
psiDegree = getPolynomialDegree<RowFeSpace>;
phiDegree = getPolynomialDegree<ColFeSpace>;
// TODO: calc quadrature degree here.
}
......@@ -52,9 +54,22 @@ namespace AMDiS
return psiDegree + phiDegree - order + maxTermDegree;
}
template <class MeshView>
template <class LocalView>
typename Operator<MeshView>::IdType
Operator<MeshView>::getElementId(LocalView const& localView)
{
auto const& element = localView.element();
auto const& grid = localView.globalBasis().gridView().grid();
auto const& idSet = grid.localIdSet();
return idSet.id(element);
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
......@@ -67,21 +82,44 @@ namespace AMDiS
if (fac == 0.0)
return false;
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, colView, elementMatrix, fac);
if (!firstOrderGrdPhi.empty())
assembleFirstOrderTermsGrdPhi(rowView, colView, elementMatrix, fac);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, colView, elementMatrix, fac);
if (!secondOrder.empty())
assembleSecondOrderTerms(rowView, colView, elementMatrix, fac);
const auto nRows = rowView.tree().finiteElement().size();
const auto nCols = colView.tree().finiteElement().size();
auto currentId = getElementId(rowView);
bool useCachedMatrix = (lastMatrixId == currentId
&& num_rows(cachedElementMatrix) == nRows
&& num_cols(cachedElementMatrix) == nCols);
bool assign = true;
if (!useCachedMatrix) {
cachedElementMatrix.change_dim(nRows, nCols);
set_to_zero(cachedElementMatrix);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, colView, cachedElementMatrix);
if (!firstOrderGrdPhi.empty())
assembleFirstOrderTermsGrdPhi(rowView, colView, cachedElementMatrix);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, colView, cachedElementMatrix);
if (!secondOrder.empty())
assembleSecondOrderTerms(rowView, colView, cachedElementMatrix);
assign = !zeroOrder.empty() || !firstOrderGrdPhi.empty() ||
!firstOrderGrdPsi.empty() || !secondOrder.empty();
}
if (assign)
elementMatrix += fac * cachedElementMatrix;
lastMatrixId = currentId;
return true;
}
template <class MeshView>
template <class RowView, class ElementVector>
template <class RowView>
bool Operator<MeshView>::getElementVector(RowView const& rowView,
ElementVector& elementVector,
double* factor)
......@@ -96,21 +134,38 @@ namespace AMDiS
AMDIS_TEST_EXIT( firstOrderGrdPhi.empty() && secondOrder.empty(),
"Derivatives on ansatz-functions not allowed on the vector-side!");
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, elementVector, fac);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, elementVector, fac);
const auto nRows = rowView.tree().finiteElement().size();
auto currentId = getElementId(rowView);
bool useCachedVector = (lastVectorId == currentId && size(cachedElementVector) == nRows);
bool assign = true;
if (!useCachedVector) {
cachedElementVector.change_dim(nRows);
set_to_zero(cachedElementVector);
if (!zeroOrder.empty())
assembleZeroOrderTerms(rowView, cachedElementVector);
if (!firstOrderGrdPsi.empty())
assembleFirstOrderTermsGrdPsi(rowView, cachedElementVector);
assign = !zeroOrder.empty() || !firstOrderGrdPsi.empty();
}
if (assign)
elementVector += fac * cachedElementVector;
lastVectorId = currentId;
return true;
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor)
ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");
......@@ -119,11 +174,11 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
......@@ -133,7 +188,7 @@ namespace AMDiS
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos) * factor;
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
......@@ -143,14 +198,13 @@ namespace AMDiS
else
colLocalBasis.evaluateFunction(quadPos, colShapeValues);
for (size_t i = 0; i < elementMatrix.N(); ++i) {
for (size_t j = 0; j < elementMatrix.M(); ++j) {
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : zeroOrder)
elementMatrix[local_i][local_j]
+= operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j])
* quad[iq].weight() * integrationElement;
+= operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j]) * factor;
}
}
}
......@@ -159,10 +213,9 @@ namespace AMDiS
template <class MeshView>
template <class RowView, class ElementVector>
template <class RowView>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
ElementVector& elementvector,
double factor)
ElementVector& elementvector)
{
AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");
......@@ -171,41 +224,40 @@ namespace AMDiS
const int dim = Element::dimension;
auto geometry = element.geometry();
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(0);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos) * factor;
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
for (size_t i = 0; i < elementvector.size(); ++i) {
for (size_t i = 0; i < size(elementvector); ++i) {
const int local_i = rowView.tree().localIndex(i);
for (auto* operatorTerm : zeroOrder)
elementvector[local_i]
+= operatorTerm->evalZot(iq, rowShapeValues[i])
* quad[iq].weight() * integrationElement;
elementvector[local_i] += operatorTerm->evalZot(iq, rowShapeValues[i]) * factor;
}
}
}
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
double factor)
ElementMatrix& elementMatrix)
{
AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");
......@@ -214,15 +266,17 @@ namespace AMDiS
auto geometry = element.geometry();
const int dim = Element::dimension;
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
auto const& rowLocalBasis = rowView.tree().finiteElement().localBasis();
auto const& colLocalBasis = colView.tree().finiteElement().localBasis();
int order = getQuadratureDegree(2);
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
auto const& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPhi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -231,7 +285,7 @@ namespace AMDiS
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
// The multiplicative factor in the integral transformation formula
const double integrationElement = geometry.integrationElement(quadPos) * factor;
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
......@@ -246,14 +300,13 @@ namespace AMDiS
for (size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
for (size_t i = 0; i < elementMatrix.N(); ++i) {
for (size_t j = 0; j < elementMatrix.M(); ++j) {
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : firstOrderGrdPhi)
elementMatrix[local_i][local_j]
+= operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j])
* quad[iq].weight() * integrationElement;
+= operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j]) * factor;
}
}
}
......@@ -261,11 +314,10 @@ namespace AMDiS
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
template <class RowView, class ColView>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,