Skip to content

D
dune-amdis
Commit 339ba6ed authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files
Options

YaspGrid now works, Operators with factor argument

parent 5b2dce34
Pipeline #124 skipped
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1759 additions and 1078 deletions
dune/CMakeLists.txt
View file @ 339ba6ed
add_subdirectory(amdis)
add_subdirectory(amdis/test)
\ No newline at end of file
dune/amdis/AdaptBase.hpp
View file @ 339ba6ed
......@@ -42,7 +42,7 @@ namespace AMDiS
return name;
}
/// Returns \ref problemIteration_
/// Returns \ref problemIteration
ProblemIterationInterface* getProblemIteration() const
{
return problemIteration;
......@@ -60,7 +60,7 @@ namespace AMDiS
return adaptInfo;
}
/// Returns \ref problemTime_
/// Returns \ref problemTime
ProblemTimeInterface* getProblemTime() const
{
return problemTime;
......@@ -72,7 +72,7 @@ namespace AMDiS
problemTime = pti;
}
/// Returns \ref initialAdaptInfo_
/// Returns \ref initialAdaptInfo
AdaptInfo& getInitialAdaptInfo() const
{
return *initialAdaptInfo;
......@@ -93,7 +93,7 @@ namespace AMDiS
/** \brief
* Adapt info for initial adapt. Will be given to
* problemTime_->solveInitialProblem().
* problemTime->solveInitialProblem().
*/
AdaptInfo* initialAdaptInfo;
......
dune/amdis/Basic.hpp
View file @ 339ba6ed
......@@ -12,169 +12,169 @@
namespace AMDiS
{
template <int I>
using int_ = std::integral_constant<int, I>;
template <int I>
using int_ = std::integral_constant<int, I>;
template <bool B>
using bool_ = std::integral_constant<bool, B>;
template <bool B>
using bool_ = std::integral_constant<bool, B>;
template <size_t I>
using index_ = std::integral_constant<size_t, I>;
template <class T>
using IdxPairList = std::map< std::pair<int, int>, std::list<std::shared_ptr<T> > >;
template <class T>
using IdxList = std::map< int, std::list<std::shared_ptr<T> > >;
template <size_t I>
using index_ = std::integral_constant<size_t, I>;
template <size_t I, class T, class A>
T const& get(std::vector<T,A> const& vec)
{
return vec[I];
}
template <class T>
using IdxPairList = std::map< std::pair<int, int>, std::list<std::shared_ptr<T> > >;
namespace Impl
{
template <class Tuple, size_t N>
struct ConstructTuple
{
// add arg to repeated constructor argument list
template <class Arg, class... Args>
static Tuple make(Arg&& arg, Args&&... args)
{
return ConstructTuple<Tuple, N-1>::make(
std::forward<Arg>(arg), std::forward<Arg>(arg), std::forward<Args>(args)...);
}
};
template <class Tuple>
struct ConstructTuple<Tuple, 1>
{
template <class... Args>
static Tuple make(Args&&... args)
{
static_assert(std::tuple_size<Tuple>::value == sizeof...(args),
"Nr. of argument != tuple-size");
return Tuple{std::forward<Args>(args)...};
}
};
template <class Tuple>
struct ConstructTuple<Tuple, 0>
{
template <class... Args>
static Tuple make(Args&&... args)
{
static_assert(std::tuple_size<Tuple>::value == 0,
"Construction of empty tuples with empty argument list only!");
return {};
}
};
template <class Tuple>
struct FoldTuples
{
// add arg to repeated constructor argument list
template <size_t... Is, class... Tuples>
static Tuple make(Seq<Is...>, Tuples&&... tuples)
{
return Tuple{make_element(index_<Is>(), std::forward<Tuples>(tuples)...)...};
}
template <size_t I, class... Tuples>
static std::tuple_element_t<I, Tuple> make_element(index_<I>, Tuples&&... tuples)
{
using AMDiS::get;
return std::tuple_element_t<I, Tuple>{get<I>(std::forward<Tuples>(tuples))...};
}
};
} // end namespace Impl
// construct a tuple with each element constructed by the same argument arg.
template <class Tuple, class Arg>
Tuple construct_tuple(Arg&& arg)
{
return Impl::ConstructTuple<Tuple, std::tuple_size<Tuple>::value>::make(
std::forward<Arg>(arg));
}
template <class Tuple, class... Args>
Tuple fold_tuples(Args&&... args)
{
return Impl::FoldTuples<Tuple>::make(MakeSeq_t<std::tuple_size<Tuple>::value>(),
std::forward<Args>(args)...);
}
// -----------
template <template<class> class Base, class Tuple, class Indices> struct MakeTuple;
template <template<class> class Base, class Tuple, size_t... I>
struct MakeTuple<Base, Tuple, Seq<I...>>
template <class T>
using IdxList = std::map< int, std::list<std::shared_ptr<T> > >;
template <size_t I, class T, class A>
T const& get(std::vector<T,A> const& vec)
{
return vec[I];
}
namespace Impl
{
template <class Tuple, size_t N>
struct ConstructTuple
{
using type = std::tuple<Base<std::tuple_element_t<I, Tuple>>...>;
// add arg to repeated constructor argument list
template <class Arg, class... Args>
static Tuple make(Arg&& arg, Args&&... args)
{
return ConstructTuple<Tuple, N-1>::make(
std::forward<Arg>(arg), std::forward<Arg>(arg), std::forward<Args>(args)...);
}
};
template <template<class> class Base, class Tuple>
using MakeTuple_t =
typename MakeTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
// -----------
template <template<class,class> class Base, class Tuple1, class Tuple2, class Indices> struct MakeTuple2;
template <template<class,class> class Base, class Tuple1, class Tuple2, size_t... I>
struct MakeTuple2<Base, Tuple1, Tuple2, Seq<I...>>
template <class Tuple>
struct ConstructTuple<Tuple, 1>
{
using type = std::tuple<Base<std::tuple_element_t<I, Tuple1>, std::tuple_element_t<I, Tuple2>>...>;
template <class... Args>
static Tuple make(Args&&... args)
{
static_assert(std::tuple_size<Tuple>::value == sizeof...(args),
"Nr. of argument != tuple-size");
return Tuple{std::forward<Args>(args)...};
}
};
template <template<class,class> class Base, class Tuple1, class Tuple2>
using MakeTuple2_t =
typename MakeTuple2<Base, Tuple1, Tuple2, MakeSeq_t<std::tuple_size<Tuple1>::value>>::type;
// -----------
template <template<class...> class Base, class Tuple, class Indices> struct ExpandTuple;
template <template<class...> class Base, class Tuple, size_t... I>
struct ExpandTuple<Base, Tuple, Seq<I...>>
template <class Tuple>
struct ConstructTuple<Tuple, 0>
{
using type = Base<std::tuple_element_t<I, Tuple>...>;
template <class... Args>
static Tuple make(Args&&... args)
{
static_assert(std::tuple_size<Tuple>::value == 0,
"Construction of empty tuples with empty argument list only!");
return {};
}
};
template <template<class...> class Base, class Tuple>
using ExpandTuple_t =
typename ExpandTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
// -----------
template <class T, class Indices> struct RepeatedTuple;
template <class T, size_t... I>
struct RepeatedTuple<T, Seq<I...>>
template <class Tuple>
struct FoldTuples
{
template <size_t, class U> using Id = U;
using type = std::tuple<Id<I, T>...>;
// add arg to repeated constructor argument list
template <size_t... Is, class... Tuples>
static Tuple make(Seq<Is...>, Tuples&&... tuples)
{
return Tuple{make_element(index_<Is>(), std::forward<Tuples>(tuples)...)...};
}
template <size_t I, class... Tuples>
static std::tuple_element_t<I, Tuple> make_element(index_<I>, Tuples&&... tuples)
{
using AMDiS::get;
return std::tuple_element_t<I, Tuple>{get<I>(std::forward<Tuples>(tuples))...};
}
};
template <size_t N, class T>
using Repeat_t =
typename RepeatedTuple<T, MakeSeq_t<N>>::type;
} // end namespace Impl
// -----------
template <class T>
using owner = T;
// construct a tuple with each element constructed by the same argument arg.
template <class Tuple, class Arg>
Tuple construct_tuple(Arg&& arg)
{
return Impl::ConstructTuple<Tuple, std::tuple_size<Tuple>::value>::make(
std::forward<Arg>(arg));
}
template <class Tuple, class... Args>
Tuple fold_tuples(Args&&... args)
{
return Impl::FoldTuples<Tuple>::make(MakeSeq_t<std::tuple_size<Tuple>::value>(),
std::forward<Args>(args)...);
}
// -----------
template <template<class> class Base, class Tuple, class Indices> struct MakeTuple;
template <template<class> class Base, class Tuple, size_t... I>
struct MakeTuple<Base, Tuple, Seq<I...>>
{
using type = std::tuple<Base<std::tuple_element_t<I, Tuple>>...>;
};
template <template<class> class Base, class Tuple>
using MakeTuple_t =
typename MakeTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
// -----------
template <template<class,class> class Base, class Tuple1, class Tuple2, class Indices> struct MakeTuple2;
template <template<class,class> class Base, class Tuple1, class Tuple2, size_t... I>
struct MakeTuple2<Base, Tuple1, Tuple2, Seq<I...>>
{
using type = std::tuple<Base<std::tuple_element_t<I, Tuple1>, std::tuple_element_t<I, Tuple2>>...>;
};
template <template<class,class> class Base, class Tuple1, class Tuple2>
using MakeTuple2_t =
typename MakeTuple2<Base, Tuple1, Tuple2, MakeSeq_t<std::tuple_size<Tuple1>::value>>::type;
// -----------
template <template<class...> class Base, class Tuple, class Indices> struct ExpandTuple;
template <template<class...> class Base, class Tuple, size_t... I>
struct ExpandTuple<Base, Tuple, Seq<I...>>
{
using type = Base<std::tuple_element_t<I, Tuple>...>;
};
template <template<class...> class Base, class Tuple>
using ExpandTuple_t =
typename ExpandTuple<Base, Tuple, MakeSeq_t<std::tuple_size<Tuple>::value>>::type;
// -----------
template <class T, class Indices> struct RepeatedTuple;
template <class T, size_t... I>
struct RepeatedTuple<T, Seq<I...>>
{
template <size_t, class U> using Id = U;
using type = std::tuple<Id<I, T>...>;
};
template <size_t N, class T>
using Repeat_t =
typename RepeatedTuple<T, MakeSeq_t<N>>::type;
// -----------
template <class T>
using owner = T;
} // end namespace AMDiS
dune/amdis/ClonablePtr.hpp
View file @ 339ba6ed
......@@ -4,108 +4,107 @@
namespace AMDiS
{
// A pointer class that deletes only when owning the pointer
template <class T>
class ClonablePtr
{
private:
struct alloc_tag {}; ///< hidden helper struct, used by \ref make
public:
using Self = ClonablePtr;
using element_type = T;
// A pointer class that deletes only when owning the pointer
template <class T>
class ClonablePtr
{
private:
struct alloc_tag {}; ///< hidden helper struct, used by \ref make
public:
using Self = ClonablePtr;
using element_type = T;
/// Constructor from pointer. Can only be used via make method,
/// Transfers ownership.
ClonablePtr(owner<T>* p, alloc_tag) noexcept
: p(p)
, is_owner(true)
{}
/// Constructor from reference
ClonablePtr(T& ref) noexcept
: p(&ref)
, is_owner(false)
{}
/// Destructor, deletes in case of owner only
~ClonablePtr() noexcept
{
if (is_owner)
delete p;
}
/// Copy constructor, creates a clone of the pointed to object
ClonablePtr(Self const& that)
noexcept( std::is_nothrow_copy_constructible<T>::value )
: p(new T(*that.p))
, is_owner(true)
{}
/// Move constructor, copies the pointer only.
ClonablePtr(Self&& that) noexcept
: p(that.p)
, is_owner(that.is_owner)
{
that.p = NULL;
that.is_owner = false;
}
/// Copy and move assignment operator, using the copy-and-swap idiom
Self& operator=(Self that) noexcept
{
swap(that);
return *this;
}
/// Factory method. creates a new Object of type T and stores the pointer.
template <class... Args>
static Self make(Args&&... args)
noexcept( std::is_nothrow_constructible<T, std::decay_t<Args>...>::value )
{
return {new T(std::forward<Args>(args)...), Self::alloc_tag()};
}
/// Constructor from pointer. Can only be used via make method,
/// Transfers ownership.
ClonablePtr(owner<T>* p, alloc_tag) noexcept
: p(p)
, is_owner(true)
{}
/// Constructor from reference
ClonablePtr(T& ref) noexcept
: p(&ref)
, is_owner(false)
{}
/// Destructor, deletes in case of owner only
~ClonablePtr() noexcept
{
if (is_owner)
delete p;
}
/// Copy constructor, creates a clone of the pointed to object
ClonablePtr(Self const& that) noexcept( std::is_nothrow_copy_constructible<T>::value )
: p(new T(*that.p))
, is_owner(true)
{}
/// Move constructor, copies the pointer only.
ClonablePtr(Self&& that) noexcept
: p(that.p)
, is_owner(that.is_owner)
{
that.p = NULL;
that.is_owner = false;
}
/// Copy and move assignment operator, using the copy-and-swap idiom
Self& operator=(Self that) noexcept
{
swap(that);
return *this;
}
/// Factory method. creates a new Object of type T and stores the pointer.
template <class... Args>
static Self make(Args&&... args)
noexcept( std::is_nothrow_constructible<T, std::decay_t<Args>...>::value )
{
return {new T(std::forward<Args>(args)...), Self::alloc_tag()};
}
/// Access-method by dereferencing
T& operator*() const noexcept
{
return *p;
}
/// Access-method by pointer access
T* operator->() const noexcept
{
return p;
}
/// retrieve the underlying pointer
T* get() const noexcept
{
return p;
}
/// Test whether pointer is NULL
operator bool() const noexcept
{
return !(p == NULL);
}
void swap(Self& that) noexcept
{
using std::swap;
swap(p, that.p);
swap(is_owner, that.is_owner);
}
private:
T* p; ///< managed pointer
bool is_owner; ///< true, if class is owner of pointer, false otherwise
};
/// Access-method by dereferencing
T& operator*() const noexcept
{
return *p;
}
/// Access-method by pointer access
T* operator->() const noexcept
{
return p;
}
/// retrieve the underlying pointer
T* get() const noexcept
{
return p;
}
/// Test whether pointer is NULL
operator bool() const noexcept
{
return !(p == NULL);
}
template <class T>
void swap(ClonablePtr<T>& a, ClonablePtr<T>& b) noexcept
void swap(Self& that) noexcept
{
a.swap(b);
using std::swap;
swap(p, that.p);
swap(is_owner, that.is_owner);
}
private:
T* p; ///< managed pointer
bool is_owner; ///< true, if class is owner of pointer, false otherwise
};
template <class T>
void swap(ClonablePtr<T>& a, ClonablePtr<T>& b) noexcept
{
a.swap(b);
}
} // end namespace AMDiS
dune/amdis/DirichletBC.inc.hpp
View file @ 339ba6ed
......@@ -33,8 +33,8 @@ namespace AMDiS
matrix.clearDirichletRows(dirichletNodes, apply);
if (apply) {
interpolate(matrix.getRowFeSpace(), rhs.getVector(), values, dirichletNodes);
interpolate(matrix.getColFeSpace(), solution.getVector(), values, dirichletNodes);
interpolate(matrix.getRowFeSpace(), rhs.getVector(), values, dirichletNodes);
interpolate(matrix.getColFeSpace(), solution.getVector(), values, dirichletNodes);
}
}
......
dune/amdis/FileWriter.hpp
View file @ 339ba6ed
......@@ -15,9 +15,11 @@
namespace AMDiS
{
template <class Traits>
class FileWriter
{
template <class Traits>
class FileWriter
{
private: // typedefs and static constants
using Mesh = typename Traits::Mesh;
using MeshView = typename Mesh::LeafGridView;
......@@ -29,85 +31,108 @@ class FileWriter
/// Number of problem components
static constexpr int nComponents = Traits::nComponents;
public:
public:
/// Constructor.
FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names_)
: meshView(meshView)
, names(names_)
: meshView(meshView)
, names(names_)
{
std::string filename = "solution";
Parameters::get(base + "->filename", filename);
std::string dir = "output";
Parameters::get(base + "->output directory", dir);
vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
filename = "solution";
Parameters::get(base + "->filename", filename);
dir = "output";
Parameters::get(base + "->output directory", dir);
vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
}
/// default write method for time-depended data
template <class SystemVectorType>
void write(double time, SystemVectorType&& solutions)
{
vtkWriter->clear();
// copy dofvector to vertex data
For<0, nComponents>::loop([this, &solutions](const auto _i)
{
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
});
vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
vtkWriter->clear();
// copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i)
{
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
});
vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
}
/// default write method for stationary data
template <class SystemVectorType>
void write(SystemVectorType&& solutions)
{
vtkWriter->clear();
// copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i)
{
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
});
vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
}
protected:
protected:
template <class DOFVector, class Vector>
void dofVector2vertexVector(DOFVector dofvector, Vector& data)
void dofVector2vertexVector(DOFVector const& dofvector, Vector& data)
{
using Geometry = typename MeshView::template Codim<0>::Geometry;
using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
data.resize(meshView.size(dim));
auto const& indexSet = meshView.indexSet();
auto const& feSpace = dofvector.getFeSpace();
auto localView = feSpace.localView();
auto localIndexSet = feSpace.localIndexSet();
// copy data to P1-vector
for (auto const& element : elements(meshView)) {
localView.bind(element);
localIndexSet.bind(localView);
auto const& localBasis = localView.tree().finiteElement().localBasis();
auto const& refElement = RefElements::general(element.type());
std::vector<Dune::FieldVector<double,1> > shapeValues;
size_t nVertices = element.subEntities(dim);
for (size_t i = 0; i < nVertices; ++i) {
auto const& v = element.template subEntity<dim>(i);
auto pos = refElement.position(i, dim);
localBasis.evaluateFunction(pos, shapeValues);
size_t idx = indexSet.index(v);
data[idx] = 0.0;
for (size_t j = 0; j < shapeValues.size(); ++j) {
const auto global_idx = localIndexSet.index(j);
data[idx] += dofvector[global_idx] * shapeValues[j];
}
}
}
using Geometry = typename MeshView::template Codim<0>::Geometry;
using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
data.resize(meshView.size(dim));
auto const& indexSet = meshView.indexSet();
auto const& feSpace = dofvector.getFeSpace();