adaptivity error found and corrected

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.0)
 project(dune-amdis CXX)
 
 set(ALBERTA_ROOT /opt/software/alberta)
-# set(ALBERTA_INCLUDE_DIR  /opt/software/alberta/include)
+set(ALBERTA_INCLUDE_DIR  /opt/software/alberta/include)
 set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
 set(UG_DIR /opt/software/dune/lib/cmake/ug)
 set(Vc_DIR /opt/software/dune/lib/cmake/Vc)

dune/amdis/Basic.hpp

@@ -27,37 +27,50 @@ namespace AMDiS
   
     namespace Impl
     {
-	// add arg to repeated constructor argument list
-	template <class Tuple, size_t N, class Arg, class... Args>
-	Tuple construct_tuple_aux(index_<N>, Arg&& arg, Args&&... args)
+	template <class Tuple, size_t N>
+	struct ConstructTuple
 	{
-	    return construct_tuple_aux<Tuple>(index_<N-1>(), 
-		std::forward<Arg>(arg), std::forward<Arg>(arg), std::forward<Args>(args)...);
-	}
+	    // 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, class... Args>
-	Tuple construct_tuple_aux(index_<1>, Args&&... args)
+	template <class Tuple>
+	struct ConstructTuple<Tuple, 1>
 	{
-	    static_assert(std::tuple_size<Tuple>::value == sizeof...(args), 
-		"Nr. of argument != tuple-size");
-	    return Tuple{std::forward<Args>(args)...};
-	}
+	    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, class... Args>
-	Tuple construct_tuple_aux(index_<0>, Args&&... args)
+	template <class Tuple>
+	struct ConstructTuple<Tuple, 0>
 	{
-	    static_assert(std::tuple_size<Tuple>::value == 0, 
-		"Construction of empty tuples with empty argument list only!");
-	    return {};
-	}
+	    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 {};
+	    }
+	};
+	
     } // 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::construct_tuple_aux<Tuple>(index_<std::tuple_size<Tuple>::value>(), 
-	std::forward<Arg>(arg));
+	return Impl::ConstructTuple<Tuple, std::tuple_size<Tuple>::value>::make( 
+	    std::forward<Arg>(arg));
     }
     
     // -----------
@@ -119,4 +132,10 @@ namespace AMDiS
 	typename RepeatedTuple<T, MakeSeq_t<N>>::type;
     
   
+    // -----------
+	
+    template <class T>
+    using owner = T;
+	
+	
 } // end namespace AMDiS

dune/amdis/ClonablePtr.hpp

+#pragma once
+
+#include "Basic.hpp"
+
+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;
+
+	/// 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
+    };
+    
+    template <class T> 
+    void swap(ClonablePtr<T>& a, ClonablePtr<T>& b) noexcept
+    {
+	a.swap(b);
+    }
+
+} // end namespace AMDiS

dune/amdis/DOFVector.hpp

@@ -6,6 +6,7 @@
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/istl/bvector.hh>
 
+#include "ClonablePtr.hpp"
 #include "Log.hpp"
 
 namespace AMDiS
@@ -22,30 +23,23 @@ namespace AMDiS
 	
 	using value_type = ValueType;
 	
-	/// Constructor.
+	/// Constructor. Constructs new BaseVector.
 	DOFVector(FeSpace const& feSpace, std::string name)
 	  : feSpace(feSpace)
 	  , name(name)
-	  , vector(new BaseVector())
-	  , allocated(true)
+	  , vector(ClonablePtr<BaseVector>::make())
 	{
 	    compress();
 	}
 	
 	/// Constructor. Takes pointer of data-vector.
 	DOFVector(FeSpace const& feSpace, std::string name, 
-		  BaseVector& vector)
+		  BaseVector& vector_ref)
 	  : feSpace(feSpace)
 	  , name(name)
-	  , vector(&vector)
+	  , vector(vector_ref)
 	{}
 	
-	~DOFVector()
-	{
-	    if (allocated)
-		delete vector;
-	}
-	
 	/// Return the basis \ref feSpace of the vector
 	FeSpace const& getFeSpace() const
 	{
@@ -79,7 +73,8 @@ namespace AMDiS
 	/// Resize the \ref vector to the size of the \ref feSpace.
 	void compress()
 	{
-	    vector->resize(getSize());
+	    if (vector->size() != getSize())
+		vector->resize(getSize());
 	}
 	
 	
@@ -111,8 +106,7 @@ namespace AMDiS
 	FeSpace const&	feSpace;
 	std::string	name;
 	
-	BaseVector* vector;
-	bool allocated = false;
+	ClonablePtr<BaseVector> vector;
 	
     // friend class declarations
     template <class, class>

dune/amdis/DirichletBC.inc.hpp

@@ -50,7 +50,7 @@ namespace AMDiS
 
     if (apply) {
 	interpolate(rowFeSpace, *rhs, values, dirichletNodes);
-	interpolate(colFeSpace, *solution, values, dirichletNodes);
+// 	interpolate(colFeSpace, *solution, values, dirichletNodes);
     }
   }
 

