Skip to content

GitLab

Commit fd9ce2da authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

segfault in FixVec corrected

parent 85d9bec2
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ fd9ce2da
......@@ -17,6 +17,8 @@ include(DuneMacros)
dune_project()
dune_enable_all_packages(MODULE_LIBRARIES dunedec json11)
add_definitions(-Wall -Wextra -pedantic)
add_subdirectory("doc")
add_subdirectory("examples" EXCLUDE_FROM_ALL)
add_subdirectory("dune")
......
dune/dec/GridConcepts.hpp 0 → 100644
View file @ fd9ce2da
#pragma once
#include <dune/dec/common/ConceptsBase.hpp>
namespace Dec
{
namespace concepts
{
#ifndef DOXYGEN
namespace definition
{
struct Grid
{
template <class G, int dim = G::dimension, int dow = G::dimensionworld>
auto requires_(G&& grid) -> decltype(
concepts::valid_expr(
grid.leafGridView(),
grid.leafIndexSet(),
grid.levelGridView(int(0)),
grid.levelIndexSet(int(0)),
grid.globalRefine(int(0)),
grid.maxLevel()
));
};
struct GridView
{
template <class GV, int dim = GV::dimension>
auto requires_(GV&& gv) -> decltype(
concepts::valid_expr(
gv.grid(),
gv.indexSet(),
gv.size(int(0)),
gv.template begin<0>(),
gv.template end<0>(),
gv.level()
));
};
} // end namespace definition
#endif
/** \addtogroup concept
* @{
**/
/// \brief Argument `G` is a grid object.
/**
* Requirement:
* - ...
**/
template <class G>
constexpr bool Grid = models<definition::Grid(G)>;
/// \brief Argument `GV` is a grid-view object.
/**
* Requirement:
* - ...
**/
template <class GV>
constexpr bool GridView = models<definition::GridView(GV)>;
/** @} */
} // end namespace concepts
} // end namespace Dec
dune/dec/common/Index.hpp
View file @ fd9ce2da
......@@ -47,9 +47,6 @@ namespace Dec
/// Returns the size of the range
static constexpr auto size() { return std::integral_constant<Int, J-I>{}; }
/// Returns whether the range is not empty
constexpr operator bool() const { return !empty(); }
};
/// Extracts the Ith element element from the range
......
dune/dec/fixmatvec/FixMat.hpp
View file @ fd9ce2da
......@@ -98,7 +98,7 @@ namespace Dec
FixMat(MatExpr<Model, S, N, M> const& expr)
: Super()
{
*this = expr;
assign_impl(expr, [](auto& m, auto const& e) { m = e; });
}
/** @} **/
......@@ -176,6 +176,19 @@ namespace Dec
return static_cast<Super&>(*this)[idx];
}
/// Access the (i,j)th element of the matrix
T const& operator()(std::size_t i, std::size_t j) const
{
std::size_t const idx = i*M + j;
return static_cast<Super const&>(*this)[idx];
}
/// Access the i'th entry of the vector
T const& operator[](std::size_t idx) const
{
return static_cast<Super const&>(*this)[idx];
}
/// \brief Generates a sub-matrix of indices [I0,I1)x[J0,J1)
/**
* Example:
......@@ -223,8 +236,9 @@ namespace Dec
return {*this};
}
using Expr::operator(); // const element access
using Expr::operator[]; // const element access
using Expr::operator[];
// using Super::operator[];
using Expr::operator();
/// @cond HIDDEN_SYMBOLS
......@@ -242,7 +256,7 @@ namespace Dec
/// @endcond
/** @} **/
using Super::begin;
using Super::end;
using Super::cbegin;
......@@ -290,7 +304,7 @@ namespace Dec
FixMat<T,N,M> store_type_impl(MatExpr<E,T,N,M> const&);
template <class E>
struct StoreType<T, Void_t<decltype(store_type_impl<E>(std::declval<E>()))>>
struct StoreType<E, Void_t<decltype(store_type_impl<E>(std::declval<E>()))>>
{
using type = decltype(store_type_impl<E>(std::declval<E>()));
};
......
dune/dec/linear_algebra/block/BlockVector.hpp
View file @ fd9ce2da
......@@ -4,6 +4,7 @@
#include <tuple>
#include <vector>
#include <dune/dec/GridConcepts.hpp>
#include <dune/dec/LinearAlgebra.hpp>
#include <dune/dec/Mapper.hpp>
#include <dune/dec/linear_algebra/Expression.hpp>
......@@ -50,18 +51,21 @@ namespace Dec
BlockVector(std::vector<std::size_t> sizes)
: mapper_(sizes)
{
msg("constructor(sizes)");
Expects( N == sizes.size() );
for (std::size_t i = 0; i < N; ++i)
(*this)[i].resize(sizes[i]);
}
template <class Grid>
BlockVector(Grid const& grid, std::vector<int> Ks)
template <class GV,
REQUIRES( concepts::GridView<GV> )>
BlockVector(GV const& gv, std::vector<int> Ks)
{
msg("constructor(gv,sizes)");
Expects( N == Ks.size() );
std::vector<std::size_t> sizes(N);
for (std::size_t i = 0; i < N; ++i) {
sizes[i] = grid.size(Grid::dimension - Ks[i]);
sizes[i] = gv.size(GV::dimension - Ks[i]);
(*this)[i].resize(sizes[i]);
}
......@@ -73,6 +77,7 @@ namespace Dec
BlockVector(Splitter<Vector_> const& splitter)
: mapper_(splitter.mapper())
{
msg("constructor(splitter)");
auto const& input = splitter.vector();
for (std::size_t i = 0; i < N; ++i) {
(*this)[i].resize(mapper_.size(i));
......@@ -85,6 +90,7 @@ namespace Dec
: mapper_(that.mapper_)
, vector_(mapper_.size())
{
msg("constructor(that,ressource)");
for (std::size_t i = 0; i < N; ++i)
(*this)[i].resize(mapper_.size(i));
}
......@@ -92,6 +98,7 @@ namespace Dec
BlockVector& operator=(DOFVector<T> const& that)
{
msg("operator=(DOFVector)");
Expects( that.size() == size_type(mapper_.size()) );
for (std::size_t i = 0; i < N; ++i)
(*this)[i] = that.segment(mapper_.shift(i),mapper_.size(i));
......@@ -102,6 +109,7 @@ namespace Dec
template <class Vector_>
BlockVector& operator=(Splitter<Vector_> const& splitter)
{
msg("operator=(splitter)");
mapper_ = splitter.mapper();
auto const& input = splitter.vector();
for (std::size_t i = 0; i < N; ++i) {
......@@ -111,7 +119,7 @@ namespace Dec
return *this;
}
using Super::operator=;
// BlockVector& operator+=(BlockVector const& that)
......@@ -120,35 +128,35 @@ namespace Dec
// (*this)[i] += that[i];
// return *this;
// }
//
//
// BlockVector& operator-=(BlockVector const& that)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] -= that[i];
// return *this;
// }
//
//
// BlockVector& operator*=(float_type factor)
// {
// for (std::size_t i = 0; i < N; ++i)
// (*this)[i] *= factor;
// return *this;
// }
//
//
// template <class M>
// BlockVector& operator=(LinearAlgebraExpression<M> const& expr)
// {
// expr.assign(*this);
// return *this;
// }
//
//
// template <class M>
// BlockVector& operator+=(LinearAlgebraExpression<M> const& expr)
// {
// expr.add_assign(*this);
// return *this;
// }
//
//
// template <class M>
// BlockVector& operator-=(LinearAlgebraExpression<M> const& expr)
// {
......@@ -163,6 +171,13 @@ namespace Dec
(*this)[i].setZero();
}
/// Sets all vectors to zero
void setConstant(T value)
{
for (std::size_t i = 0; i < N; ++i)
(*this)[i].setConstant(value);
}
DOFVector<T>& vector() { return vector_; }
DOFVector<T> const& vector() const { return vector_; }
......
dune/dec/linear_algebra/eigen/DOFVector.hpp
View file @ fd9ce2da
......@@ -3,6 +3,7 @@
#include <Eigen/Dense>
#include <dune/dec/Dec.hpp>
#include <dune/dec/GridConcepts.hpp>
#include <dune/dec/common/ConceptsBase.hpp>
#include <dune/dec/common/Output.hpp>
#include <dune/dec/linear_algebra/DOFVectorBase.hpp>
......@@ -78,6 +79,7 @@ namespace Dec
: DOFVectorBase{size, name}
, Super(size)
{
msg("DOFVector::constructor()");
this->setZero();
}
......@@ -86,22 +88,27 @@ namespace Dec
: DOFVectorBase{size, name}
, Super(size)
{
msg("DOFVector::constructor(size,value)");
this->setConstant(value);
}
/// Constructor taking the grid and the dim 'K' of the entities, the value are located on.
template <class GridView,
REQUIRES(!std::is_arithmetic<GridView>::value)>
DOFVector(GridView const& gv, int K, std::string name = "")
: DOFVector{size_type (gv.size(GridView::dimension - K)), name}
{}
template <class GV,
REQUIRES( concepts::GridView<GV>)>
DOFVector(GV const& gv, int K, std::string name = "")
: DOFVector{size_type (gv.size(GV::dimension - K)), name}
{
msg("DOFVector::constructor(gv)");
}
/// Constructor taking the grid and the dim 'K' of the entities, the value are located on.
template <class GridView,
REQUIRES(!std::is_arithmetic<GridView>::value)>
DOFVector(GridView const& gv, int K, value_type value, std::string name = "")
: DOFVector{size_type (gv.size(GridView::dimension - K)), value, name}
{}
template <class GV,
REQUIRES( concepts::GridView<GV>)>
DOFVector(GV const& gv, int K, value_type value, std::string name = "")
: DOFVector{size_type (gv.size(GV::dimension - K)), value, name}
{
msg("DOFVector::constructor(gv,value)");
}
/// Copy/Move constructor
DOFVector(DOFVector const& that) = default;
......@@ -111,7 +118,9 @@ namespace Dec
DOFVector(DOFVector const& that, tag::ressource)
: DOFVectorBase{that.size_, that.name_}
, Super(that.size_)
{}
{
msg("DOFVector::constructor(that,ressource)");
}
/// Destructor.
~DOFVector() = default;
......@@ -123,12 +132,16 @@ namespace Dec
DOFVector(Eigen::MatrixBase<ThatDerived> const& that)
: DOFVectorBase{size_type(that.size()), "expressions"}
, Super(that)
{}
{
msg("DOFVector::constructor(matrix-expr)");
}
// This method allows you to assign Eigen expressions to DOFVector
template <class ThatDerived>
DOFVector& operator=(Eigen::MatrixBase<ThatDerived> const& that)
{
msg("DOFVector::operator=(matrix-expr)");
Super::operator=(that);
return *this;
}
......@@ -138,12 +151,16 @@ namespace Dec
DOFVector(Eigen::ArrayBase<ThatDerived> const& that)
: DOFVectorBase{size_type(that.size()), "expressions"}
, Super(that)
{}
{
msg("DOFVector::constructor(array-expr)");
}
// This method allows you to assign Eigen expressions to DOFVector
template <class ThatDerived>
DOFVector& operator=(Eigen::ArrayBase<ThatDerived> const& that)
{
msg("DOFVector::operator=(array-expr)");
Super::operator=(that);
return *this;
}
......
dune/dec/linear_algebra/interface/Solver.hpp
View file @ fd9ce2da
......@@ -16,12 +16,12 @@ namespace Dec
* \note Should be implemented by derived class.
**/
template <class Matrix>
Impl& compute(Matrix const& A)
{
Impl& compute(Matrix const& /*A*/)
{
error_exit("Should be implemented by derived class.");
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// condition is zero.
/**
* Requirement:
......@@ -33,8 +33,8 @@ namespace Dec
u.setZero();
solveWithGuessImpl(b, u);
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// condition is zero using iteration.
/**
* Requirement:
......@@ -51,15 +51,15 @@ namespace Dec
/**
* Requirement:
* - \ref compute() must be called in advance
*
*
* \note Should be implemented by derived class.
**/
template <class VectorB, class VectorU>
void solveWithGuess(VectorB const& b, VectorU& u) const
void solveWithGuess(VectorB const& /*b*/, VectorU& /*u*/) const
{
error_exit("Should be implemented by derived class.");
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ using iteration.
/**
* Requirement:
......@@ -71,14 +71,14 @@ namespace Dec
for (iter.init(b); !iter.finished(); ++iter)
solveWithGuessImpl(b, u);
}
private:
// redirect to implementation
template <class VectorB, class VectorU>
void solveWithGuessImpl(VectorB const& b, VectorU& u) const
{
static_cast<Impl const&>(*this).solveWithGuess(b, u);
static_cast<Impl const&>(*this).solveWithGuess(b, u);
}
};
......
dune/dec/linear_algebra/smoother/SOR.hpp
View file @ fd9ce2da
......@@ -58,7 +58,7 @@ namespace Dec
}
template <bool B>
void apply(DOFMatrix<float_type> const& A, DOFVector<float_type>& u, DOFVector<float_type> const& b, bool_t<B>) const
void apply(DOFMatrix<float_type> const& /*A*/, DOFVector<float_type>& /*u*/, DOFVector<float_type> const& /*b*/, bool_t<B>) const
{
static_assert(B, "Not implemented!");
}
......
dune/dec/linear_algebra/smoother/SSOR.hpp
View file @ fd9ce2da
......@@ -73,7 +73,7 @@ namespace Dec
}
template <bool B>
void apply(DOFMatrix<float_type> const& A, DOFVector<float_type>& u, DOFVector<float_type> const& b, bool_t<B>) const
void apply(DOFMatrix<float_type> const& /*A*/, DOFVector<float_type>& /*u*/, DOFVector<float_type> const& /*b*/, bool_t<B>) const
{
static_assert(B, "Not implemented!");
}
......
examples/frank_oseen.cpp
View file @ fd9ce2da
......@@ -160,8 +160,8 @@ int main(int argc, char** argv)
}
}
BVector b(grid, {1,1}); // rhs vector
BVector solution(grid, {1,1});
BVector b(leafGridView, {1,1}); // rhs vector
BVector solution(leafGridView, {1,1});
DEC_MSG("size(0) = ",grid.size(0),", size(1) = ",grid.size(1),", size(2) = ",grid.size(2));
......
test/CMakeLists.txt
View file @ fd9ce2da
set(CMAKE_MEMORYCHECK_COMMAND valgrind)
function(add_memcheck_test name binary)
set(memcheck_command "${CMAKE_MEMORYCHECK_COMMAND} ${CMAKE_MEMORYCHECK_COMMAND_OPTIONS}")
separate_arguments(memcheck_command)
add_test(${name} ${binary} ${ARGN})
add_test(memcheck_${name} ${memcheck_command} ./${binary} ${ARGN})
endfunction(add_memcheck_test)
function(set_memcheck_test_properties name)
set_tests_properties(${name} ${ARGN})
set_tests_properties(memcheck_${name} ${ARGN})
endfunction(set_memcheck_test_properties)
add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND})
set(DEC_DIM 2)
set(DEC_DOW 2)
# add a test program
add_executable(test_filesystem test_filesystem.cc)
target_link_libraries(test_filesystem dunedec)
add_test(NAME RunTestFilesystem COMMAND test_filesystem)
add_memcheck_test(RunTestFilesystem test_filesystem)
if (DEC_HAS_EIGEN)
add_executable(test_bernstein test_bernstein.cc)
......@@ -17,7 +31,7 @@ endif ()
add_executable(test_simplegrid test_simplegrid.cc)
target_link_libraries(test_simplegrid dunedec)
add_test(NAME RunTestSimpleGrid COMMAND test_simplegrid)
add_memcheck_test(RunTestSimpleGrid test_simplegrid)
# add_executable(test_simplicialcomplex test_simplicialcomplex.cc)
# target_link_libraries(test_simplicialcomplex dunedec)
......@@ -25,43 +39,43 @@ add_test(NAME RunTestSimpleGrid COMMAND test_simplegrid)
add_executable(test_linearalgebra test_linearalgebra.cc)
target_link_libraries(test_linearalgebra dunedec)
add_test(NAME RunTestLinearAlgebra COMMAND test_linearalgebra)
add_memcheck_test(RunTestLinearAlgebra test_linearalgebra)
# add_dec_executable(test_gridcheck test_gridcheck.cc)
# add_test(NAME RunTestGridCheck COMMAND test_gridcheck)
add_dec_executable(test_gridfactory ALBERTA test_gridfactory.cc)
add_test(NAME RunTestGridFactory COMMAND test_gridfactory)
add_memcheck_test(RunTestGridFactory test_gridfactory)
add_dec_executable(test_half_edges test_half_edges.cc)
add_test(NAME RunTestHalfEdges COMMAND test_half_edges)
add_memcheck_test(RunTestHalfEdges test_half_edges)
add_dec_executable(test_laplace_RR test_laplace_RR.cc)
add_test(NAME RunTestLaplaceRR COMMAND test_laplace_RR)
add_memcheck_test(RunTestLaplaceRR test_laplace_RR)
add_dec_executable(test_ranges test_ranges.cc)
add_test(NAME RunTestRanges COMMAND test_ranges)
add_memcheck_test(RunTestRanges test_ranges)
add_dec_executable(test_solver test_solver.cc)
add_test(NAME RunTestSolver COMMAND test_solver)
add_memcheck_test(RunTestSolver test_solver)
add_dec_executable(test_array test_array.cc)
add_test(NAME RunTestArray COMMAND test_array)
add_memcheck_test(RunTestArray test_array)
add_dec_executable(test_capacity_vector test_capacity_vector.cc)
add_test(NAME RunTestCapacityVector COMMAND test_capacity_vector)
add_memcheck_test(RunTestCapacityVector test_capacity_vector)
add_dec_executable(test_cast test_cast.cc)
add_test(NAME RunTestCast COMMAND test_cast)
add_memcheck_test(RunTestCast test_cast)
add_dec_executable(test_fixvec test_fixvec.cc)
add_test(NAME RunTestFixVec COMMAND test_fixvec)
add_memcheck_test(RunTestFixVec test_fixvec)
add_dec_executable(test_string test_string.cc)
add_test(NAME RunTestString COMMAND test_string)
add_memcheck_test(RunTestString test_string)
add_dec_executable(test_block test_block.cc)
add_test(NAME RunTestBlock COMMAND test_block)
add_memcheck_test(RunTestBlock test_block)
add_dependencies(check
test_filesystem
......
test/test_block.cc
View file @ fd9ce2da
......@@ -38,29 +38,34 @@ namespace test {
// build up laplace matrix
LaplaceBeltrami<GridView> laplace(gv);
BlockMatrix<double, 2, 2> A;
laplace.build(A(0,0));
laplace.build(A(1,1));
A(0,0).clear_dirichlet_row(0);
A(1,1).clear_dirichlet_row(0);
// build rhs vector
BlockVector<double, 2> b(grid, {0,0});
b[0] = 1.0;
b[1] = 1.0;
msg("build rhs vector...");
BlockVector<double, 2> b(grid.leafGridView(), {0,0});
msg("b.setConstant");
b.setConstant(1.0);
// b[1].setConstant(1.0);
msg("construct x from b...");
BlockVector<double, 2> x(b, tag::ressource{});
msg("x.setZero");
x.setZero();
msg("bcgs...");
solver::bcgs(A, x, b);
auto v1 = unary_dot(b);
auto v2 = two_norm(b);
auto v3 = dot(b, x);
}
}}
......
test/test_linearalgebra.cc
View file @ fd9ce2da
......@@ -41,6 +41,8 @@ namespace test {
DEC_TEST_EQ( v2[0], v3[0] );
v1.setConstant(1.0);
#ifndef DEC_HAS_PETSC
// extract a sub-vector
auto sub_v2 = v2.segment(0, 5);
......
test/test_solver.cc
View file @ fd9ce2da
......@@ -45,43 +45,43 @@ namespace test {
A.clear_dirichlet_row(0);