Commit 7668c6b7 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

a bit cleanup, but nevertheless, valgrind gives lots of errors. Needs corrections.

parent e18609d3
cmake_minimum_required(VERSION 3.1) cmake_minimum_required(VERSION 3.1)
project(dune-amdis CXX) project(dune-amdis CXX)
set(CXX_MAX_STANDARD 14 CACHE BOOL "" FORCE)
if(NOT (dune-common_DIR OR dune-common_ROOT OR if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*")) "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
${PROJECT_BINARY_DIR}) ${PROJECT_BINARY_DIR})
endif() endif()
set(CXX_MAX_STANDARD 14)
#find dune-common and set the module path #find dune-common and set the module path
find_package(dune-common REQUIRED) find_package(dune-common REQUIRED)
......
...@@ -17,9 +17,11 @@ namespace AMDiS ...@@ -17,9 +17,11 @@ namespace AMDiS
const Flag EQUAL_BASES = 0x01L; const Flag EQUAL_BASES = 0x01L;
const Flag EQUAL_LOCALBASES = 0x02L; const Flag EQUAL_LOCALBASES = 0x02L;
template <class GlobalBasis> template <class Traits>
class Assembler class Assembler
{ {
using GlobalBasis = typename Traits::GlobalBasis;
/// The grid view the global FE basis lives on /// The grid view the global FE basis lives on
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
......
...@@ -6,9 +6,9 @@ ...@@ -6,9 +6,9 @@
namespace AMDiS { namespace AMDiS {
template <class GlobalBasis> template <class Traits>
template <class SystemMatrixType, class SystemVectorType> template <class SystemMatrixType, class SystemVectorType>
void Assembler<GlobalBasis>::assemble( void Assembler<Traits>::assemble(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
...@@ -91,9 +91,9 @@ void Assembler<GlobalBasis>::assemble( ...@@ -91,9 +91,9 @@ void Assembler<GlobalBasis>::assemble(
} }
template <class GlobalBasis> template <class Traits>
template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews> template <class ElementContainer, class Container, class Operators, class Geometry, class... LocalViews>
void Assembler<GlobalBasis>::assembleElementOperators( void Assembler<Traits>::assembleElementOperators(
ElementContainer& elementContainer, ElementContainer& elementContainer,
Container& container, Container& container,
Operators& operators, Operators& operators,
...@@ -131,9 +131,9 @@ void Assembler<GlobalBasis>::assembleElementOperators( ...@@ -131,9 +131,9 @@ void Assembler<GlobalBasis>::assembleElementOperators(
} }
template <class GlobalBasis> template <class Traits>
template <class SystemMatrixType, class SystemVectorType> template <class SystemMatrixType, class SystemVectorType>
void Assembler<GlobalBasis>::initMatrixVector( void Assembler<Traits>::initMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
...@@ -147,7 +147,7 @@ void Assembler<GlobalBasis>::initMatrixVector( ...@@ -147,7 +147,7 @@ void Assembler<GlobalBasis>::initMatrixVector(
forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath) forEachNode(localView.tree(), [&,this](auto const& rowNode, auto rowTreePath)
{ {
if (rowNode.isLeaf) if (rowNode.isLeaf)
msg(0, " DOFs for Basis[", to_string(rowTreePath), "]"); // TODO: add right values msg(globalBasis_.dimension(rowTreePath), " DOFs for Basis[", to_string(rowTreePath), "]");
auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath); auto rowBasis = Dune::Functions::subspaceBasis(globalBasis_, rowTreePath);
//if (rhsOperators_[rowNode].assemble(asmVector)) //if (rhsOperators_[rowNode].assemble(asmVector))
...@@ -169,9 +169,9 @@ void Assembler<GlobalBasis>::initMatrixVector( ...@@ -169,9 +169,9 @@ void Assembler<GlobalBasis>::initMatrixVector(
} }
template <class GlobalBasis> template <class Traits>
template <class SystemMatrixType, class SystemVectorType> template <class SystemMatrixType, class SystemVectorType>
std::size_t Assembler<GlobalBasis>::finishMatrixVector( std::size_t Assembler<Traits>::finishMatrixVector(
SystemMatrixType& matrix, SystemMatrixType& matrix,
SystemVectorType& solution, SystemVectorType& solution,
SystemVectorType& rhs, SystemVectorType& rhs,
......
...@@ -17,12 +17,14 @@ ...@@ -17,12 +17,14 @@
namespace AMDiS namespace AMDiS
{ {
template <class GlobalBasis, class TreePath> template <class Traits, class TreePath>
class FileWriter class FileWriter
: public FileWriterInterface : public FileWriterInterface
{ {
private: // typedefs and static constants private: // typedefs and static constants
using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
/// Dimension of the mesh /// Dimension of the mesh
...@@ -200,4 +202,15 @@ namespace AMDiS ...@@ -200,4 +202,15 @@ namespace AMDiS
int mode_; int mode_;
}; };
template <class Traits, class GlobalBasis, class TreePath, class Vector>
std::shared_ptr<FileWriter<Traits,TreePath>> makeFileWriterPtr(
std::string baseName,
std::shared_ptr<GlobalBasis> const& basis,
TreePath const& tp,
std::shared_ptr<Vector> const& vector)
{
return std::make_shared<FileWriter<Traits,TreePath>>(baseName, basis, tp, vector);
}
} // end namespace AMDiS } // end namespace AMDiS
...@@ -42,7 +42,7 @@ ...@@ -42,7 +42,7 @@
namespace AMDiS namespace AMDiS
{ {
template <class GlobalBasis> template <class Traits>
class ProblemStat class ProblemStat
: public ProblemStatBase : public ProblemStatBase
, public StandardProblemIteration , public StandardProblemIteration
...@@ -51,12 +51,13 @@ namespace AMDiS ...@@ -51,12 +51,13 @@ namespace AMDiS
public: // typedefs and static constants public: // typedefs and static constants
using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView; using GridView = typename GlobalBasis::GridView;
using Grid = typename GridView::Grid; using Grid = typename GridView::Grid;
using Element = typename GridView::template Codim<0>::Entity; using Element = typename GridView::template Codim<0>::Entity;
static const int nComponents = 1; // TODO: count leaf nodes in GlobalBasis static const std::size_t nComponents = 1; // TODO: count leaf nodes in GlobalBasis
/// Dimension of the mesh /// Dimension of the mesh
static constexpr int dim = Grid::dimension; static constexpr int dim = Grid::dimension;
...@@ -309,7 +310,7 @@ namespace AMDiS ...@@ -309,7 +310,7 @@ namespace AMDiS
return componentNames; return componentNames;
} }
int getNumComponents() const std::size_t getNumComponents() const
{ {
return nComponents; return nComponents;
} }
......
...@@ -10,8 +10,8 @@ ...@@ -10,8 +10,8 @@
namespace AMDiS { namespace AMDiS {
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>::initialize( void ProblemStat<Traits>::initialize(
Flag initFlag, Flag initFlag,
Self* adoptProblem, Self* adoptProblem,
Flag adoptFlag) Flag adoptFlag)
...@@ -105,8 +105,8 @@ void ProblemStat<GlobalBasis>::initialize( ...@@ -105,8 +105,8 @@ void ProblemStat<GlobalBasis>::initialize(
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>::createFileWriter() void ProblemStat<Traits>::createFileWriter()
{ {
auto localView = globalBasis->localView(); auto localView = globalBasis->localView();
forEachNode(localView.tree(), [&,this](auto const& node, auto treePath) forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
...@@ -116,8 +116,7 @@ void ProblemStat<GlobalBasis>::createFileWriter() ...@@ -116,8 +116,7 @@ void ProblemStat<GlobalBasis>::createFileWriter()
if (!Parameters::get<std::string>(componentName + "->filename")) if (!Parameters::get<std::string>(componentName + "->filename"))
return; return;
using TP = decltype(treePath); auto writer = makeFileWriterPtr<Traits>(componentName, globalBasis, treePath, solution);
auto writer = std::make_shared<FileWriter<GlobalBasis,TP>>(componentName, globalBasis, treePath, solution);
filewriter.push_back(writer); filewriter.push_back(writer);
}); });
} }
...@@ -125,9 +124,9 @@ void ProblemStat<GlobalBasis>::createFileWriter() ...@@ -125,9 +124,9 @@ void ProblemStat<GlobalBasis>::createFileWriter()
// add matrix/vector operator terms // add matrix/vector operator terms
template <class GlobalBasis> template <class Traits>
template <class RowTreePath, class ColTreePath> template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator( void ProblemStat<Traits>::addMatrixOperator(
ElementOperator const& op, ElementOperator const& op,
RowTreePath row, ColTreePath col, RowTreePath row, ColTreePath col,
double* factor, double* estFactor) double* factor, double* estFactor)
...@@ -139,9 +138,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( ...@@ -139,9 +138,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator(
} }
template <class GlobalBasis> template <class Traits>
template <class RowTreePath, class ColTreePath> template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator( void ProblemStat<Traits>::addMatrixOperator(
IntersectionOperator const& op, IntersectionOperator const& op,
RowTreePath row, ColTreePath col, RowTreePath row, ColTreePath col,
double* factor, double* estFactor) double* factor, double* estFactor)
...@@ -153,9 +152,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( ...@@ -153,9 +152,9 @@ void ProblemStat<GlobalBasis>::addMatrixOperator(
} }
template <class GlobalBasis> template <class Traits>
template <class RowTreePath, class ColTreePath> template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator( void ProblemStat<Traits>::addMatrixOperator(
BoundaryType b, BoundaryType b,
IntersectionOperator const& op, IntersectionOperator const& op,
RowTreePath row, ColTreePath col, RowTreePath row, ColTreePath col,
...@@ -168,8 +167,8 @@ void ProblemStat<GlobalBasis>::addMatrixOperator( ...@@ -168,8 +167,8 @@ void ProblemStat<GlobalBasis>::addMatrixOperator(
} }
#if 0 #if 0
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops) addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
{ {
for (auto op : ops){ for (auto op : ops){
...@@ -180,8 +179,8 @@ addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops) ...@@ -180,8 +179,8 @@ addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
} }
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops) addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
{ {
for (auto op : ops) { for (auto op : ops) {
...@@ -194,9 +193,9 @@ addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> ...@@ -194,9 +193,9 @@ addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool>
#endif #endif
template <class GlobalBasis> template <class Traits>
template <class TreePath> template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator( void ProblemStat<Traits>::addVectorOperator(
ElementOperator const& op, ElementOperator const& op,
TreePath path, TreePath path,
double* factor, double* estFactor) double* factor, double* estFactor)
...@@ -207,9 +206,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator( ...@@ -207,9 +206,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator(
} }
template <class GlobalBasis> template <class Traits>
template <class TreePath> template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator( void ProblemStat<Traits>::addVectorOperator(
IntersectionOperator const& op, IntersectionOperator const& op,
TreePath path, TreePath path,
double* factor, double* estFactor) double* factor, double* estFactor)
...@@ -220,9 +219,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator( ...@@ -220,9 +219,9 @@ void ProblemStat<GlobalBasis>::addVectorOperator(
} }
template <class GlobalBasis> template <class Traits>
template <class TreePath> template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator( void ProblemStat<Traits>::addVectorOperator(
BoundaryType b, BoundaryType b,
IntersectionOperator const& op, IntersectionOperator const& op,
TreePath path, TreePath path,
...@@ -234,8 +233,8 @@ void ProblemStat<GlobalBasis>::addVectorOperator( ...@@ -234,8 +233,8 @@ void ProblemStat<GlobalBasis>::addVectorOperator(
} }
#if 0 #if 0
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
addVectorOperator(std::map< int, ElementOperator> ops) addVectorOperator(std::map< int, ElementOperator> ops)
{ {
for (auto op : ops) { for (auto op : ops) {
...@@ -244,8 +243,8 @@ addVectorOperator(std::map< int, ElementOperator> ops) ...@@ -244,8 +243,8 @@ addVectorOperator(std::map< int, ElementOperator> ops)
} }
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops) addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
{ {
for (auto op : ops) { for (auto op : ops) {
...@@ -257,9 +256,9 @@ addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops) ...@@ -257,9 +256,9 @@ addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
// Adds a Dirichlet boundary condition // Adds a Dirichlet boundary condition
template <class GlobalBasis> template <class Traits>
template <class Predicate, class RowTreePath, class ColTreePath, class Values> template <class Predicate, class RowTreePath, class ColTreePath, class Values>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values) addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
{ {
static_assert( Concepts::Functor<Predicate, bool(WorldVector)>, static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
...@@ -279,8 +278,8 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val ...@@ -279,8 +278,8 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{ {
Timer t; Timer t;
...@@ -290,6 +289,11 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) ...@@ -290,6 +289,11 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solverInfo.setStoreMatrixData(storeMatrixData); solverInfo.setStoreMatrixData(storeMatrixData);
solution->compress(); solution->compress();
msg("matrix: ",num_rows(systemMatrix->getMatrix()), " x ", num_cols(systemMatrix->getMatrix()));
msg("x: ",num_rows(solution->getVector()));
msg("rhs: ",num_rows(rhs->getVector()));
linearSolver->solve(systemMatrix->getMatrix(), linearSolver->solve(systemMatrix->getMatrix(),
solution->getVector(), rhs->getVector(), solution->getVector(), rhs->getVector(),
solverInfo); solverInfo);
...@@ -319,13 +323,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData) ...@@ -319,13 +323,13 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector) buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
{ {
Timer t; Timer t;
Assembler<GlobalBasis> assembler(*globalBasis, matrixOperators, rhsOperators, constraints); Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
auto gv = leafGridView(); auto gv = leafGridView();
assembler.update(gv); assembler.update(gv);
...@@ -335,8 +339,8 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool ...@@ -335,8 +339,8 @@ buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool
} }
template <class GlobalBasis> template <class Traits>
void ProblemStat<GlobalBasis>:: void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force) writeFiles(AdaptInfo& adaptInfo, bool force)
{ {
Timer t; Timer t;
......
...@@ -80,18 +80,26 @@ namespace AMDiS ...@@ -80,18 +80,26 @@ namespace AMDiS
/// \brief A factory for a composite basis composed of lagrange bases of different degree. /// \brief A factory for a composite basis composed of lagrange bases of different degree.
template <class GridView, int... degrees> template <class GridView, int... degrees>
using LagrangeBasis struct LagrangeBasis {
= Dune::Functions::DefaultGlobalBasis<typename Impl::LagrangeBasisBuilder<GridView, degrees...>::type>; using GlobalBasis
= Dune::Functions::DefaultGlobalBasis<typename Impl::LagrangeBasisBuilder<GridView, degrees...>::type>;
};
/// Specialization of \ref LagrangeBasis for Grid type \ref Dune::YaspGrid for a given dimension. /// Specialization of \ref LagrangeBasis for Grid type \ref Dune::YaspGrid for a given dimension.
template <int dim, int... degrees> template <int dim, int... degrees>
using YaspGridBasis = LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...>; struct YaspGridBasis
{
using GlobalBasis = typename LagrangeBasis<typename Dune::YaspGrid<dim>::LeafGridView, degrees...>::GlobalBasis;
};
template <class GridView> template <class GridView>
using TaylorHoodBasis struct TaylorHoodBasis
= Dune::Functions::DefaultGlobalBasis<typename Impl::TaylorHoodBasisBuilder<GridView>::type>; {
using GlobalBasis
= Dune::Functions::DefaultGlobalBasis<typename Impl::TaylorHoodBasisBuilder<GridView>::type>;
};
} // end namespace AMDiS } // end namespace AMDiS
......
...@@ -83,12 +83,14 @@ namespace AMDiS ...@@ -83,12 +83,14 @@ namespace AMDiS
/// Return a reference to the data-matrix \ref matrix /// Return a reference to the data-matrix \ref matrix
BaseMatrix& getMatrix() BaseMatrix& getMatrix()
{ {
assert( !insertionMode );
return *matrix; return *matrix;
} }
/// Return a reference to the data-matrix \ref matrix /// Return a reference to the data-matrix \ref matrix
BaseMatrix const& getMatrix() const BaseMatrix const& getMatrix() const
{ {
assert( !insertionMode );
return *matrix; return *matrix;
} }
...@@ -115,12 +117,12 @@ namespace AMDiS ...@@ -115,12 +117,12 @@ namespace AMDiS
using BlockIndex = Dune::ReservedVector<size_type,1>; using BlockIndex = Dune::ReservedVector<size_type,1>;
/// \brief Returns an update-proxy of the inserter, to inserte/update a value at /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
/// position (\p r, \p c) in the matrix. Need an initialized inserter, that can /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
/// be created by \ref init and must be closed by \ref finish after insertion. /// be created by \ref init and must be closed by \ref finish after insertion.
auto operator()(FlatIndex row, FlatIndex col) auto operator()(FlatIndex row, FlatIndex col)
{ {
size_type r = row[0], c = col[0]; size_type r = row[0], c = col[0];
test_exit_dbg( initialized, "Inserter not initialized!"); test_exit_dbg( insertionMode, "Inserter not initilized!");
test_exit_dbg( r < N() && c < M() , test_exit_dbg( r < N() && c < M() ,
"Indices out of range [0,", N(), ")x[0,", M(), ")" ); "Indices out of range [0,", N(), ")x[0,", M(), ")" );
return (*inserter)[r][c]; return (*inserter)[r][c];
...@@ -129,7 +131,7 @@ namespace AMDiS ...@@ -129,7 +131,7 @@ namespace AMDiS
auto operator()(BlockIndex row, BlockIndex col) auto operator()(BlockIndex row, BlockIndex col)
{ {
size_type r = row[0], c = col[0]; size_type r = row[0], c = col[0];
test_exit_dbg( initialized, "Inserter not initialized!"); test_exit_dbg( insertionMode, "Inserter not initilized!");
test_exit_dbg( r < N() && c < M() , test_exit_dbg( r < N() && c < M() ,
"Indices out of range [0,", N(), ")x[0,", M(), ")" ); "Indices out of range [0,", N(), ")x[0,", M(), ")" );
return (*inserter)[r][c]; return (*inserter)[r][c];
...@@ -138,14 +140,14 @@ namespace AMDiS ...@@ -138,14 +140,14 @@ namespace AMDiS
/// Create inserter. Assumes that no inserter is currently active on this matrix. /// Create inserter. Assumes that no inserter is currently active on this matrix.