dune/amdis/FileWriter.hpp

+#pragma once
+
+#include <string>
+#include <vector>
+#include <array>
+#include <memory>
+
+#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
+#include <dune/geometry/genericgeometry/referenceelements.hh>
+
+#include "Loops.hpp"
+#include "Initfile.hpp"
+
+namespace AMDiS
+{
+
+template <class Traits>
+class FileWriter
+{
+    using Mesh     = typename Traits::Mesh;
+    using MeshView = typename Mesh::LeafGridView;
+    
+    /// Dimension of the mesh
+    static constexpr int dim = Traits::dim;
+    
+    /// Dimension of the world
+    static constexpr int dow = Traits::dimworld;
+    
+    /// Number of problem components
+    static constexpr int nComponents = Traits::nComponents;
+    
+    
+public:
+    FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& 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, "");
+    }
+
+    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*/);
+    }
+    
+    
+protected:
+    template <class DOFVector, class Vector>
+    void dofVector2vertexVector(DOFVector 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];
+		}
+	    }
+	}
+    }
+
+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;
+};
+
+} // end namespace AMDiS

dune/amdis/Log.hpp

@@ -6,6 +6,23 @@
 #include <iostream>
 #include <sstream>
 
+#ifdef MSG 
+#undef MSG 
+#endif
+#ifdef ERROR 
+#undef ERROR 
+#endif
+#ifdef WARNING 
+#undef WARNING 
+#endif
+#ifdef TEST 
+#undef TEST 
+#endif
+#ifdef TEST_EXIT 
+#undef TEST_EXIT 
+#endif
+
+
 #define AMDIS_STATIC_ASSERT(...) \
   static_assert(__VA_ARGS__, #__VA_ARGS__)
 

dune/amdis/Operator.inc.hpp

@@ -53,8 +53,8 @@ namespace AMDiS
     const int dim = Element::dimension;
     auto geometry = element.geometry();
     
-    const auto& rowFE = rowView.tree().finiteElement();
-    const auto& colFE = colView.tree().finiteElement();
+    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+    const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
     
     int order = getQuadratureDegree(0);
     const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
@@ -71,12 +71,12 @@ namespace AMDiS
 	const double integrationElement = geometry.integrationElement(quadPos);
 	
 	std::vector<Dune::FieldVector<double,1> > rowShapeValues, colShapeValues;
-	rowFE.localBasis().evaluateFunction(quadPos, rowShapeValues);
+	rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
 	
 	if (psiDegree == phiDegree)
 	  colShapeValues = rowShapeValues;
 	else
-	  colFE.localBasis().evaluateFunction(quadPos, colShapeValues);
+	  colLocalBasis.evaluateFunction(quadPos, colShapeValues);
 
 	for (size_t i = 0; i < elementMatrix.N(); ++i) {
 	    for (size_t j = 0; j < elementMatrix.M(); ++j) {
@@ -100,7 +100,7 @@ namespace AMDiS
     const int dim = Element::dimension;
     auto geometry = element.geometry();
     
-    const auto& rowFE = rowView.tree().finiteElement();
+    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
     
     int order = getQuadratureDegree(0);
     const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
@@ -117,7 +117,7 @@ namespace AMDiS
 	const double integrationElement = geometry.integrationElement(quadPos);
 	
 	std::vector<Dune::FieldVector<double,1> > rowShapeValues;
-	rowFE.localBasis().evaluateFunction(quadPos, rowShapeValues);
+	rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
 
 	for (size_t i = 0; i < elementvector.size(); ++i) {
 	    const int local_i = rowView.tree().localIndex(i);
@@ -138,8 +138,8 @@ namespace AMDiS
     const int dim = Element::dimension;
     auto geometry = element.geometry();
     
-    const auto& rowFE = rowView.tree().finiteElement();
-    const auto& colFE = colView.tree().finiteElement();
+    const auto& rowLocalBasis = rowView.tree().finiteElement().localBasis();
+    const auto& colLocalBasis = colView.tree().finiteElement().localBasis();
     
     int order = getQuadratureDegree(2);
     const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
@@ -163,7 +163,7 @@ namespace AMDiS
 
 	// The gradients of the shape functions on the reference element
 	std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients;
-	rowFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
+	rowLocalBasis.evaluateJacobian(quadPos, referenceGradients);
 
 	// Compute the shape function gradients on the real element
 	std::vector<Dune::FieldVector<double,dim> > gradients(referenceGradients.size());

dune/amdis/OperatorTerm.hpp

@@ -151,8 +151,9 @@ namespace AMDiS
     using value_type = typename DOFVectorType::value_type;
     
   public:
-    DOFVectorTerm(DOFVectorType const& dofvector, double factor = 1.0)
-      : vector(dofvector)
+    template <class DOFVectorType_>
+    DOFVectorTerm(DOFVectorType_&& dofvector, double factor = 1.0)
+      : vector(dofvector.getVector())
       , factor(factor)
       , localView(dofvector.getFeSpace().localView())
       , localIndexSet(dofvector.getFeSpace().localIndexSet())
@@ -195,7 +196,7 @@ namespace AMDiS
     }
     
   private:
-    DOFVectorType vector;
+    typename DOFVectorType::BaseVector const& vector;
     double factor;
     
     mutable typename Basis::LocalView localView;
@@ -208,9 +209,9 @@ namespace AMDiS
   
   template <class DOFVectorType>
   DOFVectorTerm<DOFVectorType>
-  valueOf(DOFVectorType const& vector, double factor = 1.0)
+  valueOf(DOFVectorType&& vector, double factor = 1.0)
   {
-      return {vector, factor};
+      return {std::forward<DOFVectorType>(vector), factor};
   }
   
   

dune/amdis/ProblemStat.cpp

@@ -7,6 +7,7 @@ namespace AMDiS
 {
   // explicit template instatiation
 //   template class ProblemStat<ProblemParametersBase<1>>;
-  template class ProblemStat<ProblemParametersBase<2>>;
+  template class ProblemStat<ProblemStatTraits<2,2,1>>;
+  template class ProblemStat<ProblemStatTraits<2,2,2>>;
 //   template class ProblemStat<ProblemParametersBase<3>>;
 } // end namespace AMDiS

dune/amdis/ProblemStat.hpp

@@ -15,6 +15,7 @@
 
 #include "Basic.hpp"
 #include "DirichletBC.hpp"
+#include "FileWriter.hpp"
 #include "Flag.hpp"
 #include "Initfile.hpp"
 #include "Operator.hpp"
@@ -24,13 +25,13 @@
 
 namespace AMDiS {
     
-template <class Param>
+template <class Traits>
 class ProblemStat : public ProblemStatBase
 {
 public:
       
-    using FeSpaces = typename Param::FeSpaces;
-    using Mesh     = typename Param::Mesh;
+    using FeSpaces = typename Traits::FeSpaces;
+    using Mesh     = typename Traits::Mesh;
     using MeshView = typename Mesh::LeafGridView;
     
     using Codim0   = typename MeshView::template Codim<0>;
@@ -41,13 +42,13 @@ public:
     
     
     /// Dimension of the mesh
-    static constexpr int dim = Param::dim;
+    static constexpr int dim = Traits::dim;
     
     /// Dimension of the world
-    static constexpr int dow = Param::dimworld;
+    static constexpr int dow = Traits::dimworld;
     
     /// Number of problem components
-    static constexpr int nComponents = Param::nComponents;
+    static constexpr int nComponents = Traits::nComponents;
     
     
     template <size_t I>
@@ -133,7 +134,7 @@ public:
     
     /// Return the I'th finite element space
     template <size_t I = 0>
-    decltype(auto) getFeSpace(const index_<I> = index_<I>()) const
+    FeSpace<I> const& getFeSpace(const index_<I> = index_<I>()) const
     {
 	return std::get<I>(*feSpaces);
     }
@@ -156,7 +157,7 @@ public:
     template <size_t I = 0>
     auto getSolution(const index_<I> _i = index_<I>())
     {
-	return (*solution)[_i];
+	return std::move((*solution)[_i]);
     }
     
     
@@ -220,6 +221,14 @@ public:
 		  DOFMatrix* matrix,
 		  Vector* rhs);
     
+private:
+    template <size_t I = 0>
+    typename FeSpace<I>::NodeFactory& 
+    getNodeFactory(const index_<I> = index_<I>())
+    {
+	return const_cast<typename FeSpace<I>::NodeFactory&>(std::get<I>(*feSpaces).nodeFactory());
+    }
+    
 private:
     /// Name of this problem.
     std::string name;
@@ -242,6 +251,8 @@ private:
     /// FE spaces of this problem.
     std::shared_ptr<FeSpaces> feSpaces; // eventuell const
     
+    std::shared_ptr<FileWriter<Traits>> filewriter;
+    
     SystemMatrix systemMatrix;
     
     std::shared_ptr<SystemVectorType> solution;
@@ -254,9 +265,10 @@ private:
   
   
 #ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
-  extern template class ProblemStat<ProblemParametersBase<1>>;
-  extern template class ProblemStat<ProblemParametersBase<2>>;
-  extern template class ProblemStat<ProblemParametersBase<3>>;
+//   extern template class ProblemStat<ProblemStatTraits<1>>;
+  extern template class ProblemStat<ProblemStatTraits<2,2,1>>;
+  extern template class ProblemStat<ProblemStatTraits<2,2,2>>;
+//   extern template class ProblemStat<ProblemStatTraits<3>>;
 #endif
   
 }

dune/amdis/ProblemStat.inc.hpp

@@ -4,6 +4,7 @@
 #include <dune/istl/solvers.hh>
 #include <dune/istl/preconditioners.hh>
 #include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <dune/istl/matrixmarket.hh>
 
 #include "AdaptInfo.hpp"
 #include "Loops.hpp"
@@ -23,14 +24,16 @@ 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();
 	    }
 	    
@@ -41,6 +44,9 @@ namespace AMDiS
 	    AMDIS_WARNING("no mesh created");