From 339ba6ed96c8d9f46fae8decf6795f13e5878b9b Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 27 Apr 2016 18:39:22 +0200
Subject: [PATCH] YaspGrid now works, Operators with factor argument
---
dune/CMakeLists.txt | 1 +
dune/amdis/AdaptBase.hpp | 8 +-
dune/amdis/Basic.hpp | 292 +++++++-------
dune/amdis/ClonablePtr.hpp | 195 +++++----
dune/amdis/DirichletBC.inc.hpp | 4 +-
dune/amdis/FileWriter.hpp | 149 ++++---
dune/amdis/IndexSeq.hpp | 11 +-
dune/amdis/Literals.hpp | 2 +
dune/amdis/Log.hpp | 4 +-
dune/amdis/Loops.hpp | 9 +-
dune/amdis/Math.hpp | 87 ++--
dune/amdis/Operator.hpp | 94 +++--
dune/amdis/Operator.inc.hpp | 347 +++++++++-------
dune/amdis/OperatorTerm.hpp | 381 +++++++++---------
dune/amdis/OperatorTermBase.hpp | 126 +++---
dune/amdis/ProblemInstat.hpp | 2 +-
dune/amdis/ProblemStat.hpp | 160 +++++---
dune/amdis/ProblemStat.inc.hpp | 136 ++++---
dune/amdis/Traits.hpp | 92 ++---
.../linear_algebra/LinearSolverInterface.hpp | 3 +-
.../linear_algebra/mtl/BlockMTLMatrix.hpp | 7 +-
.../linear_algebra/mtl/BlockMTLVector.hpp | 11 +-
dune/amdis/linear_algebra/mtl/DOFVector.hpp | 19 +-
.../amdis/linear_algebra/mtl/KrylovRunner.hpp | 2 -
.../amdis/linear_algebra/mtl/SystemVector.hpp | 31 +-
dune/amdis/test/CMakeLists.txt | 18 +
dune/amdis/test/kdtree.hpp | 245 +++++++++++
dune/amdis/test/macro.stand.2d | 30 ++
dune/amdis/test/test.json.2d | 29 ++
dune/amdis/test/test1.cc | 150 +++++++
init/heat.json.2d | 4 +-
init/test.json.2d | 29 ++
src/CMakeLists.txt | 2 +-
src/heat.cc | 112 ++---
src/navier_stokes.cc | 6 +-
src/pfc.cc | 35 +-
src/stokes.cc | 4 +-
37 files changed, 1759 insertions(+), 1078 deletions(-)
create mode 100644 dune/amdis/test/CMakeLists.txt
create mode 100644 dune/amdis/test/kdtree.hpp
create mode 100644 dune/amdis/test/macro.stand.2d
create mode 100644 dune/amdis/test/test.json.2d
create mode 100644 dune/amdis/test/test1.cc
create mode 100644 init/test.json.2d
diff --git a/dune/CMakeLists.txt b/dune/CMakeLists.txt
index 40be2be0..81c0d133 100644
--- a/dune/CMakeLists.txt
+++ b/dune/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory(amdis)
+add_subdirectory(amdis/test)
\ No newline at end of file
diff --git a/dune/amdis/AdaptBase.hpp b/dune/amdis/AdaptBase.hpp
index 31f642e5..5d0422f1 100644
--- a/dune/amdis/AdaptBase.hpp
+++ b/dune/amdis/AdaptBase.hpp
@@ -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;
diff --git a/dune/amdis/Basic.hpp b/dune/amdis/Basic.hpp
index 9f4f8990..121cf7bd 100644
--- a/dune/amdis/Basic.hpp
+++ b/dune/amdis/Basic.hpp
@@ -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
diff --git a/dune/amdis/ClonablePtr.hpp b/dune/amdis/ClonablePtr.hpp
index 4f2b384e..aa0076b7 100644
--- a/dune/amdis/ClonablePtr.hpp
+++ b/dune/amdis/ClonablePtr.hpp
@@ -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
diff --git a/dune/amdis/DirichletBC.inc.hpp b/dune/amdis/DirichletBC.inc.hpp
index 57174396..e4620ad1 100644
--- a/dune/amdis/DirichletBC.inc.hpp
+++ b/dune/amdis/DirichletBC.inc.hpp
@@ -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);
}
}
diff --git a/dune/amdis/FileWriter.hpp b/dune/amdis/FileWriter.hpp
index 717daef4..a411b253 100644
--- a/dune/amdis/FileWriter.hpp
+++ b/dune/amdis/FileWriter.hpp
@@ -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();
+
+ 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];
+ }
+ }
+ }
}
-private:
+ private:
MeshView const& meshView;
std::shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
std::shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
std::vector<std::string> names;
std::array<std::vector<double>, nComponents> data_vectors;
-};
+
+ std::string filename;
+ std::string dir;
+ };
} // end namespace AMDiS
diff --git a/dune/amdis/IndexSeq.hpp b/dune/amdis/IndexSeq.hpp
index 942f34fa..4f4c504d 100644
--- a/dune/amdis/IndexSeq.hpp
+++ b/dune/amdis/IndexSeq.hpp
@@ -3,9 +3,10 @@
// from http://stackoverflow.com/questions/17424477/implementation-c14-make-integer-sequence/17426611#17426611
namespace AMDiS
{
+ /// class that represents a sequence of integers
template <size_t...> struct Seq {};
- namespace detail
+ namespace Impl
{
template <size_t s, class S>
struct Concat;
@@ -28,19 +29,19 @@ namespace AMDiS
using type = Seq<Is..., sizeof...(Is)>;
};
- } // end namespace detail
+ } // end namespace Impl
template <size_t N>
struct MakeSeq
{
- using type = typename detail::IncSeq_if< (N % 2 != 0),
- typename detail::Concat<N/2, typename MakeSeq<N/2>::type>::type >::type;
+ using type = typename Impl::IncSeq_if< (N % 2 != 0),
+ typename Impl::Concat<N/2, typename MakeSeq<N/2>::type>::type >::type;
};
// break condition
template <> struct MakeSeq<0> { using type = Seq<>; };
- // alias template
+ /// Alias template to create a sequence of integers
template <size_t N>
using MakeSeq_t = typename MakeSeq<N>::type;
diff --git a/dune/amdis/Literals.hpp b/dune/amdis/Literals.hpp
index 4a011d65..32a81f5a 100644
--- a/dune/amdis/Literals.hpp
+++ b/dune/amdis/Literals.hpp
@@ -5,6 +5,7 @@
namespace AMDiS
{
// inspired by Boost.hana
+ // see also: http://blog.mattbierner.com/stupid-template-tricks-stdintegral_constant-user-defined-literal/
namespace Impl
{
@@ -33,6 +34,7 @@ namespace AMDiS
} // end namespace Impl
+ /// Literal to create integer compile-time constant, e.g. 0_c -> index_<0>
template <char... digits>
constexpr auto operator"" _c()
{
diff --git a/dune/amdis/Log.hpp b/dune/amdis/Log.hpp
index d183807c..cdabe9d7 100644
--- a/dune/amdis/Log.hpp
+++ b/dune/amdis/Log.hpp
@@ -30,9 +30,11 @@
// ===== message macros ======================================================
// ===========================================================================
+#define AMDIS_UNUSED(var) __attribute__((unused)) var
+
/// Should be the first call in every functions. It defines the current
/// function name nn for message output via MSG, WARNING, ...
-#define AMDIS_FUNCNAME(nn) const char *funcName; funcName = nn;
+#define AMDIS_FUNCNAME(nn) AMDIS_UNUSED(const char *funcName); funcName = nn;
#ifdef NDEBUG
#define AMDIS_FUNCNAME_DBG(nn)
diff --git a/dune/amdis/Loops.hpp b/dune/amdis/Loops.hpp
index 776dbad1..fc964cdd 100644
--- a/dune/amdis/Loops.hpp
+++ b/dune/amdis/Loops.hpp
@@ -4,7 +4,7 @@
namespace AMDiS
{
- // class that performes the loop over all indices [I, N), for i < N
+ /// class that performes the loop over all indices [I, N), for i < N
template <size_t I, size_t N>
struct For
{
@@ -20,4 +20,11 @@ namespace AMDiS
template <size_t N>
struct For<N, N> { template <class F> static void loop(F) {} };
+ /// Function to call static loop \ref For
+ template <size_t I, size_t N, class F>
+ std::enable_if_t< (I <= N) > for_(F&& f)
+ {
+ For<I, N>::loop(std::forward<F>(f));
+ }
+
} // end namespace AMDiS
diff --git a/dune/amdis/Math.hpp b/dune/amdis/Math.hpp
index 6d6d98ff..bc6e353f 100644
--- a/dune/amdis/Math.hpp
+++ b/dune/amdis/Math.hpp
@@ -10,64 +10,61 @@
namespace AMDiS
{
- namespace math
+ namespace math
+ {
+ template <class T>
+ std::enable_if_t<std::is_arithmetic<T>::value, T>
+ constexpr abs(T a)
{
- template <class T>
- std::enable_if_t<std::is_arithmetic<T>::value, T>
- constexpr abs(T a)
- {
- return a >= 0 ? a : -a;
- }
+ return a >= 0 ? a : -a;
+ }
- template <class T>
- std::enable_if_t<std::is_arithmetic<T>::value, T>
- constexpr sqr(T a)
- {
- return a*a;
- }
+ template <class T>
+ std::enable_if_t<std::is_arithmetic<T>::value, T>
+ constexpr sqr(T a)
+ {
+ return a*a;
+ }
- template <class T0, class T1>
- constexpr auto min(T0 a, T1 b)
- {
- return a > b ? b : a;
- }
+ template <class T0, class T1>
+ constexpr auto min(T0 a, T1 b)
+ {
+ return a > b ? b : a;
+ }
- template <class T0, class T1>
- constexpr auto max(T0 a, T1 b)
- {
- return a > b ? a : b;
- }
+ template <class T0, class T1>
+ constexpr auto max(T0 a, T1 b)
+ {
+ return a > b ? a : b;
+ }
- } // end namespace math
+ } // end namespace math
- template <class T> inline void nullify(T& a)
- {
- a = 0;
- }
+ template <class T> inline void nullify(T& a)
+ {
+ a = 0;
+ }
- inline void nullify(std::string& s)
- {
- s = "";
- }
+ inline void nullify(std::string& s)
+ {
+ s = "";
+ }
- /// Calculates factorial of i
- inline long fac(long i)
- {
- if (i <= 1)
- return 1;
- else
- return i * fac(i - 1);
- }
+ /// Calculates factorial of i
+ constexpr long factorial(long i)
+ {
+ return i <= 1 ? 1 : i * factorial(i - 1);
+ }
- /// check for inf and nan values
- inline bool isNumber(double val)
- {
- return !std::isnan(val) && !std::isinf(val);
- }
+ /// check for inf and nan values
+ inline bool isNumber(double val)
+ {
+ return !std::isnan(val) && !std::isinf(val);
+ }
// ===== some predefined values ===============================================
constexpr double m_e = 2.7182818284590452354;
diff --git a/dune/amdis/Operator.hpp b/dune/amdis/Operator.hpp
index 52fd43f7..e307c0fc 100644
--- a/dune/amdis/Operator.hpp
+++ b/dune/amdis/Operator.hpp
@@ -8,8 +8,8 @@ namespace AMDiS
{
enum FirstOrderType
{
- GRD_PSI,
- GRD_PHI
+ GRD_PSI,
+ GRD_PHI
};
@@ -25,57 +25,62 @@ namespace AMDiS
ColFeSpace const& colFeSpace);
template <class RowView, class ColView, class ElementMatrix>
- void getElementMatrix(RowView const& rowView,
+ bool getElementMatrix(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix);
+ ElementMatrix& elementMatrix,
+ double* factor = NULL);
template <class RowView, class ElementVector>
- void getElementVector(RowView const& rowView,
- ElementVector& elementVector);
+ bool getElementVector(RowView const& rowView,
+ ElementVector& elementVector,
+ double* factor = NULL);
+
template <class Term>
Self& addZOT(Term const& term)
{
- zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
- return *this;
+ zeroOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
+ return *this;
}
+
template <class Term>
Self& addFOT(Term const& term, FirstOrderType firstOrderType)
{
- using OpTerm = GenericOperatorTerm<MeshView, Term>;
- if (firstOrderType == GRD_PHI)
- firstOrderGrdPhi.push_back(new OpTerm(term));
- else
- firstOrderGrdPsi.push_back(new OpTerm(term));
- return *this;
+ using OpTerm = GenericOperatorTerm<MeshView, Term>;
+ if (firstOrderType == GRD_PHI)
+ firstOrderGrdPhi.push_back(new OpTerm(term));
+ else
+ firstOrderGrdPsi.push_back(new OpTerm(term));
+ return *this;
}
template <class Term, size_t I>
Self& addFOT(Term const& term, const index_<I>, FirstOrderType firstOrderType)
{
- using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent<I>>;
-
- if (firstOrderType == GRD_PHI)
- firstOrderGrdPhi.push_back(new OpTerm(term));
- else
- firstOrderGrdPsi.push_back(new OpTerm(term));
- return *this;
+ using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent<I>>;
+
+ if (firstOrderType == GRD_PHI)
+ firstOrderGrdPhi.push_back(new OpTerm(term));
+ else
+ firstOrderGrdPsi.push_back(new OpTerm(term));
+ return *this;
}
+
template <class Term>
Self& addSOT(Term const& term)
{
- secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
- return *this;
+ secondOrder.push_back(new GenericOperatorTerm<MeshView, Term>(term));
+ return *this;
}
template <class Term, size_t I, size_t J>
Self& addSOT(Term const& term, const index_<I>, const index_<J>)
{
- using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent<I,J>>;
- secondOrder.push_back(new OpTerm(term));
- return *this;
+ using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent<I,J>>;
+ secondOrder.push_back(new OpTerm(term));
+ return *this;
}
@@ -87,42 +92,47 @@ namespace AMDiS
template <class RowView, class ColView, class ElementMatrix>
void assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix);
+ ElementMatrix& elementMatrix,
+ double factor);
template <class RowView, class ElementVector>
void assembleZeroOrderTerms(RowView const& rowView,
- ElementVector& elementVector);
+ ElementVector& elementVector,
+ double factor);
template <class RowView, class ColView, class ElementMatrix>
void assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix);
+ ElementMatrix& elementMatrix,
+ double factor);
template <class RowView, class ColView, class ElementMatrix>
void assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix);
+ ElementMatrix& elementMatrix,
+ double factor);
template <class RowView, class ColView, class ElementMatrix>
void assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix);
+ ElementMatrix& elementMatrix,
+ double factor);
private:
- /// List of all second order terms
- std::list<OperatorTermType*> secondOrder;
+ /// List of all second order terms
+ std::list<OperatorTermType*> secondOrder;
- /// List of all first order terms derived to psi
- std::list<OperatorTermType*> firstOrderGrdPsi;
+ /// List of all first order terms derived to psi
+ std::list<OperatorTermType*> firstOrderGrdPsi;
- /// List of all first order terms derived to phi
- std::list<OperatorTermType*> firstOrderGrdPhi;
+ /// List of all first order terms derived to phi
+ std::list<OperatorTermType*> firstOrderGrdPhi;
- /// List of all zero order terms
- std::list<OperatorTermType*> zeroOrder;
-
- int psiDegree = 1;
- int phiDegree = 1;
+ /// List of all zero order terms
+ std::list<OperatorTermType*> zeroOrder;
+
+ int psiDegree = 1;
+ int phiDegree = 1;
};
} // end namespace AMDiS
diff --git a/dune/amdis/Operator.inc.hpp b/dune/amdis/Operator.inc.hpp
index 118668bb..c5a8d9e8 100644
--- a/dune/amdis/Operator.inc.hpp
+++ b/dune/amdis/Operator.inc.hpp
@@ -9,6 +9,8 @@ namespace AMDiS
void Operator<MeshView>::init(RowFeSpace const& rowFeSpace,
ColFeSpace const& colFeSpace)
{
+ AMDIS_FUNCNAME("Operator::init()");
+
// const auto& rowFE = rowView.tree().finiteElement();
// const auto& colFE = colView.tree().finiteElement();
@@ -24,28 +26,48 @@ namespace AMDiS
template <class MeshView>
template <class RowView, class ColView, class ElementMatrix>
- void Operator<MeshView>::getElementMatrix(RowView const& rowView,
+ bool Operator<MeshView>::getElementMatrix(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix)
+ ElementMatrix& elementMatrix,
+ double* factor)
{
+ AMDIS_FUNCNAME("Operator::getElementMatrix()");
+
+ double fac = factor ? *factor : 1.0;
+
+ if (fac == 0.0)
+ return false;
+
if (!zeroOrder.empty())
- assembleZeroOrderTerms(rowView, colView, elementMatrix);
+ assembleZeroOrderTerms(rowView, colView, elementMatrix, fac);
if (!firstOrderGrdPhi.empty())
- assembleFirstOrderTermsGrdPhi(rowView, colView, elementMatrix);
+ assembleFirstOrderTermsGrdPhi(rowView, colView, elementMatrix, fac);
if (!firstOrderGrdPsi.empty())
- assembleFirstOrderTermsGrdPsi(rowView, colView, elementMatrix);
+ assembleFirstOrderTermsGrdPsi(rowView, colView, elementMatrix, fac);
if (!secondOrder.empty())
- assembleSecondOrderTerms(rowView, colView, elementMatrix);
+ assembleSecondOrderTerms(rowView, colView, elementMatrix, fac);
+
+ return true;
}
template <class MeshView>
template <class RowView, class ElementVector>
- void Operator<MeshView>::getElementVector(RowView const& rowView,
- ElementVector& elementVector)
- {
+ bool Operator<MeshView>::getElementVector(RowView const& rowView,
+ ElementVector& elementVector,
+ double* factor)
+ {
+ AMDIS_FUNCNAME("Operator::getElementVector()");
+
+ double fac = factor ? *factor : 1.0;
+
+ if (fac == 0.0)
+ return false;
+
if (!zeroOrder.empty())
- assembleZeroOrderTerms(rowView, elementVector);
+ assembleZeroOrderTerms(rowView, elementVector, fac);
+
+ return true;
}
@@ -53,8 +75,11 @@ namespace AMDiS
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix)
+ ElementMatrix& elementMatrix,
+ double factor)
{
+ AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementMatrix)");
+
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
@@ -67,33 +92,33 @@ namespace AMDiS
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : zeroOrder)
- operatorTerm->init(element, quad);
+ operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
- // Position of the current quadrature point in the reference element
- const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
-
- // The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos);
-
- std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
- rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
-
- if (psiDegree == phiDegree)
- colShapeValues = rowShapeValues;
- else
- colLocalBasis.evaluateFunction(quadPos, colShapeValues);
-
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
- const int local_i = rowView.tree().localIndex(i);
- const int local_j = colView.tree().localIndex(j);
- for (auto* operatorTerm : zeroOrder)
- elementMatrix[local_i][local_j]
- += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j])
- * quad[iq].weight() * integrationElement;
- }
- }
+ // Position of the current quadrature point in the reference element
+ const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
+
+ // The multiplicative factor in the integral transformation formula
+ const double integrationElement = geometry.integrationElement(quadPos) * factor;
+
+ std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
+ rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
+
+ if (psiDegree == phiDegree)
+ colShapeValues = rowShapeValues;
+ else
+ colLocalBasis.evaluateFunction(quadPos, colShapeValues);
+
+ for (size_t i = 0; i < elementMatrix.N(); ++i) {
+ for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ const int local_i = rowView.tree().localIndex(i);
+ const int local_j = colView.tree().localIndex(j);
+ for (auto* operatorTerm : zeroOrder)
+ elementMatrix[local_i][local_j]
+ += operatorTerm->evalZot(iq, rowShapeValues[i], colShapeValues[j])
+ * quad[iq].weight() * integrationElement;
+ }
+ }
}
}
@@ -102,8 +127,11 @@ namespace AMDiS
template <class MeshView>
template <class RowView, class ElementVector>
void Operator<MeshView>::assembleZeroOrderTerms(RowView const& rowView,
- ElementVector& elementvector)
+ ElementVector& elementvector,
+ double factor)
{
+ AMDIS_FUNCNAME("Operator::assembleZeroOrderTerms(elementvector)");
+
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
@@ -118,22 +146,22 @@ namespace AMDiS
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
- // Position of the current quadrature point in the reference element
- const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
-
- // The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos);
-
- std::vector<Dune::FieldVector<double,1> > rowShapeValues;
- rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
-
- for (size_t i = 0; i < elementvector.size(); ++i) {
- const int local_i = rowView.tree().localIndex(i);
- for (auto* operatorTerm : zeroOrder)
- elementvector[local_i]
- += operatorTerm->evalZot(iq, rowShapeValues[i])
- * quad[iq].weight() * integrationElement;
- }
+ // Position of the current quadrature point in the reference element
+ const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
+
+ // The multiplicative factor in the integral transformation formula
+ const double integrationElement = geometry.integrationElement(quadPos) * factor;
+
+ std::vector<Dune::FieldVector<double,1> > rowShapeValues;
+ rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
+
+ for (size_t i = 0; i < elementvector.size(); ++i) {
+ const int local_i = rowView.tree().localIndex(i);
+ for (auto* operatorTerm : zeroOrder)
+ elementvector[local_i]
+ += operatorTerm->evalZot(iq, rowShapeValues[i])
+ * quad[iq].weight() * integrationElement;
+ }
}
}
@@ -142,12 +170,15 @@ namespace AMDiS
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleFirstOrderTermsGrdPhi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix)
+ ElementMatrix& elementMatrix,
+ double factor)
{
+ AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPhi(elementMatrix)");
+
using Element = typename RowView::Element;
auto element = rowView.element();
- const int dim = Element::dimension;
auto geometry = element.geometry();
+ const int dim = Element::dimension;
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
@@ -156,41 +187,41 @@ namespace AMDiS
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPhi)
- operatorTerm->init(element, quad);
+ operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
- // Position of the current quadrature point in the reference element
- const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
-
- // The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos);
-
- std::vector<Dune::FieldVector<double,1> > rowShapeValues;
- rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
-
- // The gradients of the shape functions on the reference element
- std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
- colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
-
- // Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,dim> > colGradients(colReferenceGradients.size());
-
- for (size_t i = 0; i < colGradients.size(); ++i)
- jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
-
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
- const int local_i = rowView.tree().localIndex(i);
- const int local_j = colView.tree().localIndex(j);
- for (auto* operatorTerm : firstOrderGrdPhi)
- elementMatrix[local_i][local_j]
- += operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j])
- * quad[iq].weight() * integrationElement;
- }
- }
+ // Position of the current quadrature point in the reference element
+ const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ // The multiplicative factor in the integral transformation formula
+ const double integrationElement = geometry.integrationElement(quadPos) * factor;
+
+ std::vector<Dune::FieldVector<double,1> > rowShapeValues;
+ rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
+
+ // The gradients of the shape functions on the reference element
+ std::vector<Dune::FieldMatrix<double,1,dim> > colReferenceGradients;
+ colLocalBasis.evaluateJacobian(quadPos, colReferenceGradients);
+
+ // Compute the shape function gradients on the real element
+ std::vector<Dune::FieldVector<double,dim> > colGradients(colReferenceGradients.size());
+
+ for (size_t i = 0; i < colGradients.size(); ++i)
+ jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
+
+ for (size_t i = 0; i < elementMatrix.N(); ++i) {
+ for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ const int local_i = rowView.tree().localIndex(i);
+ const int local_j = colView.tree().localIndex(j);
+ for (auto* operatorTerm : firstOrderGrdPhi)
+ elementMatrix[local_i][local_j]
+ += operatorTerm->evalFot1(iq, rowShapeValues[i], colGradients[j])
+ * quad[iq].weight() * integrationElement;
+ }
+ }
}
}
@@ -199,12 +230,15 @@ namespace AMDiS
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleFirstOrderTermsGrdPsi(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix)
+ ElementMatrix& elementMatrix,
+ double factor)
{
+ AMDIS_FUNCNAME("Operator::assembleFirstOrderTermsGrdPsi(elementMatrix)");
+
using Element = typename RowView::Element;
auto element = rowView.element();
- const int dim = Element::dimension;
auto geometry = element.geometry();
+ const int dim = Element::dimension;
const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
@@ -213,41 +247,41 @@ namespace AMDiS
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : firstOrderGrdPsi)
- operatorTerm->init(element, quad);
+ operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
- // Position of the current quadrature point in the reference element
- const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
-
- // The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos);
-
- // The gradients of the shape functions on the reference element
- std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
- rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
-
- // Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,dim> > rowGradients(rowReferenceGradients.size());
-
- for (size_t i = 0; i < rowGradients.size(); ++i)
- jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
-
- std::vector<Dune::FieldVector<double,1> > colShapeValues;
- colLocalBasis.evaluateFunction(quadPos, colShapeValues);
-
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
- const int local_i = rowView.tree().localIndex(i);
- const int local_j = colView.tree().localIndex(j);
- for (auto* operatorTerm : firstOrderGrdPsi)
- elementMatrix[local_i][local_j]
- += operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j])
- * quad[iq].weight() * integrationElement;
- }
- }
+ // Position of the current quadrature point in the reference element
+ const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ // The multiplicative factor in the integral transformation formula
+ const double integrationElement = geometry.integrationElement(quadPos) * factor;
+
+ // The gradients of the shape functions on the reference element
+ std::vector<Dune::FieldMatrix<double,1,dim> > rowReferenceGradients;
+ rowLocalBasis.evaluateJacobian(quadPos, rowReferenceGradients);
+
+ // Compute the shape function gradients on the real element
+ std::vector<Dune::FieldVector<double,dim> > rowGradients(rowReferenceGradients.size());
+
+ for (size_t i = 0; i < rowGradients.size(); ++i)
+ jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
+
+ std::vector<Dune::FieldVector<double,1> > colShapeValues;
+ colLocalBasis.evaluateFunction(quadPos, colShapeValues);
+
+ for (size_t i = 0; i < elementMatrix.N(); ++i) {
+ for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ const int local_i = rowView.tree().localIndex(i);
+ const int local_j = colView.tree().localIndex(j);
+ for (auto* operatorTerm : firstOrderGrdPsi)
+ elementMatrix[local_i][local_j]
+ += operatorTerm->evalFot2(iq, rowGradients[i], colShapeValues[j])
+ * quad[iq].weight() * integrationElement;
+ }
+ }
}
}
@@ -256,8 +290,11 @@ namespace AMDiS
template <class RowView, class ColView, class ElementMatrix>
void Operator<MeshView>::assembleSecondOrderTerms(RowView const& rowView,
ColView const& colView,
- ElementMatrix& elementMatrix)
+ ElementMatrix& elementMatrix,
+ double factor)
{
+ AMDIS_FUNCNAME("Operator::assembleSecondOrderTerms(elementMatrix)");
+
using Element = typename RowView::Element;
auto element = rowView.element();
const int dim = Element::dimension;
@@ -270,41 +307,41 @@ namespace AMDiS
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
for (auto* operatorTerm : secondOrder)
- operatorTerm->init(element, quad);
+ operatorTerm->init(element, quad);
// TODO: currently only the implementation for equal fespaces
assert( psiDegree == phiDegree );
for (size_t iq = 0; iq < quad.size(); ++iq) {
- // Position of the current quadrature point in the reference element
- const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
-
- // The multiplicative factor in the integral transformation formula
- const double integrationElement = geometry.integrationElement(quadPos);
-
- // The gradients of the shape functions on the reference element
- std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
- rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);
-
- // Compute the shape function gradients on the real element
- std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
-
- for (size_t i = 0; i < gradients.size(); ++i)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
-
- for (size_t i = 0; i < elementMatrix.N(); ++i) {
- for (size_t j = 0; j < elementMatrix.M(); ++j) {
- const int local_i = rowView.tree().localIndex(i);
- const int local_j = colView.tree().localIndex(j);
- for (auto* operatorTerm : secondOrder)
- elementMatrix[local_i][local_j]
- += operatorTerm->evalSot(iq, gradients[i], gradients[j])
- * quad[iq].weight() * integrationElement;
- }
- }
+ // Position of the current quadrature point in the reference element
+ const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
+
+ // The transposed inverse Jacobian of the map from the reference element to the element
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ // The multiplicative factor in the integral transformation formula
+ const double integrationElement = geometry.integrationElement(quadPos) * factor;
+
+ // The gradients of the shape functions on the reference element
+ std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
+ rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);
+
+ // Compute the shape function gradients on the real element
+ std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());
+
+ for (size_t i = 0; i < gradients.size(); ++i)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ for (size_t i = 0; i < elementMatrix.N(); ++i) {
+ for (size_t j = 0; j < elementMatrix.M(); ++j) {
+ const int local_i = rowView.tree().localIndex(i);
+ const int local_j = colView.tree().localIndex(j);
+ for (auto* operatorTerm : secondOrder)
+ elementMatrix[local_i][local_j]
+ += operatorTerm->evalSot(iq, gradients[i], gradients[j])
+ * quad[iq].weight() * integrationElement;
+ }
+ }
}
}
@@ -332,7 +369,7 @@ namespace AMDiS
int maxTermDegree = 0;
for (OperatorTermType* term : *terms)
- maxTermDegree = std::max(maxTermDegree, term->getDegree());
+ maxTermDegree = std::max(maxTermDegree, term->getDegree());
return psiDegree + phiDegree - order + maxTermDegree;
}
diff --git a/dune/amdis/OperatorTerm.hpp b/dune/amdis/OperatorTerm.hpp
index d67f824f..b7bfbc48 100644
--- a/dune/amdis/OperatorTerm.hpp
+++ b/dune/amdis/OperatorTerm.hpp
@@ -11,106 +11,108 @@
#include "OperatorTermBase.hpp"
-namespace AMDiS {
-
-template <class MeshView>
-class OperatorTerm
+namespace AMDiS
{
-protected:
- using Codim0 = typename MeshView::template Codim<0>;
- using Element = typename Codim0::Entity;
+
+ template <class MeshView>
+ class OperatorTerm
+ {
+ protected:
+ using Codim0 = typename MeshView::template Codim<0>;
+ using Element = typename Codim0::Entity;
static constexpr int dim = Element::dimension;
using QuadratureRule = Dune::QuadratureRule<double, dim>;
using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
-
-public:
+
+ public:
virtual void init(Element const& element,
- PointList const& points) = 0;
-
+ PointList const& points) = 0;
+
virtual double evalZot(size_t iq,
- Dune::FieldVector<double,1> const& test,
- Dune::FieldVector<double,1> const trial = 1.0) const = 0;
-
+ Dune::FieldVector<double,1> const& test,
+ Dune::FieldVector<double,1> const trial = 1.0) const = 0;
+
virtual double evalFot1(size_t iq,
- Dune::FieldVector<double,1> const& test,
- Dune::FieldVector<double,dim> const& grad_trial) const = 0;
-
+ Dune::FieldVector<double,1> const& test,
+ Dune::FieldVector<double,dim> const& grad_trial) const = 0;
+
virtual double evalFot2(size_t iq,
- Dune::FieldVector<double,dim> const& grad_test,
- Dune::FieldVector<double,1> const trial = 1.0) const = 0;
-
+ Dune::FieldVector<double,dim> const& grad_test,
+ Dune::FieldVector<double,1> const trial = 1.0) const = 0;
+
virtual double evalSot(size_t iq,
- Dune::FieldVector<double,dim> const& grad_test,
- Dune::FieldVector<double,dim> const& grad_trial) const = 0;
-
+ Dune::FieldVector<double,dim> const& grad_test,
+ Dune::FieldVector<double,dim> const& grad_trial) const = 0;
+
virtual int getDegree() const = 0;
-};
+ };
-template <class MeshView, class Term, class Traits = _none>
-class GenericOperatorTerm
- : public OperatorTerm<MeshView>
- , public OperatorEvaluation
-{
+ /// Base class for all operators based on expressions
+ template <class MeshView, class Term, class Traits = _none>
+ class GenericOperatorTerm
+ : public OperatorTerm<MeshView>
+ , public OperatorEvaluation
+ {
using Super = OperatorTerm<MeshView>;
using Element = typename Super::Element;
using PointList = typename Super::PointList;
static constexpr int dim = Element::dimension;
-public:
+ public:
GenericOperatorTerm(Term const& term)
- : term(term)
+ : term(term)
{}
virtual void init(Element const& element,
- PointList const& points) override
+ PointList const& points) override
{
- term.init(element, points);
-
- // cache term evaluation
- values.resize(points.size());
- for (size_t iq = 0; iq < points.size(); ++iq)
- values[iq] = term[iq];
+ term.init(element, points);
+
+ // cache term evaluation
+ values.resize(points.size());
+ for (size_t iq = 0; iq < points.size(); ++iq)
+ values[iq] = term[iq];
}
-
+
virtual double evalZot(size_t iq,
- Dune::FieldVector<double,1> const& test,
- Dune::FieldVector<double,1> const trial = 1.0) const override
+ Dune::FieldVector<double,1> const& test,
+ Dune::FieldVector<double,1> const trial = 1.0) const override
{
- return this->evalZotImpl(_cat{}, _traits{}, values[iq], test, trial);
+ return this->evalZotImpl(_cat{}, _traits{}, values[iq], test, trial);
}
-
+
virtual double evalFot1(size_t iq,
- Dune::FieldVector<double,1> const& test,
- Dune::FieldVector<double,dim> const& grad_trial) const override
+ Dune::FieldVector<double,1> const& test,
+ Dune::FieldVector<double,dim> const& grad_trial) const override
{
- return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
+ return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_trial, test);
}
-
+
virtual double evalFot2(size_t iq,
- Dune::FieldVector<double,dim> const& grad_test,
- Dune::FieldVector<double,1> const trial = 1.0) const override
+ Dune::FieldVector<double,dim> const& grad_test,
+ Dune::FieldVector<double,1> const trial = 1.0) const override
{
- return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
+ return this->evalFotImpl(_cat{}, _traits{}, values[iq], grad_test, trial);
}
-
+
virtual double evalSot(size_t iq,
- Dune::FieldVector<double,dim> const& grad_test,
- Dune::FieldVector<double,dim> const& grad_trial) const override
+ Dune::FieldVector<double,dim> const& grad_test,
+ Dune::FieldVector<double,dim> const& grad_trial) const override
{
- return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
+ return this->evalSotImpl(_cat{}, _traits{}, values[iq], grad_test, grad_trial);
}
-
+
virtual int getDegree() const override
{
- return term.getDegree();
+ return term.getDegree();
}
-private:
+ private:
Term term;
using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
@@ -118,17 +120,17 @@ private:
using _traits = Traits;
std::vector<value_type> values;
-};
+ };
// some example terms
// -----------------------------------------------------------------------------
-
-template <class ValueType>
-class ConstantTerm
-{
-public:
+ /// An expression representing a constant (arithmetic) value
+ template <class ValueType>
+ class ConstantTerm
+ {
+ public:
using value_type = ValueType;
ConstantTerm(value_type value)
@@ -137,34 +139,38 @@ public:
template <class Element, class PointList>
void init(Element const& element,
- PointList const& points) {}
+ PointList const& points) { /* do nothing */ }
value_type operator[](size_t iq) const
{
- return value;
+ return value;
}
int getDegree() const
{
- return 0;
+ return 0;
}
-
-private:
+
+ private:
value_type value;
-};
-
+ };
+
+
+ /// generator function for \ref ConstantTerm expressions
+ template <class T>
+ std::enable_if_t< std::is_arithmetic<T>::value, ConstantTerm<T> >
+ constant(T value) { return {value}; }
-// generator function for coordinate expressions
-template <class T>
-std::enable_if_t< std::is_arithmetic<T>::value, ConstantTerm<T> >
-constant(T value) { return {value}; }
+// -----------------------------------------------------------------------------
-template <class Functor>
-class CoordsTerm
-{
-public:
+ /// An expression that evaluates to the current coordinate of a dof or
+ /// quadrature point with given index.
+ template <class Functor>
+ class CoordsTerm
+ {
+ public:
using value_type = typename std::result_of<Functor(Dune::FieldVector<double, 2>)>::type;
template <class F,
@@ -176,52 +182,55 @@ public:
template <class Element, class PointList>
void init(Element const& element,
- PointList const& points)
+ PointList const& points)
{
- values.resize(points.size());
- for (size_t iq = 0; iq < points.size(); ++iq)
- values[iq] = fct(element.geometry().global(points[iq].position()));
+ AMDIS_FUNCNAME("CoordsTerm::init()");
+ values.resize(points.size());
+ for (size_t iq = 0; iq < points.size(); ++iq)
+ values[iq] = fct(element.geometry().global(points[iq].position()));
}
value_type const& operator[](size_t iq) const
{
- return values[iq];
+ return values[iq];
}
int getDegree() const
{
- return degree;
+ return degree;
}
-
-private:
+
+ private:
Functor fct;
int degree;
std::vector<value_type> values;
-};
+ };
-// generator function for coordinate expressions
-template <class F>
-CoordsTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
+ /// generator function for \ref CoordsTerm expressions
+ template <class F>
+ CoordsTerm< std::decay_t<F> > eval(F&& f) { return {std::forward<F>(f)}; }
-template <class FeSpace> struct GetDegree : int_<1> {};
-template <class GV, int k, class ST>
-struct GetDegree<Dune::Functions::PQkNodalBasis<GV, k, ST> > : int_<k> {};
+// -----------------------------------------------------------------------------
-template <class DOFVectorType>
-class DOFVectorTerm
-{
- using Basis = typename DOFVectorType::FeSpace;
- using field_type = typename DOFVectorType::field_type;
+ // helper class to extract the polynomial degree of a pqk nodal basis
+ template <class FeSpace> struct GetDegree : int_<1> {};
+ template <class GV, int k, class ST>
+ struct GetDegree<Dune::Functions::PQkNodalBasis<GV, k, ST> > : int_<k> {};
-public:
+
+ /// An expression that evalues a DOFVector at a given element, either on the
+ /// dofs or on the quadrature points
+ template <class DOFVectorType>
+ class DOFVectorTerm
+ {
+ public:
using value_type = typename DOFVectorType::value_type;
- template <class DOFVectorType_>
- DOFVectorTerm(DOFVectorType_&& dofvector, double factor = 1.0)
+ DOFVectorTerm(DOFVectorType const& dofvector, double factor = 1.0)
: vector(dofvector.getVector())
, factor(factor)
, localView(dofvector.getFeSpace().localView())
@@ -230,69 +239,80 @@ public:
template <class Element, class PointList>
void init(Element const& element,
- PointList const& points)
+ PointList const& points)
{
- localView.bind(element);
- localIndexSet.bind(localView);
-
- const auto& localFiniteElem = localView.tree().finiteElement();
- const size_t nBasisFct = localFiniteElem.size();
-
- std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
- std::vector<value_type> localVec(nBasisFct);
-
- for (size_t j = 0; j < nBasisFct; ++j) {
- const auto global_idx = localIndexSet.index(j);
- localVec[j] = vector[global_idx];
- }
-
- values.resize(points.size());
- for (size_t iq = 0; iq < points.size(); ++iq) {
- localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
- value_type result = 0;
- for (size_t j = 0; j < shapeValues.size(); ++j)
- result += localVec[j] * (factor * shapeValues[j]);
- values[iq] = result;
- }
+ AMDIS_FUNCNAME("DOFVectorTerm::init()");
+ localView.bind(element);
+ localIndexSet.bind(localView);
+
+ const auto& localFiniteElem = localView.tree().finiteElement();
+ const size_t nBasisFct = localFiniteElem.size();
+
+ std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
+ std::vector<value_type> localVec(nBasisFct);
+
+ for (size_t j = 0; j < nBasisFct; ++j) {
+ const auto global_idx = localIndexSet.index(j);
+ localVec[j] = factor * vector[global_idx];
+ }
+
+ values.resize(points.size());
+ for (size_t iq = 0; iq < points.size(); ++iq) {
+ localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
+ value_type result = 0;
+ for (size_t j = 0; j < shapeValues.size(); ++j)
+ result += localVec[j] * shapeValues[j];
+ values[iq] = result;
+ }
}
value_type const& operator[](size_t iq) const
{
- return values[iq];
+ return values[iq];
}
int getDegree() const
{
- return degree;
+ return degree;
}
-
-private:
+
+ private:
typename DOFVectorType::BaseVector const& vector;
double factor;
+ using Basis = typename DOFVectorType::FeSpace;
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
int degree = GetDegree<Basis>::value;
std::vector<value_type> values;
-};
+ };
+
+ /// generator function for \ref DOFVector expressions
+ template <class DOFVectorType>
+ DOFVectorTerm<std::decay_t<DOFVectorType>>
+ valueOf(DOFVectorType&& vector, double factor = 1.0)
+ {
+ return {std::forward<DOFVectorType>(vector), factor};
+ }
+// -----------------------------------------------------------------------------
-template <class DOFVectorType, class F>
-class DOFVectorFuncTerm
-{
- using Basis = typename DOFVectorType::FeSpace;
- using field_type = typename DOFVectorType::field_type;
-
-public:
+
+ /// An expression that evalues a DOFVector at a given element, either on the
+ /// dofs or on the quadrature points, and applies a functor to the value
+ template <class DOFVectorType, class Func>
+ class DOFVectorFuncTerm
+ {
+ public:
using value_type = typename DOFVectorType::value_type;
- template <class DOFVectorType_, class F_>
- DOFVectorFuncTerm(DOFVectorType_&& dofvector, F_&& f, int f_deg)
+ template <class F_>
+ DOFVectorFuncTerm(DOFVectorType const& dofvector, F_&& f, int f_deg)
: vector(dofvector.getVector())
- , f(f)
+ , f(std::forward<F_>(f))
, localView(dofvector.getFeSpace().localView())
, localIndexSet(dofvector.getFeSpace().localIndexSet())
, f_deg(f_deg)
@@ -300,47 +320,49 @@ public:
template <class Element, class PointList>
void init(Element const& element,
- PointList const& points)
+ PointList const& points)
{
- localView.bind(element);
- localIndexSet.bind(localView);
-
- const auto& localFiniteElem = localView.tree().finiteElement();
- const size_t nBasisFct = localFiniteElem.size();
-
- std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
- std::vector<value_type> localVec(nBasisFct);
-
- for (size_t j = 0; j < nBasisFct; ++j) {
- const auto global_idx = localIndexSet.index(j);
- localVec[j] = vector[global_idx];
- }
-
- values.resize(points.size());
- for (size_t iq = 0; iq < points.size(); ++iq) {
- localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
- value_type result = 0;
- for (size_t j = 0; j < shapeValues.size(); ++j)
- result += localVec[j] * shapeValues[j];
-
- values[iq] = f(result);
- }
+ AMDIS_FUNCNAME("DOFVectorFuncTerm::init()");
+ localView.bind(element);
+ localIndexSet.bind(localView);
+
+ const auto& localFiniteElem = localView.tree().finiteElement();
+ const size_t nBasisFct = localFiniteElem.size();
+
+ std::vector<Dune::FieldVector<double,1> > shapeValues(nBasisFct);
+ std::vector<value_type> localVec(nBasisFct);
+
+ for (size_t j = 0; j < nBasisFct; ++j) {
+ const auto global_idx = localIndexSet.index(j);
+ localVec[j] = vector[global_idx];
+ }
+
+ values.resize(points.size());
+ for (size_t iq = 0; iq < points.size(); ++iq) {
+ localFiniteElem.localBasis().evaluateFunction(points[iq].position(), shapeValues);
+ value_type result = 0;
+ for (size_t j = 0; j < shapeValues.size(); ++j)
+ result += localVec[j] * shapeValues[j];
+
+ values[iq] = f(result);
+ }
}
value_type const& operator[](size_t iq) const
{
- return values[iq];
+ return values[iq];
}
int getDegree() const
{
- return degree * f_deg;
+ return degree * f_deg;
}
-
-private:
+
+ private:
typename DOFVectorType::BaseVector const& vector;
- F f;
+ Func f;
+ using Basis = typename DOFVectorType::FeSpace;
typename Basis::LocalView localView;
typename Basis::LocalIndexSet localIndexSet;
@@ -348,23 +370,16 @@ private:
int f_deg;
std::vector<value_type> values;
-};
-
-
-template <class DOFVectorType>
-DOFVectorTerm<std::decay_t<DOFVectorType>>
-valueOf(DOFVectorType&& vector, double factor = 1.0)
-{
- return {std::forward<DOFVectorType>(vector), factor};
-}
+ };
-
-template <class DOFVectorType, class F>
-DOFVectorFuncTerm<std::decay_t<DOFVectorType>, std::decay_t<F>>
-valueOfFunc(DOFVectorType&& vector, F&& f, int deg = 1)
-{
+
+ /// generator function for \ref DOFVectorFuncTerm expressions
+ template <class DOFVectorType, class F>
+ DOFVectorFuncTerm<std::decay_t<DOFVectorType>, std::decay_t<F>>
+ valueOfFunc(DOFVectorType&& vector, F&& f, int deg = 1)
+ {
return {std::forward<DOFVectorType>(vector), std::forward<F>(f), deg};
-}
+ }
} // end namespace AMDiS
diff --git a/dune/amdis/OperatorTermBase.hpp b/dune/amdis/OperatorTermBase.hpp
index a0e4e521..c402cbcc 100644
--- a/dune/amdis/OperatorTermBase.hpp
+++ b/dune/amdis/OperatorTermBase.hpp
@@ -8,77 +8,77 @@ namespace AMDiS {
template <class K, int SIZE>
K sum(Dune::FieldVector<K, SIZE> const& vector)
{
- K result = 0;
- for (int i = 0; i < SIZE; ++i)
- result += vector[i];
- return result;
+ K result = 0;
+ for (int i = 0; i < SIZE; ++i)
+ result += vector[i];
+ return result;
}
class OperatorEvaluation
{
- using Scalar = Dune::FieldVector<double, 1>;
+ using Scalar = Dune::FieldVector<double, 1>;
protected:
- template <class... Args>
- auto evalZotImpl(Args...) const { return 0; }
-
- template <class ValueType>
- auto evalZotImpl(_scalar, _none, ValueType const& v,
- Scalar const& test, Scalar const& trial) const
- {
- return v * test * trial;
- }
-
- template <class... Args>
- auto evalFotImpl(Args...) const { return 0; }
-
- template <size_t I, class ValueType, class Gradient>
- auto evalFotImpl(_scalar, VectorComponent<I>, ValueType const& v,
- Gradient const& grad_test, Scalar const& trial) const
- {
- return v * grad_test[I] * trial;
- }
-
- template <class ValueType, class Gradient>
- auto evalFotImpl(_vector, _none, ValueType const& v,
- Gradient const& grad_test, Scalar const& trial) const
- {
- return (v * grad_test) * trial;
- }
-
- template <class ValueType, class Gradient>
- auto evalFotImpl(_scalar, _none, ValueType const& v,
- Gradient const& grad_test, Scalar const& trial) const
- {
- return v * sum(grad_test) * trial;
- }
-
-
- template <class... Args>
- auto evalSotImpl(Args...) const { return 0; }
-
- template <size_t R, size_t C, class ValueType, class Gradient>
- auto evalSotImpl(_scalar, MatrixComponent<R, C>, ValueType const& v,
- Gradient const& grad_test, Gradient const& grad_trial) const
- {
- return v * (grad_test[R] * grad_trial[C]);
- }
-
- template <class ValueType, class Gradient>
- auto evalSotImpl(_matrix, _none, ValueType const& v,
- Gradient const& grad_test, Gradient const& grad_trial) const
- {
- return grad_test * (v * grad_trial);
- }
-
- template <class ValueType, class Gradient>
- auto evalSotImpl(_scalar, _none, ValueType const& v,
- Gradient const& grad_test, Gradient const& grad_trial) const
- {
- return v * (grad_test * grad_trial);
- }
+ template <class... Args>
+ auto evalZotImpl(Args...) const { assert( false ); return 0; }
+
+ template <class ValueType>
+ auto evalZotImpl(_scalar, _none, ValueType const& v,
+ Scalar const& test, Scalar const& trial) const
+ {
+ return v * test * trial;
+ }
+
+ template <class... Args>
+ auto evalFotImpl(Args...) const { assert( false ); return 0; }
+
+ template <size_t I, class ValueType, class Gradient>
+ auto evalFotImpl(_scalar, VectorComponent<I>, ValueType const& v,
+ Gradient const& grad_test, Scalar const& trial) const
+ {
+ return v * grad_test[I] * trial;
+ }
+
+ template <class ValueType, class Gradient>
+ auto evalFotImpl(_vector, _none, ValueType const& v,
+ Gradient const& grad_test, Scalar const& trial) const
+ {
+ return (v * grad_test) * trial;
+ }
+
+ template <class ValueType, class Gradient>
+ auto evalFotImpl(_scalar, _none, ValueType const& v,
+ Gradient const& grad_test, Scalar const& trial) const
+ {
+ return v * sum(grad_test) * trial;
+ }
+
+
+ template <class... Args>
+ auto evalSotImpl(Args...) const { assert( false ); return 0; }
+
+ template <size_t R, size_t C, class ValueType, class Gradient>
+ auto evalSotImpl(_scalar, MatrixComponent<R, C>, ValueType const& v,
+ Gradient const& grad_test, Gradient const& grad_trial) const
+ {
+ return v * (grad_test[R] * grad_trial[C]);
+ }
+
+ template <class ValueType, class Gradient>
+ auto evalSotImpl(_matrix, _none, ValueType const& M,
+ Gradient const& grad_test, Gradient const& grad_trial) const
+ {
+ return grad_test * (M * grad_trial);
+ }
+
+ template <class ValueType, class Gradient>
+ auto evalSotImpl(_scalar, _none, ValueType const& v,
+ Gradient const& grad_test, Gradient const& grad_trial) const
+ {
+ return v * (grad_test * grad_trial);
+ }
};
diff --git a/dune/amdis/ProblemInstat.hpp b/dune/amdis/ProblemInstat.hpp
index 5720660a..cb3a18df 100644
--- a/dune/amdis/ProblemInstat.hpp
+++ b/dune/amdis/ProblemInstat.hpp
@@ -40,7 +40,7 @@ namespace AMDiS
{}
/// Initialisation of the problem.
- virtual void initialize(Flag initFlag);
+ virtual void initialize(Flag initFlag = INIT_NOTHING);
/// Used in \ref initialize().
virtual void createUhOld();
diff --git a/dune/amdis/ProblemStat.hpp b/dune/amdis/ProblemStat.hpp
index 34579b9c..56ea8013 100644
--- a/dune/amdis/ProblemStat.hpp
+++ b/dune/amdis/ProblemStat.hpp
@@ -29,7 +29,8 @@ namespace AMDiS
class ProblemStatSeq : public ProblemStatBase
{
using Self = ProblemStatSeq;
- public:
+
+ public: // typedefs and static constants
using FeSpaces = typename Traits::FeSpaces;
using Mesh = typename Traits::Mesh;
@@ -65,11 +66,12 @@ namespace AMDiS
typename SystemVectorType::MultiVector>;
public:
- /// \brief Constructor. Takes the name of the problem that is used to
- /// access values correpsonding to this püroblem in the parameter file.
- /**
+ /**
+ * \brief Constructor. Takes the name of the problem that is used to
+ * access values correpsonding to this püroblem in the parameter file.
+ *
* Parameters read by ProblemStatSeq, with name 'PROB'
- * PROB->names: \ref componentNames
+ * PROB->names: \ref componentNames
**/
ProblemStatSeq(std::string name)
: name(name)
@@ -80,8 +82,9 @@ namespace AMDiS
}
- /// \brief Initialisation of the problem.
/**
+ * \brief Initialisation of the problem.
+ *
* Parameters read in initialize() for problem with name 'PROB'
* MESH[0]->global refinements: nr of initial global refinements
**/
@@ -125,10 +128,9 @@ namespace AMDiS
/// Writes output files.
- void writeFiles(AdaptInfo& adaptInfo, bool force);
+ void writeFiles(AdaptInfo& adaptInfo, bool force = false);
- // get-methods
- // -------------------------------------------------------------------------
+ public: // get-methods
/// Returns nr of components \ref nComponents
static constexpr size_t getNumComponents()
@@ -136,10 +138,12 @@ namespace AMDiS
return nComponents;
}
+
/// Return the \ref systemMatrix
auto getSystemMatrix() { return systemMatrix; }
auto getSystemMatrix() const { return systemMatrix; }
+
/// Return the \ref solution
auto getSolution() { return solution; }
auto getSolution() const { return solution; }
@@ -151,19 +155,24 @@ namespace AMDiS
return (*solution)[_i];
}
+
/// Return the \ref rhs
auto getRhs() { return rhs; }
auto getRhs() const { return rhs; }
+
/// Return the \ref linearSolver
auto getSolver() { return linearSolver; }
/// Return the \ref mesh
auto getMesh() { return mesh; }
+ auto getMeshView() { return meshView; }
+
/// Return the \ref feSpaces
auto getFeSpaces() { return feSpaces; }
+
/// Return the I'th \ref feSpaces component
template <size_t I = 0>
FeSpace<I> const& getFeSpace(const index_<I> = index_<I>()) const
@@ -171,22 +180,21 @@ namespace AMDiS
return std::get<I>(*feSpaces);
}
- /// Return the \ref name of the problem. Implementation of \ref ProblemStatBase::getName
+ /// Implementation of \ref ProblemStatBase::getName
virtual std::string getName() const override
{
return name;
}
- protected:
+ std::vector<std::string> getComponentNames() const
+ {
+ return componentNames;
+ }
- // initialization methods
- // -------------------------------------------------------------------------
+ protected: // initialization methods
- // Reads a filename from the init file and creates the mesh object with
- // respect to this file.
void createMesh()
{
- // TODO: Creation from meshname must be generalized to other meshes than AlbertaGrid.
meshName = "";
Parameters::get(name + "->mesh", meshName);
AMDIS_TEST_EXIT(!meshName.empty(),
@@ -194,23 +202,28 @@ namespace AMDiS
mesh = MeshCreator<Mesh>::create(meshName);
meshView = std::make_shared<MeshView>(mesh->leafGridView());
+
+ AMDIS_MSG("Create mesh:");
+ AMDIS_MSG("#elements = " << mesh->size(0));
+ AMDIS_MSG("#faces/edges = " << mesh->size(1));
+ AMDIS_MSG("#vertices = " << mesh->size(dim));
+ AMDIS_MSG("");
}
- //
void createFeSpaces()
{
feSpaces = std::make_shared<FeSpaces>(construct_tuple<FeSpaces>(*meshView));
}
- //
void createMatricesAndVectors()
{
systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
- rhs = std::make_shared<SystemVectorType>(*feSpaces, std::vector<std::string>(nComponents, "rhs"));
+
+ auto rhsNames = std::vector<std::string>(nComponents, "rhs");
+ rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
}
- //
void createSolver()
{
using Creator = SolverCreator<typename SystemMatrixType::MultiMatrix,
@@ -221,33 +234,32 @@ namespace AMDiS
linearSolver = Creator::create(solverName, name + "->solver");
}
- //
void createFileWriter()
{
- filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", *meshView, componentNames);
+ filewriter = std::make_shared<FileWriter<Traits>>(name + "->output",
+ *meshView,
+ componentNames);
}
- // sub-methods to assemble DOFMatrix
- // -------------------------------------------------------------------------
+ protected: // sub-methods to assemble DOFMatrix
- //
template <class Matrix, class Vector>
void assemble(std::pair<int, int> row_col,
Matrix* matrix,
Vector* rhs);
- //
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
- std::list<std::shared_ptr<OperatorType>>& operators);
+ std::list<std::shared_ptr<OperatorType>>& operators,
+ std::list<double*> const& factors);
- //
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
- std::list<std::shared_ptr<OperatorType>>& operators);
+ std::list<std::shared_ptr<OperatorType>>& operators,
+ std::list<double*> const& factors);
template <class Matrix, class RowIndexSet, class ColIndexSet>
@@ -261,7 +273,8 @@ namespace AMDiS
IndexSet const& indexSet,
ElementVector const& elementvector);
- private:
+ private: // some internal methods
+
template <size_t I = 0>
typename FeSpace<I>::NodeFactory&
getNodeFactory(const index_<I> = index_<I>())
@@ -279,11 +292,12 @@ namespace AMDiS
std::vector<std::string> componentNames;
/// Mesh of this problem.
- // TODO: generalize to multi-mesh problems
- std::shared_ptr<Mesh> mesh;
+ std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
+
+ /// Name of the mesh
std::string meshName;
- ///
+ /// A gridView object
std::shared_ptr<MeshView> meshView;
/// Pointer to the meshes for the different problem components
@@ -292,16 +306,35 @@ namespace AMDiS
/// FE spaces of this problem.
std::shared_ptr<FeSpaces> feSpaces; // eventuell const
+ /// A FileWriter object
std::shared_ptr<FileWriter<Traits>> filewriter;
+
+ /// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver;
+ /// A block-matrix that is filled during assembling
std::shared_ptr<SystemMatrixType> systemMatrix;
+
+ /// A block-vector with the solution components
std::shared_ptr<SystemVectorType> solution;
+
+ /// A block-vector (load-vector) corresponding to the right.hand side
+ /// of the equation, filled during assembling
std::shared_ptr<SystemVectorType> rhs;
- IdxPairList<DirichletBC<WorldVector>> dirichletBc;
- IdxPairList<OperatorType> matrixOperators;
- IdxList<OperatorType> vectorOperators;
+ /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
+ /// each matrix block
+ IdxPairList<DirichletBC<WorldVector>> dirichletBc;
+
+ /// A map (i,j) -> list<OperatorType> string the differential operators for
+ /// each matrix block
+ IdxPairList<OperatorType> matrixOperators;
+ std::map<std::pair<int,int>, std::list<double*>> matrixFactors;
+
+ /// A map (i) -> list<OperatorType> string the differential operators for
+ /// each rhs block
+ IdxList<OperatorType> vectorOperators;
+ std::map<int, std::list<double*>> vectorFactors;
};
@@ -312,10 +345,7 @@ namespace AMDiS
// extern template class ProblemStatSeq<ProblemStatTraits<3>>;
#endif
-
-
-
- namespace detail
+ namespace Impl
{
template <class ProblemStatType>
struct ProblemStat : public ProblemStatType,
@@ -325,14 +355,19 @@ namespace AMDiS
/// Constructor
explicit ProblemStat(std::string nameStr)
- : ProblemStatType(nameStr),
- StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
+ : ProblemStatType(nameStr)
+ , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
{}
- /// Determines the execution order of the single adaption steps. If adapt is
- /// true, mesh adaption will be performed. This allows to avoid mesh adaption,
- /// e.g. in timestep adaption loops of timestep adaptive strategies.
- // implements StandardProblemIteration::oneIteration(AdaptInfo&, Flag)
+ /**
+ * \brief Determines the execution order of the single adaption steps.
+ *
+ * If adapt is true, mesh adaption will be performed. This allows to avoid
+ * mesh adaption, e.g. in timestep adaption loops of timestep adaptive
+ * strategies.
+ *
+ * Implementation of \ref StandardProblemIteration::oneIteration.
+ **/
virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
{
for (int i = 0; i < ProblemStatType::getNumComponents(); i++)
@@ -344,32 +379,31 @@ namespace AMDiS
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
- /// Implementation of \ref ProblemStatBase::buildBeforeRefine().
- virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::buildBeforeRefine.
+ virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }
- /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen().
- virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
+ virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }
- /// Implementation of \ref ProblemStatBase::estimate().
- virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::estimate.
+ virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
- /// Implementation of \ref ProblemStatBase::refineMesh().
- virtual Flag refineMesh(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::refineMesh.
+ virtual Flag refineMesh(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
- /// Implementation of \ref ProblemStatBase::coarsenMesh().
- virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::coarsenMesh.
+ virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
- /// Implementation of \ref ProblemStatBase::markElements().
- virtual Flag markElements(AdaptInfo& adaptInfo) override { /* Does nothing here. */ }
+ /// Implementation of \ref ProblemStatBase::markElements.
+ virtual Flag markElements(AdaptInfo& adaptInfo) override { /* Does nothing. */ }
};
- }
+
+ } // end namespace Impl
template <class Traits>
- using ProblemStat = detail::ProblemStat<ProblemStatSeq<Traits>>;
-
-
+ using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;
-}
+} // end namespace AMDiS
#include "ProblemStat.inc.hpp"
diff --git a/dune/amdis/ProblemStat.inc.hpp b/dune/amdis/ProblemStat.inc.hpp
index 0e2c0b91..711b5d71 100644
--- a/dune/amdis/ProblemStat.inc.hpp
+++ b/dune/amdis/ProblemStat.inc.hpp
@@ -18,16 +18,14 @@ namespace AMDiS
else {
if (initFlag.isSet(CREATE_MESH) ||
(!adoptFlag.isSet(INIT_MESH) &&
- (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
- {
+ (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
createMesh();
}
if (adoptProblem &&
(adoptFlag.isSet(INIT_MESH) ||
adoptFlag.isSet(INIT_SYSTEM) ||
- adoptFlag.isSet(INIT_FE_SPACE)))
- {
+ adoptFlag.isSet(INIT_FE_SPACE))) {
mesh = adoptProblem->getMesh();
}
@@ -40,7 +38,16 @@ namespace AMDiS
int globalRefinements = 0;
Parameters::get(meshName + "->global refinements", globalRefinements);
- mesh->globalRefine(globalRefinements);
+ if (globalRefinements > 0) {
+ mesh->globalRefine(globalRefinements);
+ AMDIS_MSG("After refinement:");
+ AMDIS_MSG("#elements = " << mesh->size(0));
+ AMDIS_MSG("#faces/edges = " << mesh->size(1));
+ AMDIS_MSG("#vertices = " << mesh->size(dim));
+ AMDIS_MSG("");
+ }
+
+
// === create fespace ===
if (feSpaces) {
@@ -48,14 +55,12 @@ namespace AMDiS
}
else {
if (initFlag.isSet(INIT_FE_SPACE) ||
- (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
- {
+ (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
createFeSpaces();
}
if (adoptProblem &&
- (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM)))
- {
+ (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
feSpaces = adoptProblem->getFeSpaces();
}
}
@@ -63,37 +68,35 @@ namespace AMDiS
if (!feSpaces)
AMDIS_WARNING("no feSpaces created\n");
+ AMDIS_MSG("#DOFs = " << this->getFeSpace<0>().size());
// === create system ===
if (initFlag.isSet(INIT_SYSTEM))
createMatricesAndVectors();
- if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM))
- {
+ if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->getSolution();
rhs = adoptProblem->getRhs();
systemMatrix = adoptProblem->getSystemMatrix();
}
// === create solver ===
- if (linearSolver)
- {
+ if (linearSolver) {
AMDIS_WARNING("solver already created\n");
}
- else
- {
+ else {
if (initFlag.isSet(INIT_SOLVER))
createSolver();
- if (adoptProblem && adoptFlag.isSet(INIT_SOLVER))
- {
+ if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
AMDIS_TEST_EXIT(!linearSolver, "solver already created\n");
linearSolver = adoptProblem->getSolver();
}
}
- if (!linearSolver)
+ if (!linearSolver) {
AMDIS_WARNING("no solver created\n");
+ }
// === create file writer ===
if (initFlag.isSet(INIT_FILEWRITER))
@@ -110,9 +113,10 @@ namespace AMDiS
double* estFactor)
{
// TODO: currently the factors are ignored
- assert( factor == NULL && estFactor == NULL );
+ assert( estFactor == NULL );
- matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
+ matrixOperators[{i,j}].emplace_back(new OperatorType{op});
+ matrixFactors[{i,j}].push_back(factor);
}
@@ -122,10 +126,11 @@ namespace AMDiS
double* factor,
double* estFactor)
{
- // TODO: currently the factors are ignored
- assert( factor == NULL && estFactor == NULL );
+ // TODO: currently the est-factors are ignored
+ assert( estFactor == NULL );
- vectorOperators[i].push_back(std::make_shared<OperatorType>(op));
+ vectorOperators[i].emplace_back(new OperatorType{op});
+ vectorFactors[i].push_back(factor);
}
@@ -148,10 +153,8 @@ namespace AMDiS
// boundaryConditionSet = true;
- using BcType = DirichletBC<WorldVector>;
- std::shared_ptr<BcType> dirichlet = std::make_shared<BcType>(predicate, values);
-
- dirichletBc[std::make_pair(row, col)].push_back( dirichlet );
+ using BcType = DirichletBC<WorldVector>;
+ dirichletBc[{row, col}].emplace_back( new BcType{predicate, values} );
}
@@ -193,7 +196,6 @@ namespace AMDiS
tol_str << "relTol = " << solverInfo.getRelTolerance() << " ";
AMDIS_ERROR_EXIT("Tolerance " << tol_str.str() << " could not be reached!");
}
-
}
@@ -203,8 +205,7 @@ namespace AMDiS
{
Timer t;
- For<0, nComponents>::loop([this](auto const _r)
- {
+ For<0, nComponents>::loop([this](auto const _r) {
this->getNodeFactory(_r).initializeIndices();
});
@@ -212,7 +213,7 @@ namespace AMDiS
size_t nnz = 0;
For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector](auto const _r)
{
- AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for " << "FeSpace[" << _r << "]");
+ AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
if (asmVector) {
rhs->compress(_r);
@@ -224,16 +225,13 @@ namespace AMDiS
auto const& rowFeSpace = this->getFeSpace(_r);
auto const& colFeSpace = this->getFeSpace(_c);
- int r = int(_r), c = int(_c);
- auto row_col = std::make_pair(r, c);
-
// The DOFMatrix which should be assembled
auto& dofmatrix = (*systemMatrix)(_r, _c);
- auto& solution_vec = solution->getDOFVector(_c);
- auto& rhs_vec = rhs->getDOFVector(_r);
+ auto& solution_vec = (*solution)[_c];
+ auto& rhs_vec = (*rhs)[_r];
+ int r = 0, c = 0;
if (asmMatrix) {
- dofmatrix.init();
// init boundary condition
for (auto bc_list : dirichletBc) {
@@ -245,7 +243,7 @@ namespace AMDiS
}
// assemble the DOFMatrix block and the corresponding rhs vector, of r==c
- this->assemble(row_col, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
+ this->assemble({int(_r), int(_c)}, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));
// finish boundary condition
for (auto bc_list : dirichletBc) {
@@ -255,7 +253,6 @@ namespace AMDiS
bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
}
}
- dofmatrix.finish();
nnz += dofmatrix.getNnz();
}
@@ -286,9 +283,15 @@ namespace AMDiS
auto const& rowFeSpace = dofmatrix->getRowFeSpace();
auto const& colFeSpace = dofmatrix->getColFeSpace();
- for (auto op : matrixOperators[row_col])
+ dofmatrix->init();
+ auto matrixOp = matrixOperators[row_col];
+ auto matrixOpFac = matrixFactors[row_col];
+ auto vectorOp = vectorOperators[row_col.first];
+ auto vectorOpFac = vectorFactors[row_col.first];
+
+ for (auto op : matrixOp)
op->init(rowFeSpace, colFeSpace);
- for (auto op : vectorOperators[row_col.first])
+ for (auto op : vectorOp)
op->init(rowFeSpace, colFeSpace);
auto rowLocalView = rowFeSpace.localView();
@@ -306,18 +309,21 @@ namespace AMDiS
if (dofmatrix) {
ElementMatrix elementMatrix;
- bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, matrixOperators[row_col]);
+ bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
+ matrixOp, matrixOpFac);
if (add)
addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
}
if (rhs) {
ElementVector elementVector;
- bool add = getElementVector(rowLocalView, elementVector, vectorOperators[row_col.first]);
+ bool add = getElementVector(rowLocalView, elementVector,
+ vectorOp, vectorOpFac);
if (add)
addElementVector(*rhs, rowIndexSet, elementVector);
}
}
+ dofmatrix->finish();
}
@@ -326,40 +332,61 @@ namespace AMDiS
bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
ColView const& colLocalView,
ElementMatrix& elementMatrix,
- std::list<std::shared_ptr<OperatorType>>& operators)
+ std::list<std::shared_ptr<OperatorType>>& operators,
+ std::list<double*> const& factors)
{
const auto nRows = rowLocalView.tree().finiteElement().size();
const auto nCols = colLocalView.tree().finiteElement().size();
+ AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
+ "Nr. of operators and factors must be consistent!");
+
+ // fills the entire matrix with zeroes
elementMatrix.setSize(nRows, nCols);
- elementMatrix = 0.0; // fills the entire matrix with zeroes
+ elementMatrix = 0.0;
- for (auto op : operators)
- op->getElementMatrix(rowLocalView, colLocalView, elementMatrix);
+ bool add = false;
+
+ auto op_it = operators.begin();
+ auto fac_it = factors.begin();
+ for (; op_it != operators.end(); ++op_it, ++fac_it) {
+ bool add_op = (*op_it)->getElementMatrix(rowLocalView, colLocalView, elementMatrix, *fac_it);
+ add = add || add_op;
+ }
- return !operators.empty();
+ return add;
}
-
template <class Traits>
template <class RowView>
bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
ElementVector& elementVector,
- std::list<std::shared_ptr<OperatorType>>& operators)
+ std::list<std::shared_ptr<OperatorType>>& operators,
+ std::list<double*> const& factors)
{
const auto nRows = rowLocalView.tree().finiteElement().size();
+ AMDIS_TEST_EXIT_DBG(operators.size() == factors.size(),
+ "Nr. of operators and factors must be consistent!");
+
// Set all entries to zero
elementVector.resize(nRows);
elementVector = 0.0;
-
- for (auto op : operators)
- op->getElementVector(rowLocalView, elementVector);
- return !operators.empty();
+ bool add = false;
+
+ auto op_it = operators.begin();
+ auto fac_it = factors.begin();
+ for (; op_it != operators.end(); ++op_it, ++fac_it) {
+ bool add_op = (*op_it)->getElementVector(rowLocalView, elementVector, *fac_it);
+ add = add || add_op;
+ }
+
+ return add;
}
+
template <class Traits>
template <class Matrix, class RowIndexSet, class ColIndexSet>
void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
@@ -378,6 +405,7 @@ namespace AMDiS
}
}
+
template <class Traits>
template <class Vector, class IndexSet>
void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
diff --git a/dune/amdis/Traits.hpp b/dune/amdis/Traits.hpp
index 5418b2d9..30327fc0 100644
--- a/dune/amdis/Traits.hpp
+++ b/dune/amdis/Traits.hpp
@@ -3,67 +3,67 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
-namespace AMDiS {
-
-// some tags representing Traits classes or categories
-struct _none {};
-struct _scalar {};
-struct _vector {};
-struct _matrix {};
+namespace AMDiS
+{
+ // some tags representing Traits classes or categories
+ struct _none {};
+ struct _scalar {};
+ struct _vector {};
+ struct _matrix {};
-template <size_t I>
-struct VectorComponent
-{
+ template <size_t I>
+ struct VectorComponent
+ {
static constexpr size_t index = I;
-};
+ };
-template <size_t R, size_t C>
-struct MatrixComponent
-{
- static constexpr size_t row = R;
- static constexpr size_t col = C;
-};
+ template <size_t R, size_t C>
+ struct MatrixComponent
+ {
+ static constexpr size_t row = R;
+ static constexpr size_t col = C;
+ };
-template <class T>
-struct ValueCategory;
+ template <class T>
+ struct ValueCategory;
-template <class T>
-using ValueCategory_t = typename ValueCategory<T>::type;
+ template <class T>
+ using ValueCategory_t = typename ValueCategory<T>::type;
-template <class K, int SIZE>
-struct ValueCategory< Dune::FieldVector<K, SIZE> >
-{
+ template <class K, int SIZE>
+ struct ValueCategory< Dune::FieldVector<K, SIZE> >
+ {
using type = _vector;
-};
+ };
-template <class K>
-struct ValueCategory< Dune::FieldVector<K, 1> >
-{
+ template <class K>
+ struct ValueCategory< Dune::FieldVector<K, 1> >
+ {
using type = _scalar;
-};
+ };
-template <class K, int ROWS, int COLS>
-struct ValueCategory< Dune::FieldMatrix<K, ROWS, COLS> >
-{
+ template <class K, int ROWS, int COLS>
+ struct ValueCategory< Dune::FieldMatrix<K, ROWS, COLS> >
+ {
using type = _matrix;
-};
+ };
-template <class K>
-struct ValueCategory< Dune::FieldMatrix<K, 1, 1> >
-{
+ template <class K>
+ struct ValueCategory< Dune::FieldMatrix<K, 1, 1> >
+ {
using type = _scalar;
-};
+ };
-template <> struct ValueCategory<float> { using type = _scalar; };
-template <> struct ValueCategory<double> { using type = _scalar; };
-template <> struct ValueCategory<long double> { using type = _scalar; };
-template <> struct ValueCategory<int> { using type = _scalar; };
-template <> struct ValueCategory<long> { using type = _scalar; };
-template <> struct ValueCategory<long long> { using type = _scalar; };
-template <> struct ValueCategory<unsigned int> { using type = _scalar; };
-template <> struct ValueCategory<unsigned long> { using type = _scalar; };
-template <> struct ValueCategory<unsigned long long> { using type = _scalar; };
+ template <> struct ValueCategory<float> { using type = _scalar; };
+ template <> struct ValueCategory<double> { using type = _scalar; };
+ template <> struct ValueCategory<long double> { using type = _scalar; };
+ template <> struct ValueCategory<int> { using type = _scalar; };
+ template <> struct ValueCategory<long> { using type = _scalar; };
+ template <> struct ValueCategory<long long> { using type = _scalar; };
+ template <> struct ValueCategory<unsigned int> { using type = _scalar; };
+ template <> struct ValueCategory<unsigned long> { using type = _scalar; };
+ template <> struct ValueCategory<unsigned long long> { using type = _scalar; };
} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/LinearSolverInterface.hpp b/dune/amdis/linear_algebra/LinearSolverInterface.hpp
index 0f9f7192..fb188340 100644
--- a/dune/amdis/linear_algebra/LinearSolverInterface.hpp
+++ b/dune/amdis/linear_algebra/LinearSolverInterface.hpp
@@ -48,7 +48,7 @@ namespace AMDiS
}
// return the runner/worker corresponding to this solver (optional)
- virtual std::shared_ptr<RunnerBase> getRunner() {};
+ virtual std::shared_ptr<RunnerBase> getRunner() { /* Does nothing */ };
private:
/// main methods that all solvers must implement
@@ -61,4 +61,3 @@ namespace AMDiS
struct SolverCreator;
} // end namespace AMDiS
-
diff --git a/dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp b/dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
index 6acb2e0b..5fe88c34 100644
--- a/dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
+++ b/dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
@@ -57,12 +57,12 @@ namespace AMDiS
std::array<mtl::irange, _M> r_cols;
getRanges(r_rows, r_cols);
- For<0, _N>::loop([this, &r_rows, &r_cols, &b, &x](const auto _i) {
+ For<0, _N>::loop([&](const auto _i) {
bool first = true;
// a reference to the i'th block of x
VectorOut x_i(x[r_rows[_i]]);
- For<0, _M>::loop([this, &b, &x_i, &r_cols, &first, _i](const auto _j) {
+ For<0, _M>::loop([&](const auto _j) {
auto const& A_ij = this->operator()(_i, _j);
if (num_rows(A_ij) > 0) {
// a reference to the j'th block of b
@@ -114,7 +114,8 @@ namespace AMDiS
/// Fill two arrays of irange corresponding to row and column sizes.
/// \see getRowRanges() and \see getColRanges()
- void getRanges(std::array<mtl::irange, _N>& r_rows, std::array<mtl::irange, _M>& r_cols) const
+ void getRanges(std::array<mtl::irange, _N>& r_rows,
+ std::array<mtl::irange, _M>& r_cols) const
{
getRowRanges(r_rows);
getColRanges(r_cols);
diff --git a/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp b/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
index b39e55f0..15d3149a 100644
--- a/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp
@@ -62,18 +62,19 @@ namespace AMDiS
inline void set_to_zero(BlockMTLVector<MTLVector, _N>& vec)
{
For<0, _N>::loop([&](const auto _i) {
- set_to_zero(std::get<_i>(vec));
+ set_to_zero(std::get<_i>(vec));
});
}
/// A wrapper, that creates a contiguos vector corresponding to a block-vector
- /// and copy the value on construction ans possibly back on destruction.
+ /// and copy the value on construction and eventually back on destruction, if
+ /// required.
template <class BlockVector, class Vector>
struct BlockVectorWrapper
{
BlockVectorWrapper(BlockVector& blockVector,
- bool copyBack = std::is_const<BlockVector>::value)
+ bool copyBack = !std::is_const<BlockVector>::value)
: blockVector(blockVector)
, vector(num_rows(blockVector))
, copyBack(copyBack)
@@ -109,7 +110,7 @@ namespace AMDiS
/// do not copy vector to block-vector since this is const
template <bool Copy>
- std::enable_if_t<!Copy> assignFrom(bool_<Copy>) {}
+ std::enable_if_t<!Copy> assignFrom(bool_<Copy>) { /* Do nothing */ }
/// copy from vector to block-vector
template <bool Copy>
@@ -128,7 +129,7 @@ namespace AMDiS
BlockVector& blockVector;
Vector vector;
- /// Copy data back to block-vector of destruction
+ /// Copy data back to block-vector on destruction
bool copyBack;
};
diff --git a/dune/amdis/linear_algebra/mtl/DOFVector.hpp b/dune/amdis/linear_algebra/mtl/DOFVector.hpp
index 6d6f89ec..00a6c42c 100644
--- a/dune/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -92,16 +92,16 @@ namespace AMDiS
/// Access the entry \p i of the \ref vector with read-access.
value_type const& operator[](size_type i) const
{
- AMDIS_TEST_EXIT_DBG( i < vector->size() ,
- "Index " << i << " out of range [0," << vector->size() << ")" );
+ AMDIS_TEST_EXIT_DBG( i < num_rows(*vector) ,
+ "Index " << i << " out of range [0," << num_rows(*vector) << ")" );
return (*vector)[i];
}
/// Access the entry \p i of the \ref vector with write-access.
value_type& operator[](size_type i)
{
- AMDIS_TEST_EXIT_DBG( i < vector->size() ,
- "Index " << i << " out of range [0," << vector->size() << ")" );
+ AMDIS_TEST_EXIT_DBG( i < num_rows(*vector) ,
+ "Index " << i << " out of range [0," << num_rows(*vector) << ")" );
return (*vector)[i];
}
@@ -128,11 +128,20 @@ namespace AMDiS
/// The name of the DOFVector (used in file-writing)
std::string name;
- /// The data-vector (can hold a new BAseVector or a pointer to external data
+ /// The data-vector (can hold a new BaseVector or a pointer to external data
ClonablePtr<BaseVector> vector;
// friend class declarations
template <class, class> friend class SystemVector;
};
+
+
+ /// Constructor a dofvector from given feSpace and name
+ template <class FeSpaceType, class ValueType = double>
+ DOFVector<FeSpaceType, ValueType>
+ make_dofvector(FeSpaceType const& feSpace, std::string name)
+ {
+ return {feSpace, name};
+ }
} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp b/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp
index fab31a65..99f5dae7 100644
--- a/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp
+++ b/dune/amdis/linear_algebra/mtl/KrylovRunner.hpp
@@ -69,7 +69,6 @@ namespace AMDiS
L->init(A);
R->init(A);
}
-
/// Implementation of \ref RunnerInterface::exit()
virtual void exit() override
@@ -78,7 +77,6 @@ namespace AMDiS
R->exit();
}
-
/// Implementation of \ref RunnerInterface::solve()
virtual int solve(Matrix const& A, Vector& x, Vector const& b,
SolverInfo& solverInfo) override
diff --git a/dune/amdis/linear_algebra/mtl/SystemVector.hpp b/dune/amdis/linear_algebra/mtl/SystemVector.hpp
index 7bae7721..604b5443 100644
--- a/dune/amdis/linear_algebra/mtl/SystemVector.hpp
+++ b/dune/amdis/linear_algebra/mtl/SystemVector.hpp
@@ -49,11 +49,11 @@ namespace AMDiS
std::vector<std::string> const& names,
MultiVector&& multiVector)
{
- return BuildDOFVectors<DOFVectors>::make(
- MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
- std::forward<FeSpaces>(feSpaces),
- names,
- std::forward<MultiVector>(multiVector));
+ return BuildDOFVectors<DOFVectors>::make(
+ MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
+ std::forward<FeSpaces>(feSpaces),
+ names,
+ std::forward<MultiVector>(multiVector));
}
} // end namespace Impl
@@ -232,6 +232,27 @@ namespace AMDiS
/// a tuple of dofvectors referencing \ref vector
DOFVectors dofvectors;
};
+
+
+ /// Construct a systemvector from given feSpace-tuple and names-vector
+ template <class FeSpacesType, class ValueType = double>
+ SystemVector<FeSpacesType, ValueType>
+ make_systemvector(FeSpacesType const& feSpaces, std::vector<std::string> names)
+ {
+ return {feSpaces, names};
+ }
+
+ /// Construct a systemvector from given feSpace-tuple. Create vector of names
+ /// from base_name, i.e. {base_name_0, base_name_1, ...}
+ template <class FeSpaces, class ValueType = double>
+ SystemVector<FeSpaces, ValueType>
+ make_systemvector(FeSpaces const& feSpaces, std::string base_name)
+ {
+ std::vector<std::string> names;
+ for (size_t i = 0; i < std::tuple_size<FeSpaces>::value; ++i)
+ names.push_back(base_name + "_" + std::to_string(i));
+ return {feSpaces, names};
+ }
#ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
diff --git a/dune/amdis/test/CMakeLists.txt b/dune/amdis/test/CMakeLists.txt
new file mode 100644
index 00000000..46b9a735
--- /dev/null
+++ b/dune/amdis/test/CMakeLists.txt
@@ -0,0 +1,18 @@
+# maybe unncessary (?)
+dune_add_test(NAME test_duneamdis
+ TARGET duneamdis
+ COMPILE_ONLY)
+
+find_package(GTest REQUIRED)
+
+set(projects "test1")
+foreach(project ${projects})
+ dune_add_test(NAME ${project}
+ SOURCES ${project}.cc
+ LINK_LIBRARIES duneamdis ${GTEST_BOTH_LIBRARIES}
+ COMPILE_DEFINITIONS "DIM=2;DOW=2"
+ CMD_ARGS "init/test.json.2d")
+ add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project})
+ target_include_directories(${project} PRIVATE ${GTEST_INCLUDE_DIRS})
+ set_tests_properties(${project} PROPERTIES WORKING_DIRECTORY ${CMAKE_SOURCE_DIR})
+endforeach()
\ No newline at end of file
diff --git a/dune/amdis/test/kdtree.hpp b/dune/amdis/test/kdtree.hpp
new file mode 100644
index 00000000..fd22dca2
--- /dev/null
+++ b/dune/amdis/test/kdtree.hpp
@@ -0,0 +1,245 @@
+#pragma once
+
+#include <nanoflann.hpp>
+
+#include <cstdlib>
+#include <iostream>
+
+
+namespace AMDiS { namespace extensions {
+
+using namespace nanoflann;
+
+// =====================================================================================
+typedef WorldVector<double> PointType;
+typedef boost::tuple<const MacroElement*, int, unsigned long> DataType;
+typedef std::vector<PointType> VectorOfPointsType;
+typedef std::vector<DataType> VectorOfDataType;
+// =====================================================================================
+
+
+ /** A simple vector-of-vectors adaptor for nanoflann, without duplicating the storage.
+ * The i'th vector represents a point in the state space.
+ *
+ * \tparam dim If set to >0, it specifies a compile-time fixed dimensionality for the points in the data set, allowing more compiler optimizations.
+ * \tparam value_type The type of the point coordinates (typically, double or float).
+ * \tparam Distance The distance metric to use: nanoflann::metric_L1, nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc.
+ * \tparam index_type The type for indices in the KD-tree index (typically, size_t of int)
+ */
+ template < class FeSpace >
+ struct KDTreeMeshAdaptor
+ {
+ using MeshView = typename FeSpace::GridView;
+ using SizeType = typename FeSpace::LocalIndexSet::size_type;
+ using ValueType = double;
+
+ static constexpr dimension = MeshView::dimension;
+ static constexpr dimensionworld = MeshView::dimensionworld;
+
+ using Codim0 = typename GridView::template Codim<0>;
+ using PointType = typename Codim0::Geometry::GlobalCoordinate;
+ using VectorOfPointsType = std::vector<PointType> ;
+
+ using DataType = std::vector<Dune::FieldVector<double,1> >;
+ using VectorOfDataType = std::vector<std::pair<std::hash<DataType>, DataType>>>;
+
+ using Distance = nanoflann::metric_L2_Simple;
+
+ using Self = KDTreeMeshAdaptor<VectorOfPointsType, VectorOfDataType, ValueType, dimensionworld, Distance, SizeType>;
+ using metric_t = typename Distance::template traits<ValueType, Self>::distance_t;
+ using index_t = KDTreeSingleIndexAdaptor<metric_t, Self, dimensionworld, SizeType>;
+
+ public:
+ std::unique_ptr<index_t> index; //! The kd-tree index for the user to call its methods as usual with any other FLANN index.
+ FeSpace const& feSpace;
+
+ VectorOfPointsType_ m_points;
+ VectorOfDataType_ m_data;
+
+ std::map<std::hash<DataType>, DataType> shapeValues;
+
+ public:
+ /// Constructor: takes a const ref to the vector of vectors object with the data points
+ KDTreeMeshAdaptor(FeSpace const& feSpace_,
+ const int leaf_max_size = 10)
+ : feSpace(feSpace_)
+ {
+ init();
+
+ index = std::make_unique<index_t>(dimensionworld, *this /* adaptor */,
+ nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
+ index->buildIndex();
+ }
+
+ /// Destructor.
+ ~KDTreeMeshAdaptor() {}
+
+ void init()
+ {
+ auto const& localView = feSpace.localView();
+ auto const& indexSet = feSpace.localIndexSet();
+
+ for (auto const& element : elements(feSpace.gridView())) {
+ localView.bind(element);
+ localIndexSet.bind(localView);
+
+ m_points.push_back(element.geometry().center());
+
+ auto const& localBasis = localView.tree().finiteElement().localBasis();
+ std::vector<SizeType> indices(localBasis.size());
+ for (size_t j = 0; j < localBasis.size(); ++j)
+ indices[i] = localIndexSet.index(j);
+
+ m_data.push_back(indices);
+ }
+ }
+
+ void reinit(const int leaf_max_size = 10)
+ {
+ m_points.clear();
+ m_data.clear();
+ init();
+ index.reset(new index_t(dimensionworld, *this /* adaptor */,
+ nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
+ index->buildIndex();
+ }
+
+ /** Query for the \a num_closest closest points to a given point (entered as query_point[0:dim-1]).
+ * Note that this is a short-cut method for index->findNeighbors().
+ * The user can also call index->... methods as desired.
+ * \note nChecks_IGNORED is ignored but kept for compatibility with the original FLANN interface.
+ */
+ inline void query(const ValueType *query_point,
+ const size_t num_closest,
+ SizeType *out_indices,
+ ValueType *out_distances_sq,
+ const int nChecks_IGNORED = 10) const
+ {
+ nanoflann::KNNResultSet<ValueType, SizeType> resultSet(num_closest);
+ resultSet.init(out_indices, out_distances_sq);
+ index->findNeighbors(resultSet, query_point, nanoflann::SearchParams());
+ }
+
+ /** @name Interface expected by KDTreeSingleIndexAdaptor
+ * @{ */
+
+ Self const& derived() const
+ {
+ return *this;
+ }
+ Self& derived()
+ {
+ return *this;
+ }
+
+ // Must return the number of data points
+ inline size_t kdtree_get_point_count() const
+ {
+ return m_points.size();
+ }
+
+ // Returns the distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class:
+ inline value_type kdtree_distance(const double *p1,
+ const SizeType idx_p2,
+ size_t size) const
+ {
+ double s = 0.0;
+ for (size_t i = 0; i < size; i++) {
+ const double d = p1[i] - m_points[idx_p2][i];
+ s += d*d;
+ }
+ return s;
+ }
+
+ // Returns the dim'th component of the idx'th point in the class:
+ inline value_type kdtree_get_pt(const SizeType idx, size_t dim) const
+ {
+ return m_points[idx][dim];
+ }
+
+ // Optional bounding-box computation: return false to default to a standard bbox computation loop.
+ // Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it again.
+ // Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds)
+ template <class BBOX>
+ bool kdtree_get_bbox(BBOX &bb) const
+ {
+ return false;
+ }
+
+ /** @} */
+
+
+ /**
+ * strategy and nrOfPoints are ignored
+ * == evalAtPoint_simple
+ **/
+ template <class DOFVector>
+ typename DOFVector::value_type
+ evalAtPoint(DOFVector const& vec, PointType const& x, int strategy = 0, int nrOfPoints = -1)
+ {
+ using T = typename DOFVector::value_type;
+ T value = 0;
+
+ std::vector<SizeType> indices;
+ if (getNearestElement(x, indices)) {
+ // get lambda of point...
+
+ return localEval(vec, indices, lambda);
+ }
+ return value;
+ }
+
+ private:
+ /// evaluate DOFVector at barycentric coordinates \p lambda in element that
+ /// is bind to instance in \ref bind().
+ template <class FE, class ValueType, class Indices, class LocalCoord>
+ ValueType localEval(DOFVector<FE, ValueType> const& vec,
+ Indices const& indices, LocalCoord const& lambda)
+ {
+ auto const& localBasis = localView.tree().finiteElement().localBasis();
+
+ // find a way to evaluate basisFcts at lambda!...
+ auto const& shape = shapeValues[indices.first];
+ localBasis.evaluateFunction(lambda, shapeValues);
+
+ assert( shape.size() == indices.size() );
+
+ ValueType data = 0;
+ for (size_t j = 0; j < shape.size(); ++j)
+ data += vec[indices.second[j]] * shape[j];
+
+ return data;
+ }
+
+ template <class Indices>
+ bool getNearestElement(PointType& x, Indices& indices)
+ {
+ const size_t nnp = 1;
+ std::vector<DegreeOfFreedom> ret_indexes(nnp);
+ std::vector<double> out_dists_sqr(nnp);
+ nanoflann::KNNResultSet<double, DegreeOfFreedom> resultSet(nnp);
+
+ resultSet.init(&ret_indexes[0], &out_dists_sqr[0]);
+ index->findNeighbors(resultSet, x.begin(), nanoflann::SearchParams());
+
+ indices = m_data[ret_indexes[0]];
+ return true;
+ }
+
+ }; // end of KDTreeMeshAdaptor
+
+
+
+ template <class FeSpace>
+ auto provideKDTree(FeSpace const& gridView)
+ {
+ using GridView = typename FeSpace::GridView;
+ static constexpr int dow = GridView::dimensionworld;
+
+ using KD_Tree = KDTreeMeshAdaptor<FeSpace, VectorOfPointsType, VectorOfDataType, double, dow>;
+
+ return std::make_unique<KD_Tree_M>(dow, gridView);
+ }
+
+
+} }
diff --git a/dune/amdis/test/macro.stand.2d b/dune/amdis/test/macro.stand.2d
new file mode 100644
index 00000000..7b0c2db0
--- /dev/null
+++ b/dune/amdis/test/macro.stand.2d
@@ -0,0 +1,30 @@
+DIM: 2
+DIM_OF_WORLD: 2
+
+number of elements: 4
+number of vertices: 5
+
+element vertices:
+0 1 4
+1 2 4
+2 3 4
+3 0 4
+
+element boundaries:
+0 0 1
+0 0 2
+0 0 3
+0 0 4
+
+vertex coordinates:
+ 0.0 0.0
+ 1.0 0.0
+ 1.0 1.0
+ 0.0 1.0
+ 0.5 0.5
+
+element neighbours:
+1 3 -1
+2 0 -1
+3 1 -1
+0 2 -1
diff --git a/dune/amdis/test/test.json.2d b/dune/amdis/test/test.json.2d
new file mode 100644
index 00000000..7242fe2b
--- /dev/null
+++ b/dune/amdis/test/test.json.2d
@@ -0,0 +1,29 @@
+{
+ "dimension of world": 2,
+
+ "elliptMesh": {
+ "macro file name": "macro.stand.2d",
+ "global refinements": 10
+ },
+
+ "ellipt": {
+ "mesh": "elliptMesh",
+ "dim": 2,
+ "components": 1,
+ "polynomial degree[0]": 1,
+
+ "solver": "cg",
+ "solver" : {
+ "max iteration": 1000,
+ "tolerance": 1e-8,
+ "info": 10,
+ "left precon": "diag"
+ },
+
+ "output": {
+ "filename": "ellipt.2d",
+ "ParaView format": 1,
+ "ParaView mode": 1
+ }
+ }
+}
diff --git a/dune/amdis/test/test1.cc b/dune/amdis/test/test1.cc
new file mode 100644
index 00000000..62f13ee2
--- /dev/null
+++ b/dune/amdis/test/test1.cc
@@ -0,0 +1,150 @@
+#include <dune/amdis/AMDiS.hpp>
+#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Literals.hpp>
+
+#include <gtest/gtest.h>
+
+#ifndef DIM
+#define DIM 2
+#endif
+#ifndef DOW
+#define DOW 2
+#endif
+
+using namespace AMDiS;
+namespace {
+
+ using TestParam = ProblemStatTraits<DIM, DOW, 1>;
+ using TestProblem = ProblemStat<TestParam>;
+
+ using MeshView = typename TestProblem::MeshView;
+ using Geometry = typename MeshView::template Codim<0>::Geometry;
+ using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
+
+ class AMDiSTester : public ::testing::Test
+ {
+ protected:
+ // create a problemStat
+ virtual void SetUp() override
+ {
+ prob = std::make_shared<TestProblem>("test");
+ prob->initialize(INIT_ALL);
+ }
+
+ std::shared_ptr<TestProblem> prob;
+ };
+
+
+ template <class FeSpace>
+ class LocalEvaluation
+ {
+ public:
+ LocalEvaluation(FeSpace const& feSpace)
+ : localView(feSpace.localView())
+ , localIndexSet(feSpace.localIndexSet())
+ {}
+
+ template <class Element>
+ void bind(Element const& element)
+ {
+ localView.bind(element);
+ localIndexSet.bind(localView);
+ }
+
+ /// evaluate DOFVector at barycentric coordinates \p lambda in element that
+ /// is bind to instance in \ref bind().
+ template <class ValueType, class LocalCoord>
+ ValueType eval(DOFVector<FeSpace, ValueType> const& vec, LocalCoord const& lambda)
+ {
+ auto const& localBasis = localView.tree().finiteElement().localBasis();
+ localBasis.evaluateFunction(lambda, shapeValues);
+
+ ValueType data = 0;
+ for (size_t j = 0; j < shapeValues.size(); ++j) {
+ const auto global_idx = localIndexSet.index(j);
+ data += vec[global_idx] * shapeValues[j];
+ }
+
+ return data;
+ }
+
+ private:
+ typename FeSpace::LocalView localView;
+ typename FeSpace::LocalIndexSet localIndexSet;
+
+ std::vector<Dune::FieldVector<double,1> > shapeValues;
+ };
+
+
+ TEST_F(AMDiSTester, CreateDOFVector)
+ {
+ using FeSpace = typename TestProblem::template FeSpace<0>;
+
+ // construction of DOFVectors
+ DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
+ auto vec2 = make_dofvector(prob->getFeSpace(0_c), "vec2");
+
+ // retrieving DOFVectors from problemStat
+ auto vec3 = prob->getSolution(0_c);
+ auto vec5 = prob->getSolution<0>();
+ auto vec6 = prob->getSolution(index_<0>());
+
+ auto vec7 = prob->getSolution()->getDOFVector(0_c);
+ auto vec8 = prob->getSolution()->getDOFVector<0>();
+ auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
+
+ auto& sys_vec = *prob->getSolution();
+ auto vec10 = sys_vec[0_c];
+ auto vec11 = sys_vec[index_<0>()];
+
+ // construction of SystemVector
+ using FeSpaces = typename TestProblem::FeSpaces;
+ SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
+ auto sys_vec2 = make_systemvector(*prob->getFeSpaces(), prob->getComponentNames());
+ auto sys_vec3 = make_systemvector(*prob->getFeSpaces(), "sys_vec");
+
+ // retrieving systemVector from problemStat
+ auto sys_vec4 = *prob->getSolution();
+ }
+
+
+ TEST_F(AMDiSTester, FillDOFVector)
+ {
+ auto& vec = prob->getSolution(0_c);
+ auto const& feSpace = vec.getFeSpace();
+
+ // interpolate function to dofvector
+ vec.interpol([](auto const& x) { return x*x; });
+
+ // test whether interpolation is fine
+ using FeSpace = std::decay_t< decltype(feSpace) >;
+ LocalEvaluation<FeSpace> evaluator(feSpace);
+
+ for (auto const& element : elements(feSpace.gridView())) {
+ evaluator.bind(element);
+
+ auto const& refElement = RefElements::general(element.type());
+
+ size_t nVertices = element.subEntities(DIM);
+ for (size_t i = 0; i < nVertices; ++i) {
+ auto pos = refElement.position(i, DIM);
+ auto data = evaluator.eval(vec, pos);
+
+ auto x = element.geometry().global(pos);
+ EXPECT_EQ(data, x*x);
+ }
+ }
+ }
+
+} // end namespace
+
+int main(int argc, char** argv)
+{
+ AMDiS::init(argc, argv);
+
+ ::testing::InitGoogleTest(&argc, argv);
+ int error_code = RUN_ALL_TESTS();
+
+ AMDiS::finalize();
+ return error_code;
+}
\ No newline at end of file
diff --git a/init/heat.json.2d b/init/heat.json.2d
index f5607ad3..9af04026 100644
--- a/init/heat.json.2d
+++ b/init/heat.json.2d
@@ -3,7 +3,7 @@
"heatMesh": {
"macro file name": "./macro/macro.stand.2d",
- "global refinements": 15,
+ "global refinements": 4,
"min corner": "0,0",
"max corner": "1,1",
"num cells": "2,2",
@@ -35,6 +35,6 @@
"adapt": {
"timestep": 0.1,
"start time": 0.0,
- "end time": 0.1
+ "end time": 1.0
}
}
diff --git a/init/test.json.2d b/init/test.json.2d
new file mode 100644
index 00000000..8eb8d951
--- /dev/null
+++ b/init/test.json.2d
@@ -0,0 +1,29 @@
+{
+ "dimension of world": 2,
+
+ "testMesh": {
+ "macro file name": "macro/macro.stand.2d",
+ "global refinements": 10
+ },
+
+ "test": {
+ "mesh": "testMesh",
+ "dim": 2,
+ "components": 1,
+ "polynomial degree[0]": 1,
+
+ "solver": "cg",
+ "solver" : {
+ "max iteration": 1000,
+ "tolerance": 1e-8,
+ "info": 10,
+ "left precon": "diag"
+ },
+
+ "output": {
+ "filename": "output/test.2d",
+ "ParaView format": 1,
+ "ParaView mode": 1
+ }
+ }
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 40d2807a..5b26fd10 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,5 @@
-set(projects "heat" "pfc" "stokes" "navier_stokes")
+set(projects "heat" "heat_old" "pfc" "stokes" "navier_stokes")
foreach(project ${projects})
add_executable(${project} ${project}.cc)
diff --git a/src/heat.cc b/src/heat.cc
index 44bd162b..ba29cef3 100644
--- a/src/heat.cc
+++ b/src/heat.cc
@@ -22,67 +22,69 @@
using namespace AMDiS;
-// class HeatParam
-// {
-// template <class... Bs>
-// using FeSpaceTuple = std::tuple<Bs...>;
-//
-// template <class Mesh, int deg>
-// using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
-//
-// public:
-// static constexpr int dim = DIM;
-// static constexpr int dimworld = DOW;
-// static constexpr int nComponents = 1;
-//
-// // default values
-// using Mesh = Dune::YaspGrid<dim>;
-// using FeSpaces = FeSpaceTuple<Lagrange<Mesh, 2>>;
-// };
+class HeatParam
+{
+ template <class... Bs>
+ using FeSpaceTuple = std::tuple<Bs...>;
+
+ template <class Mesh, int deg>
+ using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
+
+public:
+ static constexpr int dim = DIM;
+ static constexpr int dimworld = DOW;
+ static constexpr int nComponents = 1;
+
+ // default values
+ using Mesh = Dune::YaspGrid<dim>;
+ using FeSpaces = FeSpaceTuple<Lagrange<Mesh, 2>>;
+};
// 1 component with polynomial degree 1
-using HeatParam = ProblemStatTraits<DIM, DOW, 2>;
+// using HeatParam = ProblemStatTraits<DIM, DOW, 1>;
using HeatProblem = ProblemStat<HeatParam>;
using HeatProblemInstat = ProblemInstat<HeatParam>;
int main(int argc, char** argv)
{
- Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
- AMDiS::init(argc, argv);
-
- HeatProblem prob("heat");
- prob.initialize(INIT_ALL);
-
- HeatProblemInstat probInstat("heat", prob);
- probInstat.initialize(INIT_UH_OLD);
-
- AdaptInfo adaptInfo("adapt");
-
- double tau = adaptInfo.getTimestep();
-
- using Op = HeatProblem::OperatorType;
- Op opLhs;
- opLhs.addZOT( constant(1.0/tau) );
- opLhs.addSOT( constant(1.0) );
-
- Op opRhs;
- opRhs.addZOT( valueOf(prob.getSolution(0_c), 1.0/tau) );
- opRhs.addZOT( eval([](auto const& x) { return -1.0; }) );
-
- prob.addMatrixOperator(opLhs, 0, 0);
- prob.addVectorOperator(opRhs, 0);
-
- // set boundary condition
- auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
- auto dbcValues = [](auto const& p){ return 0.0; };
- prob.addDirichletBC(predicate, 0, 0, dbcValues);
-
- *prob.getSolution() = 0.0;
-
- AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
- adapt.adapt();
-
- AMDiS::finalize();
- return 0;
+ Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+ AMDiS::init(argc, argv);
+
+ HeatProblem prob("heat");
+ prob.initialize(INIT_ALL);
+
+ HeatProblemInstat probInstat("heat", prob);
+ probInstat.initialize(INIT_UH_OLD);
+
+ AdaptInfo adaptInfo("adapt");
+
+ using Op = HeatProblem::OperatorType;
+ Op opTimeLhs, opL, opTimeRhs, opForce;
+
+ opTimeLhs.addZOT( constant(1.0) );
+ prob.addMatrixOperator(opTimeLhs, 0, 0, probInstat.getInvTau());
+
+ opL.addSOT( constant(1.0) );
+ prob.addMatrixOperator(opL, 0, 0);
+
+ opTimeRhs.addZOT( valueOf(prob.getSolution(0_c)) );
+ prob.addVectorOperator(opTimeRhs, 0, probInstat.getInvTau());
+
+ opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
+ prob.addVectorOperator(opForce, 0);
+
+
+ // set boundary condition
+ auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
+ auto dbcValues = [](auto const& p){ return 0.0; };
+ prob.addDirichletBC(predicate, 0, 0, dbcValues);
+
+ *prob.getSolution() = 0.0;
+
+ AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
+ adapt.adapt();
+
+ AMDiS::finalize();
+ return 0;
}
\ No newline at end of file
diff --git a/src/navier_stokes.cc b/src/navier_stokes.cc
index 414e2124..7669cb23 100644
--- a/src/navier_stokes.cc
+++ b/src/navier_stokes.cc
@@ -98,8 +98,10 @@ int main(int argc, char** argv)
t < adaptInfo.getEndTime();
t+= adaptInfo.getTimestep())
{
+ adaptInfo.setTime(t);
+
AMDIS_MSG("Timestep t = " << t);
- prob.writeFiles(adaptInfo, t);
+ prob.writeFiles(adaptInfo);
// assemble and solve system
prob.buildAfterCoarsen(adaptInfo, Flag(0));
@@ -107,7 +109,7 @@ int main(int argc, char** argv)
}
// output solution
- prob.writeFiles(adaptInfo, t);
+ prob.writeFiles(adaptInfo);
AMDiS::finalize();
return 0;
diff --git a/src/pfc.cc b/src/pfc.cc
index 46582dd7..33bea8aa 100644
--- a/src/pfc.cc
+++ b/src/pfc.cc
@@ -8,8 +8,10 @@
#include <ctime>
#include <cmath>
+#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/ProblemInstat.hpp>
#ifndef DIM
#define DIM 2
@@ -25,15 +27,18 @@ using namespace AMDiS;
// 3 component with polynomial degree 1
using PfcParam = ProblemStatTraits<DIM, DOW, 1, 1, 1>;
using PfcProblem = ProblemStat<PfcParam>;
+using PfcProblemInstat = ProblemInstat<PfcParam>;
-
-int main(int argc, char** argv)// try
+int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
PfcProblem prob("pfc");
prob.initialize(INIT_ALL);
+ PfcProblemInstat probInstat("pfc", prob);
+ probInstat.initialize();
+
AdaptInfo adaptInfo("adapt");
double tau = adaptInfo.getTimestep();
@@ -74,6 +79,7 @@ int main(int argc, char** argv)// try
prob.addMatrixOperator(opLhs11, 1, 1);
prob.addMatrixOperator(opLhs21, 2, 1);
prob.addMatrixOperator(opLhs22, 2, 2);
+
prob.addVectorOperator(opRhs0, 0);
prob.addVectorOperator(opRhs1, 1);
@@ -82,31 +88,14 @@ int main(int argc, char** argv)// try
Nu.getVector() = 0.0;
std::srand( time(0) );
for (size_t i = 0; i < Psi.getSize(); ++i)
- Psi[i] = (std::rand()/double(RAND_MAX) - 0.5) + psi_mean;
+ Psi[i] = (std::rand()/double(RAND_MAX) - 0.5) + psi_mean;
// using PreconType = PfcPrecon<typename PfcProblem::SystemMatrixType::MultiMatrix, typename PfcProblem::SystemVectorType::MultiVector>;
// PreconType precon(prob.getSystemMatrix().getMatrix(), &tau, M0);
- double t = 0.0;
- for (t = adaptInfo.getStartTime();
- t < adaptInfo.getEndTime();
- t+= adaptInfo.getTimestep())
- {
- prob.writeFiles(adaptInfo, t);
-
- AMDIS_MSG("Timestep t = " << t);
- AMDIS_MSG("----------------------------------");
- prob.buildAfterCoarsen(adaptInfo, Flag(0));
-// prob.solveAdvanced(adaptInfo, precon);
- }
-
- prob.writeFiles(adaptInfo, t);
+ AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
+ adapt.adapt();
+
AMDiS::finalize();
return 0;
}
-// catch (Dune::Exception &e){
-// std::cerr << "Dune reported error: " << e << std::endl;
-// }
-// catch (...){
-// std::cerr << "Unknown exception thrown!" << std::endl;
-// }
diff --git a/src/stokes.cc b/src/stokes.cc
index a1a569aa..8cefc5c6 100644
--- a/src/stokes.cc
+++ b/src/stokes.cc
@@ -34,7 +34,7 @@ int main(int argc, char** argv)
// define the differential operators
using Op = StokesProblem::OperatorType;
- For<0,DOW>::loop([&](auto const _i)
+ for_<0,DOW>([&](auto _i)
{
// <viscosity*grad(u_i), grad(v_i)>
Op* opL = new Op;
@@ -84,7 +84,7 @@ int main(int argc, char** argv)
prob.solve(adaptInfo);
// output solution
- prob.writeFiles(adaptInfo, 0.0);
+ prob.writeFiles(adaptInfo);
AMDiS::finalize();
return 0;
--
GitLab