Commit 061cbabb authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/dofvector_basis_constructor' into 'master'

Changed raw pointer to shared_ptr in DOFVector and DOFMatrix

See merge request !37
parents 3ca7eac6 aa77a167
...@@ -29,6 +29,6 @@ add_dependencies(examples ...@@ -29,6 +29,6 @@ add_dependencies(examples
cahn_hilliard.2d) cahn_hilliard.2d)
if (dune-foamgrid_FOUND) if (dune-foamgrid_FOUND)
add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3) add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3 ALBERTA_GRID)
add_dependencies(examples surface.2d) add_dependencies(examples surface.2d)
endif () endif ()
\ No newline at end of file
#include "config.h"
#include <amdis/AMDiS.hpp> #include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp> #include <amdis/Integrate.hpp>
#include <amdis/LocalOperators.hpp> #include <amdis/LocalOperators.hpp>
...@@ -14,7 +15,7 @@ using namespace AMDiS; ...@@ -14,7 +15,7 @@ using namespace AMDiS;
using namespace Dune::Indices; using namespace Dune::Indices;
using Grid = Dune::FoamGrid<2,3>; using Grid = Dune::FoamGrid<2,3>;
using Basis = LagrangeBasis<typename Grid::LeafGridView, 1>; using Basis = LagrangeBasis<Grid, 1>;
struct UnitRadius struct UnitRadius
{ {
...@@ -22,6 +23,9 @@ struct UnitRadius ...@@ -22,6 +23,9 @@ struct UnitRadius
static double eval(T&&) { return 1.0; } static double eval(T&&) { return 1.0; }
}; };
#if !HAVE_ALBERTA
#error "Need Alberta!!"
#endif
// solve the equation -laplace(u) - k^2 u = f on the sphere // solve the equation -laplace(u) - k^2 u = f on the sphere
int main(int argc, char** argv) int main(int argc, char** argv)
......
...@@ -35,7 +35,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R ...@@ -35,7 +35,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R
Dune::Hybrid::ifElse(std::is_same<RangeType_t<RN>, R>{}, Dune::Hybrid::ifElse(std::is_same<RangeType_t<RN>, R>{},
[&](auto id) { [&](auto id) {
auto subBasis = Dune::Functions::subspaceBasis(rhs.basis(), rowTreePath); auto subBasis = Dune::Functions::subspaceBasis(*rhs.basis(), rowTreePath);
Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_); Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_);
}); });
...@@ -45,7 +45,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R ...@@ -45,7 +45,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R
Dune::Hybrid::ifElse(std::is_same<RangeType_t<CN>, R>{}, Dune::Hybrid::ifElse(std::is_same<RangeType_t<CN>, R>{},
[&](auto id) { [&](auto id) {
auto subBasis = Dune::Functions::subspaceBasis(solution.basis(), colTreePath); auto subBasis = Dune::Functions::subspaceBasis(*solution.basis(), colTreePath);
Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_); Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_);
}); });
} }
......
...@@ -107,7 +107,7 @@ namespace AMDiS ...@@ -107,7 +107,7 @@ namespace AMDiS
protected: protected:
GridView const& gridView() const GridView const& gridView() const
{ {
return discreteFct_.basis().gridView(); return discreteFct_.basis()->gridView();
} }
private: private:
......
...@@ -241,7 +241,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R ...@@ -241,7 +241,7 @@ fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, R
rhs[a.second] += rhs[a.first]; rhs[a.second] += rhs[a.first];
solution[a.second] = solution[a.first]; solution[a.second] = solution[a.first];
} }
Dune::Functions::interpolate(rhs.basis(), rhs, [](auto const&) { return 0.0; }, periodicNodes_); Dune::Functions::interpolate(*rhs.basis(), rhs, [](auto const&) { return 0.0; }, periodicNodes_);
} }
} // end namespace AMDiS } // end namespace AMDiS
...@@ -45,7 +45,7 @@ void ProblemInstat<Traits>::createUhOld() ...@@ -45,7 +45,7 @@ void ProblemInstat<Traits>::createUhOld()
if (oldSolution_) if (oldSolution_)
warning("oldSolution already created\n"); warning("oldSolution already created\n");
else // create oldSolution else // create oldSolution
oldSolution_.reset(new SystemVector(*problemStat_->globalBasis(), DataTransferOperation::INTERPOLATE)); oldSolution_.reset(new SystemVector(problemStat_->globalBasis(), DataTransferOperation::INTERPOLATE));
} }
......
...@@ -174,9 +174,9 @@ void ProblemStat<Traits>::initGlobalBasis() ...@@ -174,9 +174,9 @@ void ProblemStat<Traits>::initGlobalBasis()
template <class Traits> template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors() void ProblemStat<Traits>::createMatricesAndVectors()
{ {
systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_); systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
solution_ = std::make_shared<SystemVector>(*globalBasis_, DataTransferOperation::INTERPOLATE); solution_ = std::make_shared<SystemVector>(globalBasis_, DataTransferOperation::INTERPOLATE);
rhs_ = std::make_shared<SystemVector>(*globalBasis_, DataTransferOperation::NO_OPERATION); rhs_ = std::make_shared<SystemVector>(globalBasis_, DataTransferOperation::NO_OPERATION);
auto localView = globalBasis_->localView(); auto localView = globalBasis_->localView();
for_each_node(localView.tree(), [&,this](auto const& node, auto treePath) for_each_node(localView.tree(), [&,this](auto const& node, auto treePath)
......
...@@ -31,11 +31,29 @@ namespace AMDiS ...@@ -31,11 +31,29 @@ namespace AMDiS
using type = remove_cvref_t<T>; using type = remove_cvref_t<T>;
}; };
template <class T>
struct UnderlyingType<T*>
{
using type = remove_cvref_t<T>;
};
template <class T> template <class T>
struct UnderlyingType<std::reference_wrapper<T>> struct UnderlyingType<std::reference_wrapper<T>>
{ {
using type = remove_cvref_t<T>; using type = remove_cvref_t<T>;
}; };
template <class T>
struct UnderlyingType<std::shared_ptr<T>>
{
using type = remove_cvref_t<T>;
};
template <class T>
struct UnderlyingType<std::unique_ptr<T>>
{
using type = remove_cvref_t<T>;
};
} }
template <class T> template <class T>
......
...@@ -25,11 +25,6 @@ namespace AMDiS ...@@ -25,11 +25,6 @@ namespace AMDiS
using TreePath = TP; using TreePath = TP;
public: public:
/// Constructor forwards to the treePath constructor, with empty TreePath
DOFVectorView(DOFVector<GB,VT>& dofVector)
: DOFVectorView{dofVector, Dune::TypeTree::hybridTreePath()}
{}
/// Constructor. Stores a pointer to the mutable `dofvector`. /// Constructor. Stores a pointer to the mutable `dofvector`.
template <class PreTreePath> template <class PreTreePath>
DOFVectorView(DOFVector<GB,VT>& dofVector, PreTreePath const& preTreePath) DOFVectorView(DOFVector<GB,VT>& dofVector, PreTreePath const& preTreePath)
...@@ -37,6 +32,11 @@ namespace AMDiS ...@@ -37,6 +32,11 @@ namespace AMDiS
, mutableDofVector_(&dofVector) , mutableDofVector_(&dofVector)
{} {}
/// Constructor forwards to the treePath constructor, with empty TreePath
DOFVectorView(DOFVector<GB,VT>& dofVector)
: DOFVectorView(dofVector, Dune::TypeTree::hybridTreePath())
{}
public: public:
/// \brief Interpolation of GridFunction to DOFVector, assuming that there is no /// \brief Interpolation of GridFunction to DOFVector, assuming that there is no
/// reference to this DOFVector in the expression. /// reference to this DOFVector in the expression.
...@@ -50,7 +50,7 @@ namespace AMDiS ...@@ -50,7 +50,7 @@ namespace AMDiS
template <class Expr, class Tag = tag::average> template <class Expr, class Tag = tag::average>
void interpolate_noalias(Expr&& expr, Tag strategy = {}) void interpolate_noalias(Expr&& expr, Tag strategy = {})
{ {
auto const& basis = this->basis(); auto const& basis = *this->basis();
auto const& treePath = this->treePath(); auto const& treePath = this->treePath();
auto&& gridFct = makeGridFunction(FWD(expr), basis.gridView()); auto&& gridFct = makeGridFunction(FWD(expr), basis.gridView());
......
...@@ -66,17 +66,17 @@ namespace AMDiS ...@@ -66,17 +66,17 @@ namespace AMDiS
class LocalFunction; class LocalFunction;
public: public:
/// Constructor forwards to the treePath constructor, with empty TreePath
DiscreteFunction(DOFVector<GB,VT> const& dofVector)
: DiscreteFunction{dofVector, Dune::TypeTree::hybridTreePath()}
{}
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath. /// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath) DiscreteFunction(DOFVector<GB,VT> const& dofVector, TP const& treePath)
: dofVector_(&dofVector) : dofVector_(&dofVector)
, treePath_(treePath) , treePath_(treePath)
, entitySet_(dofVector.basis().gridView()) , entitySet_(dofVector.basis()->gridView())
, nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(dofVector.basis(), treePath)) , nodeToRangeEntry_(Dune::Functions::makeDefaultNodeToRangeMap(*dofVector.basis(), treePath))
{}
/// Constructor forwards to the treePath constructor, with empty TreePath
DiscreteFunction(DOFVector<GB,VT> const& dofVector)
: DiscreteFunction(dofVector, Dune::TypeTree::hybridTreePath())
{} {}
/// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive /// \brief Evaluate DiscreteFunction in global coordinates. NOTE: expensive
...@@ -96,7 +96,7 @@ namespace AMDiS ...@@ -96,7 +96,7 @@ namespace AMDiS
public: public:
/// \brief Return global basis bound to the DOFVector /// \brief Return global basis bound to the DOFVector
GlobalBasis const& basis() const std::shared_ptr<GlobalBasis const> basis() const
{ {
return dofVector_->basis(); return dofVector_->basis();
} }
......
...@@ -23,15 +23,17 @@ private: ...@@ -23,15 +23,17 @@ private:
using Geometry = typename Element::Geometry; using Geometry = typename Element::Geometry;
public: public:
/// Constructor. Stores a copy of the DiscreteFunction.
LocalFunction(DiscreteFunction const& globalFunction) LocalFunction(DiscreteFunction const& globalFunction)
: globalFunction_(globalFunction) : globalFunction_(globalFunction)
, localView_(globalFunction_.basis().localView()) , localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath())) , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{} {}
/// Copy constructor.
LocalFunction(LocalFunction const& other) LocalFunction(LocalFunction const& other)
: globalFunction_(other.globalFunction_) : globalFunction_(other.globalFunction_)
, localView_(globalFunction_.basis().localView()) , localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath())) , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{} {}
...@@ -100,15 +102,17 @@ private: ...@@ -100,15 +102,17 @@ private:
using Geometry = typename Element::Geometry; using Geometry = typename Element::Geometry;
public: public:
/// Constructor. Stores a copy of the DiscreteFunction.
GradientLocalFunction(DiscreteFunction const& globalFunction) GradientLocalFunction(DiscreteFunction const& globalFunction)
: globalFunction_(globalFunction) : globalFunction_(globalFunction)
, localView_(globalFunction_.basis().localView()) , localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath())) , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{} {}
/// Copy constructor.
GradientLocalFunction(GradientLocalFunction const& other) GradientLocalFunction(GradientLocalFunction const& other)
: globalFunction_(other.globalFunction_) : globalFunction_(other.globalFunction_)
, localView_(globalFunction_.basis().localView()) , localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath())) , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{} {}
...@@ -260,7 +264,7 @@ operator()(Domain const& x) const ...@@ -260,7 +264,7 @@ operator()(Domain const& x) const
using Grid = typename GlobalBasis::GridView::Grid; using Grid = typename GlobalBasis::GridView::Grid;
using IS = typename GlobalBasis::GridView::IndexSet; using IS = typename GlobalBasis::GridView::IndexSet;
auto const& gv = this->basis().gridView(); auto const& gv = this->basis()->gridView();
Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()}; Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
auto element = hsearch.findEntity(x); auto element = hsearch.findEntity(x);
......
...@@ -38,31 +38,37 @@ namespace AMDiS ...@@ -38,31 +38,37 @@ namespace AMDiS
/// The index/size - type /// The index/size - type
using size_type = typename RowBasis::size_type; using size_type = typename RowBasis::size_type;
/// The type of the data matrix used in the backend
using BaseMatrix = typename Backend::BaseMatrix; using BaseMatrix = typename Backend::BaseMatrix;
/// The type of the matrix filled on an element with local contributions
using ElementMatrix = Dune::DynamicMatrix<double>; using ElementMatrix = Dune::DynamicMatrix<double>;
public: public:
/// Constructor. Constructs new BaseVector. /// Constructor. Stores the shared_ptr to the bases.
DOFMatrixBase(RowBasis const& rowBasis, ColBasis const& colBasis) DOFMatrixBase(std::shared_ptr<RowBasis> rowBasis, std::shared_ptr<ColBasis> colBasis)
: rowBasis_(&rowBasis) : rowBasis_(rowBasis)
, colBasis_(&colBasis) , colBasis_(colBasis)
{ {
operators_.init(rowBasis, colBasis); operators_.init(*rowBasis_, *colBasis_);
} }
DOFMatrixBase(DOFMatrixBase const&) = delete; /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into a new shared_ptr.
DOFMatrixBase(DOFMatrixBase&&) = delete; template <class RB_, class CB_>
DOFMatrixBase(RB_&& rowBasis, CB_&& colBasis)
: DOFMatrixBase(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
{}
/// Return the row-basis \ref rowBasis of the matrix /// Return the row-basis \ref rowBasis of the matrix
RowBasis const& rowBasis() const std::shared_ptr<RowBasis const> rowBasis() const
{ {
return *rowBasis_; return rowBasis_;
} }
/// Return the col-basis \ref colBasis of the matrix /// Return the col-basis \ref colBasis of the matrix
ColBasis const& colBasis() const std::shared_ptr<ColBasis const> colBasis() const
{ {
return *colBasis_; return colBasis_;
} }
/// Return the data-matrix /// Return the data-matrix
...@@ -97,6 +103,7 @@ namespace AMDiS ...@@ -97,6 +103,7 @@ namespace AMDiS
void finish(bool asmMatrix); void finish(bool asmMatrix);
/// Insert a block of values into the matrix (add to existing values) /// Insert a block of values into the matrix (add to existing values)
/// The global matrix indices are determined by the corresponding localviews.
void insert(RowLocalView const& rowLocalView, void insert(RowLocalView const& rowLocalView,
ColLocalView const& colLocalView, ColLocalView const& colLocalView,
ElementMatrix const& elementMatrix); ElementMatrix const& elementMatrix);
...@@ -126,6 +133,7 @@ namespace AMDiS ...@@ -126,6 +133,7 @@ namespace AMDiS
* [[expects: row is valid tree-path in RowBasis]] * [[expects: row is valid tree-path in RowBasis]]
* [[expects: col is valid tree-path in ColBasis]] * [[expects: col is valid tree-path in ColBasis]]
**/ **/
// TODO: add method without contextTag.
template <class ContextTag, class Expr, template <class ContextTag, class Expr,
class RowTreePath = RootTreePath, class ColTreePath = RootTreePath> class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
void addOperator(ContextTag contextTag, Expr const& expr, void addOperator(ContextTag contextTag, Expr const& expr,
...@@ -138,22 +146,18 @@ namespace AMDiS ...@@ -138,22 +146,18 @@ namespace AMDiS
/// Assemble all matrix operators, TODO: incooperate boundary conditions /// Assemble all matrix operators, TODO: incooperate boundary conditions
void assemble(); void assemble();
/// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds /// Number of nonzeros in the matrix
/// a one on the diagonal of the matrix.
auto applyDirichletBC(std::vector<char> const& dirichletNodes)
{
return backend_.applyDirichletBC(dirichletNodes);
}
size_type nnz() const size_type nnz() const
{ {
return backend_.nnz(); return backend_.nnz();
} }
protected: protected:
/// The finite element space / basis associated with the rows/columns /// The finite element space / basis associated with the rows
RowBasis const* rowBasis_; std::shared_ptr<RowBasis> rowBasis_;
ColBasis const* colBasis_;
/// The finite element space / basis associated with the columns
std::shared_ptr<ColBasis> colBasis_;
/// Data backend /// Data backend
Backend backend_; Backend backend_;
......
...@@ -10,12 +10,36 @@ ...@@ -10,12 +10,36 @@
#include <amdis/GridTransferManager.hpp> #include <amdis/GridTransferManager.hpp>
#include <amdis/OperatorList.hpp> #include <amdis/OperatorList.hpp>
#include <amdis/common/Math.hpp> #include <amdis/common/Math.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/linearalgebra/DOFVectorInterface.hpp> #include <amdis/linearalgebra/DOFVectorInterface.hpp>
#include <amdis/typetree/MultiIndex.hpp> #include <amdis/typetree/MultiIndex.hpp>
#include <amdis/typetree/TreePath.hpp> #include <amdis/typetree/TreePath.hpp>
namespace AMDiS namespace AMDiS
{ {
/// The container that stores a data-vector and a corresponding basis, should be
/// derived from \ref DOFVectorBase.
template <class GlobalBasis, class ValueType = double>
class DOFVector;
/// \brief Create a DOFVector from a basis.
/**
* This generator function accepts the basis as reference, temporary, or
* shared_ptr. Internally the reference is wrapped into a non-destroying
* shared_ptr and the temporary is moved into a new shared_ptr.
*
* The DataTransferOperation controls what is done during grid changes with the
* DOFVector. The default is interpolation of the data to the new grid. See
* \ref DataTransferOperation for more options.
**/
template <class ValueType = double, class GlobalBasis>
DOFVector<Underlying_t<GlobalBasis>, ValueType>
makeDOFVector(GlobalBasis&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
{
return {FWD(basis), op};
}
/// \brief The basic container that stores a base vector and a corresponding basis /// \brief The basic container that stores a base vector and a corresponding basis
/** /**
* Basis implementation of DOFVector, i.e. a vector storing all the * Basis implementation of DOFVector, i.e. a vector storing all the
...@@ -47,6 +71,8 @@ namespace AMDiS ...@@ -47,6 +71,8 @@ namespace AMDiS
/// The type of the data vector used in the backend /// The type of the data vector used in the backend
using BaseVector = typename Backend::BaseVector; using BaseVector = typename Backend::BaseVector;
/// The type of the vector filled on an element with local contributions
using ElementVector = Dune::DynamicVector<double>; using ElementVector = Dune::DynamicVector<double>;
/// Defines an interface to transfer the data during grid adaption /// Defines an interface to transfer the data during grid adaption
...@@ -56,16 +82,23 @@ namespace AMDiS ...@@ -56,16 +82,23 @@ namespace AMDiS
using DataTransferFactory = AMDiS::DataTransferFactory<Self>; using DataTransferFactory = AMDiS::DataTransferFactory<Self>;
public: public:
/// Constructor. Constructs new BaseVector. /// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
DOFVectorBase(GlobalBasis& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE) DOFVectorBase(std::shared_ptr<GlobalBasis> basis, DataTransferOperation op)
: basis_(&basis) : basis_(basis)
, dataTransfer_(DataTransferFactory::create(op, basis)) , dataTransfer_(DataTransferFactory::create(op, *basis_))
{ {
compress(); compress();
attachToGridTransfer(); attachToGridTransfer();
operators_.init(basis); operators_.init(*basis_);
} }
/// Constructor. Wraps the reference into a non-destroying shared_ptr or moves
/// the basis into a new shared_ptr.
template <class GB_>
DOFVectorBase(GB_&& basis, DataTransferOperation op)
: DOFVectorBase(Dune::wrap_or_move(FWD(basis)), op)
{}
/// Copy constructor /// Copy constructor
DOFVectorBase(Self const& that) DOFVectorBase(Self const& that)
: basis_(that.basis_) : basis_(that.basis_)
...@@ -111,7 +144,7 @@ namespace AMDiS ...@@ -111,7 +144,7 @@ namespace AMDiS
/// Sets each DOFVector to the scalar \p value. /// Sets each DOFVector to the scalar \p value.
template <class Number, template <class Number,
std::enable_if_t<Dune::IsNumber<Number>::value, int> = 0> REQUIRES(Dune::IsNumber<Number>::value)>
Self& operator=(Number value) Self& operator=(Number value)
{ {
backend_.set(value); backend_.set(value);
...@@ -119,9 +152,9 @@ namespace AMDiS ...@@ -119,9 +152,9 @@ namespace AMDiS
} }
/// Return the basis \ref basis_ associated to the vector /// Return the basis \ref basis_ associated to the vector
GlobalBasis const& basis() const std::shared_ptr<GlobalBasis const> basis() const
{ {
return *basis_; return basis_;
} }
/// Return the data-vector /// Return the data-vector
...@@ -142,6 +175,7 @@ namespace AMDiS ...@@ -142,6 +175,7 @@ namespace AMDiS
return basis_->dimension(); return basis_->dimension();
} }