Commit 601c4cb0 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Initial commit

Pipeline #1407 failed with stages
in 9 seconds
# We require version CMake version 3.1 to prevent issues
# with dune_enable_all_packages and older CMake versions.
cmake_minimum_required(VERSION 3.1)
project(dune-scalarbasis CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
#find dune-common and set the module path
find_package(dune-common REQUIRED)
#include the dune macros
# start a dune project with information from dune.module
# finalize the dune project, e.g. generating config.h etc.
Preparing the Sources
Additional to the software mentioned in README you'll need the
following programs installed on your system:
cmake >= 2.8.12
Getting started
If these preliminaries are met, you should run
dunecontrol all
which will find all installed dune modules as well as all dune modules
(not installed) which sources reside in a subdirectory of the current
directory. Note that if dune is not installed properly you will either
have to add the directory where the dunecontrol script resides (probably
./dune-common/bin) to your path or specify the relative path of the script.
Most probably you'll have to provide additional information to dunecontrol
(e. g. compilers, configure options) and/or make options.
The most convenient way is to use options files in this case. The files
define four variables:
CMAKE_FLAGS flags passed to cmake (during configure)
An example options file might look like this:
#use this options to configure and make if no other options are given
-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
If you save this information into example.opts you can pass the opts file to
dunecontrol via the --opts option, e. g.
dunecontrol --opts=example.opts all
More info
dunecontrol --help
for further options.
The full build system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!
set(modules "DuneScalarbasisMacros.cmake")
# File for module specific CMake tests.
/* begin dune-scalarbasis
put the definitions for config.h specific to
your project here. Everything above will be
/* begin private */
/* Name of package */
/* Define to the address where bug reports for this package should be sent. */
/* Define to the full name of this package. */
/* Define to the full name and version of this package. */
/* Define to the one symbol short name of this package. */
/* Define to the home page for this package. */
/* Define to the version of this package. */
/* end private */
/* Define to the version of dune-scalarbasis */
/* Define to the major version of dune-scalarbasis */
/* Define to the minor version of dune-scalarbasis */
/* Define to the revision of dune-scalarbasis */
/* end dune-scalarbasis
Everything below here will be overwritten
# shortcut for creating the and Doxyfile
# This file contains local changes to the doxygen configuration
# please us '+=' to add file/directories to the lists
# 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/
# see e.g. dune-grid for the examples of mainpage and 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/scalarbasis/test
# 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
# 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
# the \image command).
# IMAGE_PATH += @top_srcdir@/dune/scalarbasis/pics
Version: @VERSION@
Description: dune-scalarbasis module
Requires: dune-functions
Libs: -L${libdir}
Cflags: -I${includedir}
# Dune module information file #
#Name of the module
Module: dune-scalarbasis
Version: 0.1
#depending on
Depends: dune-functions
#install headers
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/functions/functionspacebases)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/lagrange/p0.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
namespace Dune {
namespace Functions {
// forward declarations
template<typename GV, typename TP>
class ScalarNode;
template<typename GV, class MI, class TP>
class ScalarNodeIndexSet;
template<typename GV, class MI>
class ScalarPreBasis;
template<typename GV, class MI>
class ScalarPreBasis
static const int dim = GV::dimension;
template<typename, class, class>
friend class ScalarNodeIndexSet;
template<class TP>
using Node = ScalarNode<GV, TP>;
template<class TP>
using IndexSet = ScalarNodeIndexSet<GV, MI, TP>;
using GridView = GV;
using size_type = std::size_t;
using MultiIndex = MI;
using SizePrefix = Dune::ReservedVector<size_type, 1>;
ScalarPreBasis(const GridView& gv)
: gridView_(gv)
/// \brief Initialize the global indices. No initialization necessary, since
/// there is only 1 global index.
void initializeIndices() {}
/// \brief Obtain the grid view that the basis is defined on
const GridView& gridView() const
return gridView_;
/// \brief Update the stored grid view, to be called if the grid has changed
void update (const GridView& gv)
gridView_ = gv;
* \brief Create tree node with given root tree path
* \tparam TP Type of root tree path
* \param tp Root tree path
template<class TP>
Node<TP> node(const TP& tp) const
return Node<TP>{tp};
* \brief Create tree node index set with given root tree path
* \tparam TP Type of root tree path
* \param tp Root tree path
template<class TP>
IndexSet<TP> indexSet() const
return IndexSet<TP>{};
//! Same as size(prefix) with empty prefix
size_type size() const
return 1;
//! Return number of possible values for next position in multi index
size_type size(const SizePrefix prefix) const
assert(prefix.size() == 0 || prefix.size() == 1);
return (prefix.size() == 0) ? size() : 0;
/// \brief Get the total dimension of the space spanned by this basis
size_type dimension() const
return size();
/// \brief Get the maximal number of DOFs associated to node for any element
size_type maxNodeSize() const
return 1;
GridView gridView_;
template<typename GV, typename TP>
class ScalarNode
: public LeafBasisNode<std::size_t, TP>
using Base = LeafBasisNode<std::size_t,TP>;
// required typedefs
using size_type = std::size_t;
using TreePath = TP;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = P0LocalFiniteElement<typename GV::ctype,double,GV::dimension>;
ScalarNode(const TreePath& treePath)
: Base(treePath)
/// \brief Return current element, assert if unbound
const Element& element() const
return *element_;
/// \brief Return the current bound local finite-element
const FiniteElement& finiteElement() const
return *finiteElement_;
/// \brief Bind to element `e`. Stores a pointer to the element and
/// creates a local finite-element.
void bind(const Element& e)
element_ = &e;
Std::optional<FiniteElement> finiteElement_;
const Element* element_ = nullptr;
template<typename GV, class MI, class TP>
class ScalarNodeIndexSet
static const int dim = GV::dimension;
// required typedefs
using size_type = std::size_t;
using MultiIndex = MI;
using PreBasis = ScalarPreBasis<GV, MI>;
using Node = typename PreBasis::template Node<TP>;
void bind(const Node& node) {}
void unbind() {}
/// \brief Size of subtree rooted in this node (element-local)
size_type size() const
return 1;
/// \brief The scalar basis represents a 1 dimensional space, so there is
/// only one index and that is 0
template<typename It>
It indices(It it) const
*it = {{ 0u }};
return it;
namespace BasisFactory {
namespace Imp {
class ScalarPreBasisFactory
static const std::size_t requiredMultiIndexSize = 1;
template<class MultiIndex, class GridView>
auto makePreBasis(const GridView& gridView) const
return ScalarPreBasis<GridView, MultiIndex>(gridView);
} // end namespace BasisFactory::Imp
* \brief Create a pre-basis factory that can create a scalar pre-basis.
* See \ref ScalarBasis for details about a scalar basis.
* \ingroup FunctionSpaceBasesImplementations
inline auto scalar()
return Imp::ScalarPreBasisFactory();
} // end namespace BasisFactory
* \brief A scalar basis represents a 1 dimensional function space `space{1}`
* of constants. Although bound to a GridView `GV` it is essential independent
* of any grid elements. On all bound elements the global index of any local DOF
* is the index 0.
* A ScalarBasis can be created with the factory \ref scalar().
* See the following classes for details:
* - \ref ScalarPreBasis
* - \ref ScalarNodeIndexSet
* - \ref ScalarNode
template<typename GV>
using ScalarBasis = DefaultGlobalBasis<ScalarPreBasis<GV, FlatMultiIndex<std::size_t>> >;
} // end namespace Functions
} // end namespace Dune
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <array>
#include <vector>
#include <dune/common/function.hh>
#include <dune/common/indices.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/functionspacebases/scalarbasis.hh>
* This example implements the stokes equation using an average integral
* constraint to eliminate the pressure nullspace:
* ```
* -\Delta u - \nabla p = 0 in domain, \\
* \nabla\cdot u = 0 in domain, \\
* \int p = 0, \\
* u = g on boundary
* ```
using namespace Dune;
template <class LocalView, class ElementMatrix>
void getLocalMatrix(const LocalView& localView, ElementMatrix& elementMatrix)
using Element = typename LocalView::Element;
const Element element = localView.element();
const int dim = Element::dimension;
auto geometry = element.geometry();
elementMatrix = 0;
using namespace Indices;
const auto& vLocalFE = localView.tree().child(_0,0).finiteElement();
const auto& pLocalFE = localView.tree().child(_1).finiteElement();
const auto& sLocalFE = localView.tree().child(_2).finiteElement();
int order = 2*(dim*vLocalFE.localBasis().order()-1);
const auto& quad = QuadratureRules<double, dim>::rule(element.type(), order);
for (const auto& qp : quad)
const auto JinvT = geometry.jacobianInverseTransposed(qp.position());
const auto dx = geometry.integrationElement(qp.position()) * qp.weight();
thread_local std::vector<FieldMatrix<double,1,dim> > refGrads;
vLocalFE.localBasis().evaluateJacobian(qp.position(), refGrads);
thread_local std::vector<FieldVector<double,dim> > gradients;
for (std::size_t i=0; i<gradients.size(); i++)[i][0], gradients[i]);
// laplace(u)
for (std::size_t i=0; i<vLocalFE.size(); i++)
for (std::size_t j=0; j<vLocalFE.size(); j++)
for (std::size_t k=0; k<dim; k++)
std::size_t row = localView.tree().child(_0,k).localIndex(i);
std::size_t col = localView.tree().child(_0,k).localIndex(j);
elementMatrix[row][col] += (gradients[i] * gradients[j]) * dx;
thread_local std::vector<FieldVector<double,1> > pValues;
pLocalFE.localBasis().evaluateFunction(qp.position(), pValues);
// grad(p) and div(u)
for (std::size_t i=0; i<vLocalFE.size(); i++)
for (std::size_t j=0; j<pLocalFE.size(); j++)
for (std::size_t k=0; k<dim; k++)
std::size_t vIndex = localView.tree().child(_0,k).localIndex(i);
std::size_t pIndex = localView.tree().child(_1).localIndex(j);
elementMatrix[vIndex][pIndex] += gradients[i][k] * pValues[j] * dx;
elementMatrix[pIndex][vIndex] += gradients[i][k] * pValues[j] * dx;
// int(p) == 0 condition
for (std::size_t i=0; i<pLocalFE.size(); i++)
for (std::size_t j=0; j<sLocalFE.size(); j++)
std::size_t pIndex = localView.tree().child(_1).localIndex(i);
std::size_t sIndex = localView.tree().child(_2).localIndex(j);
elementMatrix[sIndex][pIndex] += pValues[i] * dx;
elementMatrix[pIndex][sIndex] += pValues[i] * dx;
template <class Basis, class MatrixType>
void setOccupationPattern(const Basis& basis, MatrixType& matrix)
MatrixIndexSet nb(basis.dimension(), basis.dimension());
auto localView = basis.localView();
for(const auto& element : elements(basis.gridView()))
for (std::size_t i=0; i<localView.size(); i++) {
auto row = localView.index(i);
for (std::size_t j=0; j<localView.size(); j++) {
auto col = localView.index(j);